Определение доверительных интервалов для оценки максимального правдоподобия

Похоже, ты в порядке с первой частью, поэтому я займусь вторым и третьим очками.

Существует множество способов сочетания гладких кривых с scipy.interpolate и сплайнами или с scipy.optimize.curve_fit . Лично я предпочитаю curve_fit , потому что вы можете предоставить свою собственную функцию и позволить ей соответствовать параметрам для вас.

В качестве альтернативы, если вы не хотите изучать параметрическую функцию, вы можете сделать простой сглаживание с помощью numpy.convolve .

Что касается качества кода: вы не пользуетесь скоростью numpy, потому что вы делаете что-то в чистом питоне. Я бы написал ваш (существующий) код следующим образом:

 from __future__ import division import numpy as np import matplotlib.pyplot as plt # N is the true number of books. # t is the number of weeks. # unk is the true number of repeats found t = 30 unk = 3 def numberrepeats(N, t, iters): rand = np.random.randint(0, N, size=(t, iters)) return t - np.array([len(set(r)) for r in rand]) iters = 1000 ydata = np.empty(500-10) for N in xrange(10,500): sampledunk = np.count_nonzero(numberrepeats(N,t,iters) == unk) ydata[N-10] = sampledunk/iters print "MLE is", np.argmax(ydata) xdata = range(10, 500) print len(xdata), len(ydata) plt.plot(xdata,ydata) plt.show() 

Вероятно, это возможно еще больше оптимизировать, но это изменение приводит к тому, что время работы вашего кода составляет от ~ 30 секунд до ~ 2 секунд на моей машине.

Простой (числовой) способ получить доверительный интервал – просто запустить сценарий много раз и посмотреть, насколько сильно ваша оценка меняется. Вы можете использовать это стандартное отклонение для вычисления доверительного интервала.

В интересах времени другой вариант состоит в том, чтобы запустить кучу испытаний при каждом значении N (я использовал 2000), а затем использовать случайную подвыборку этих испытаний, чтобы получить оценку стандартного отклонения оценки. В основном это предполагает выбор поднабора проб, генерирование кривой вероятности с использованием этого подмножества, а затем поиск максимума этой кривой для получения вашей оценки. Вы делаете это во многих подмножествах, и это дает вам кучу оценок, которые вы можете использовать, чтобы найти доверительный интервал для вашей оценки. Мой полный сценарий выглядит следующим образом:

 import numpy as np t = 30 k = 3 def trial(N): return t - len(np.unique(np.random.randint(0, N, size=t))) def trials(N, n_trials): return np.asarray([trial(N) for i in xrange(n_trials)]) n_trials = 2000 Ns = np.arange(1, 501) results = np.asarray([trials(N, n_trials=n_trials) for N in Ns]) def likelihood(results): L = (results == 3).mean(-1) # boxcar filtering n = 10 L = np.convolve(L, np.ones(n) / float(n), mode='same') return L def max_likelihood_estimate(Ns, results): i = np.argmax(likelihood(results)) return Ns[i] def max_likelihood(Ns, results): # calculate mean from all trials mean = max_likelihood_estimate(Ns, results) # randomly subsample results to estimate std n_samples = 100 sample_frac = 0.25 estimates = np.zeros(n_samples) for i in xrange(n_samples): mask = np.random.uniform(size=results.shape[1]) < sample_frac estimates[i] = max_likelihood_estimate(Ns, results[:,mask]) std = estimates.std() sterr = std * np.sqrt(sample_frac) # is this mathematically sound? ci = (mean - 1.96*sterr, mean + 1.96*sterr) return mean, std, sterr, ci mean, std, sterr, ci = max_likelihood(Ns, results) print "Max likelihood estimate: ", mean print "Max likelihood 95% ci: ", ci 

Для этого метода есть два недостатка. Во-первых, поскольку вы принимаете много подвыборки из одного и того же набора проб, ваши оценки не являются независимыми. Чтобы ограничить эффект этого, я использовал только 25% результатов для каждого подмножества. Еще один недостаток заключается в том, что каждая подвыборка является лишь частью ваших данных, поэтому оценки, полученные из этих подмножеств, будут иметь больше дисперсии, чем оценки, полученные из полного сценария много раз. Чтобы учесть это, я вычислил стандартную ошибку как стандартное отклонение, деленное на квадратный корень из 4, так как у меня было в четыре раза больше данных в моем полном наборе данных, чем в одной из подвыборки. Тем не менее, я недостаточно разбираюсь в теории Монте-Карло, чтобы узнать, математически ли это звучит. Выполнение моего скрипта несколько раз показалось, что мои результаты были разумными.

Наконец, я использовал фильтр вагонов на кривых правдоподобия, чтобы немного сгладить их. В идеале это должно улучшить результаты, но даже при фильтрации все еще значительная изменчивость результатов. При вычислении значения для общей оценки я не был уверен, будет ли лучше вычислять одну кривую правдоподобия из всех результатов и использовать максимум этого (это то, что я закончил), или использовать среднее значение для всех оценки подмножества. Использование среднего для подмножеств оценок может помочь сгладить некоторую шероховатость кривых, оставшихся после фильтрации, но я не уверен в этом.

Вот ответ на ваш первый вопрос и указатель на решение для второго:

 plot(xdata,ydata) # calculate the cumulative distribution function cdf = np.cumsum(ydata)/sum(ydata) # get the left and right boundary of the interval that contains 95% of the probability mass right=argmax(cdf>0.975) left=argmax(cdf>0.025) # indicate confidence interval with vertical lines vlines(xdata[left], 0, ydata[left]) vlines(xdata[right], 0, ydata[right]) # hatch confidence interval fill_between(xdata[left:right], ydata[left:right], facecolor='blue', alpha=0.5) 

Это приводит к следующему рисунку:

Я постараюсь ответить на вопрос 3, когда у меня будет больше времени 🙂