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

Я пытаюсь написать код для создания доверительных интервалов для количества разных книг в библиотеке (а также создать информативный сюжет).

Мой двоюродный брат учится в начальной школе, и каждую неделю ему дают книгу своего учителя. Затем он читает его и возвращает его вовремя, чтобы получить еще одну на следующей неделе. Через некоторое время мы начали замечать, что он получал книги, которые он читал раньше, и со временем это стало более распространенным.

Скажем, истинное количество книг в библиотеке равно N, и учитель выбирает один случайным образом (с заменой), чтобы дать вам каждую неделю. Если в течение недели t число случаев, когда вы получили прочитанную книгу, равно x, я могу получить оценку максимального правдоподобия для количества книг в библиотеке, следующих за https://math.stackexchange.com/questions/ 615464 / how-many-books-are-in-an-library .


Пример: рассмотрим библиотеку с пятью книгами A, B, C, D и E. Если вы получите книги [A, B, A, C, B, B, D] в семь последовательных недель, то значение для x ( количество дубликатов) будет [0, 0, 1, 1, 2, 3, 3] после каждой из этих недель, то есть через семь недель вы получили книгу, которую вы уже читали три раза.


Чтобы визуализировать функцию правдоподобия (предполагая, что я понял, что правильно), я написал следующий код, который, я считаю, отображает функцию правдоподобия. Максимум составляет около 135, что действительно является максимальной оценкой правдоподобия в соответствии с вышеприведенной ссылкой MSE.

from __future__ import division import random import matplotlib.pyplot as plt import numpy as np #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): return t - len(set([random.randint(0,N) for i in xrange(t)])) iters = 1000 ydata = [] for N in xrange(10,500): sampledunk = [numberrepeats(N,t) for i in xrange(iters)].count(unk) ydata.append(sampledunk/iters) print "MLE is", np.argmax(ydata) xdata = range(10, 500) print len(xdata), len(ydata) plt.plot(xdata,ydata) plt.show() 

Результат выглядит

введите описание изображения здесь

Мои вопросы таковы:

  • Есть ли простой способ получить доверительный интервал 95% и построить его на диаграмме?
  • Как вы можете наложить сглаженную кривую на сюжет?
  • Есть ли лучший способ, которым мой код должен был быть написан? Он не очень изящный и довольно медленный.

Поиск 95% -ного доверительного интервала означает поиск диапазона оси х, так что в 95% случаев эмпирическая оценка максимального правдоподобия, которую мы получаем по выборке (которая теоретически должна составлять 135 в этом примере), попадет в нее. Ответ @mbatchkarov дал в настоящее время не так правильно.


Теперь есть математический ответ на странице https://math.stackexchange.com/questions/656101/how-to-find-a-confidence-interval-for-a-maximum-likelihood-estimate .

3 Solutions collect form web for “Определение доверительных интервалов для оценки максимального правдоподобия”

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

Существует множество способов сочетания гладких кривых с 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, когда у меня будет больше времени 🙂

  • Как рассчитать вероятность скручивания кривой в scipy?
  • Любая библиотека Python создает таблицы регрессии стиля публикации
  • Разница в статистических моделях Python OLS и R's lm
  • ImportError: Ошибка загрузки DLL: при импорте statsmodels
  • Statsmodels Python - взвешенный GLM
  • Пользовательские приоритеты в PyMC
  • Statsmodels: Рассчитать установленные значения и R-квадрат
  • ANOVA в python с использованием рамки данных pandas с statsmodels или scipy?
  • Python - лучший язык программирования в мире.