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

Ответ 1

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

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

Ответ 2

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

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

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% результатов для каждого подмножества. Другим недостатком является то, что каждая подвыборка является лишь частью ваших данных, поэтому оценки, полученные из этих подмножеств, будут иметь больше дисперсии, чем оценки, полученные из многократного запуска полного script. Чтобы учесть это, я вычислил стандартную ошибку как стандартное отклонение, деленное на квадратный корень из 4, так как у меня было в четыре раза больше данных в моем полном наборе данных, чем в одной из подвыборки. Тем не менее, я недостаточно разбираюсь в теории Монте-Карло, чтобы узнать, математически ли это звучит. Выполнение моего script несколько раз показало, что мои результаты были разумными.

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

Ответ 3

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

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, когда у меня будет больше времени:)