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

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上有一个数学答案。

- Simd

---

你应该在“numberrepeats”中设置一个种子，以便每个N使用相同的随机样本。 - Josef

http://python-for-signal-processing.blogspot.co.uk/2012/10/maximum-likelihood-estimation-maximum.html 看起来很相关，尽管我还没有完全理解它。 - Simd

你的问题描述与你在脚本中实际实现的问题非常不同。只有通过访问你在math.stackexchange.com上的问题链接，我才能找出你真正的意思。你应该考虑重写你的问题，以反映math.stackexchange.com上的评论讨论。短语“如果在第t周你收到了一本你之前读过x次的书”让我觉得你必须收到相同的书'x'次，但显然情况并非如此。 - hunse

@hunse 谢谢。增加了澄清。这样更清楚吗？ - Simd

是的，现在这个问题中的短语更加清晰了。我也为了清晰度修改了例子。 - hunse

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

- perimosocordiae

---

谢谢。虽然第一部分没有回答。但是当前的答案不正确。我在问题中有一条评论。 - Simd

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滤波器来使它们更加平滑。理想情况下，这应该改善结果，但即使进行了过滤，结果仍然存在相当大的变异性。在计算总体估计器的值时，我不确定是计算所有结果的一个可能性曲线并使用其中的最大值（这就是我最终做的），还是使用所有子集估计器的平均值。使用子集估计器的平均值可能有助于消除过滤后剩余曲线中的一些粗糙度，但我不确定。

- hunse

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

这将产生以下图像： 。 我会在有更多时间时尝试回答第三个问题 :)。

- mbatchkarov

---

我使用 np.argmax、np.sum 和修正打字错误 vlines(xdata[right], 0, ydata[right]) 成功运行了你的代码。 - Simd

是的，我在ipython中运行它，它会自动导入numpy函数。抱歉 :) - mbatchkarov

啊...这是似然函数置信区间的错误公式。我相信正确的方法是取对数，找到最大值，然后向左右步进2个单位。我已经在问题中添加了一些内容。或者，通过模拟来完成。 - Simd