13

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

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

假设图书馆的真实书籍数量为 N，老师每周随机（有替换）随机挑选一本给你。如果在第 t 周，您收到一本书的次数是 x，那么我可以按照https://math.stackexchange.com/questions/对图书馆中的图书数量进行最大似然估计615464/图书馆中有多少书。

示例：考虑一个有五本书 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有一个数学答案。

4

3 回答 3

8

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

有很多方法可以拟合平滑曲线，使用scipy.interpolate和 splines，或者使用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 秒。

于 2014-02-01T16:56:31.613 回答
6

获得置信区间的一种简单（数字）方法是简单地运行脚本多次，并查看您的估计值有多少变化。您可以使用该标准差来计算置信区间。

为了节省时间，另一种选择是在每个 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 滤波器来平滑它们。理想情况下，这应该会改善结果，但即使经过过滤，结果仍然存在相当大的可变性。在计算总体估计量的值时，我不确定是否会更好地从所有结果中计算一条似然曲线并使用其中的最大值（这就是我最终要做的），还是使用所有结果的平均值子集估计器。使用子集估计器的平均值可能有助于消除过滤后剩余的曲线中的一些粗糙度，但我不确定这一点。

于 2014-02-01T19:37:30.710 回答
5

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

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

于 2014-01-30T11:04:46.600 回答