绘制最大似然估计的置信区间

Question

我正在尝试编写代码来生成图书馆中不同书籍数量的置信区间（以及生成信息图）。

我表弟上小学，他的老师每周都会给我一本书。然后，他阅读并及时归还，以便下周再获得一份。过了一段时间，我们开始注意到他收到了他以前读过的书，随着时间的推移，这种情况逐渐变得越来越普遍。

假设图书馆的真实图书数量为 N，老师每周随机（有替换）统一挑选一本给你。如果在第 t 周，你收到一本你读过的书的次数是 x，那么我可以对 https://math.stackexchange.com/questions/615464/how-many-books-are-in-a-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% 置信区间意味着找到 x 轴的范围，以便我们通过抽样得到的经验最大似然估计值（在本例中理论上应该为 135）有 95% 的时间落在该范围内。 @mbatchkarov 给出的答案目前没有正确执行此操作。

---

现在https://math.stackexchange.com/questions/656101/how-to-find-a-confidence-interval-for-a-maximum-likelihood-estimate 有一个数学答案。

Answer 1

看起来你在第一部分没问题，所以我会处理你的第二和第三点。

有很多方法可以拟合平滑曲线，使用scipy.interpolate 和样条曲线，或者使用scipy.optimize.curve_fit。就个人而言，我更喜欢curve_fit，因为您可以提供自己的函数并让它适合您的参数。

或者，如果您不想学习参数函数，可以使用 numpy.convolve 进行简单的滚动窗口平滑。

至于代码质量：您没有利用 numpy 的速度，因为您是在纯 python 中做事。我会像这样编写您的（现有）代码：

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 秒。

Answer 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 的平方根，因为我的完整数据集中的数据是其中一个子样本中的四倍。但是，我对蒙特卡洛理论还不够熟悉，无法知道这在数学上是否合理。多次运行我的脚本似乎表明我的结果是合理的。

最后，我确实在似然曲线上使用了 boxcar 过滤器来稍微平滑它们。理想情况下，这应该会改善结果，但即使进行了过滤，结果仍然存在相当大的可变性。在计算总体估计量的值时，我不确定是否会更好地从所有结果中计算一条似然曲线并使用其中的最大值（这就是我最终要做的），还是使用所有结果的平均值子集估计器。使用子集估计器的平均值可能有助于消除过滤后剩余的曲线中的一些粗糙度，但我不确定这一点。

Answer 3

这是您第一个问题的答案和第二个问题的解决方案pointer：

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 :)