8.9 火车头问题

火车头问题是一个非常经典的估计问题,又叫“德国坦克问题”。下面是Mosteller在Fifty Challenging Problems in Probability中提到的版本:

铁路公司将它所有的火车头都进行了编号,从1到N。有一天你看见一个编号为60的火车头,那该铁路公司总共有多少火车头呢?

在接着往下讨论之前,我们先想想如下问题。

  1. 对于一个给定的估计N^\hat{N},观测的似然函数是什么?极大似然估计量是什么?
  2. 假如我们看到编号为i的火车,我们有理由猜测一个乘数a,并用N^=ai\hat{N}=ai来估计总的火车数。那么我们该怎么样选择a使得估计的结果能使均方误差最小?
  3. 仍然假设N^=ai\hat{N}=ai, 我们能找到一个使得N^\hat{N}为无偏估计量的a吗?
  4. N是多少时60会是观测到的平均值?
  5. 在假设先验概率为1到200上的离散均匀分布的条件下,贝叶斯后验分布是什么样子的?

对于一个给定的估计量N^\hat{N},观测到编号为i(i≤ N^\hat{N})的火车的概率为1/N^\hat{N},i>N^\hat{N}的概率为0。所以N的极大似然估计量是 N^\hat{N}=i。换言之,如果观测到的火车的编号是60,而且我们要以最大的概率保证结果的正确性,那么我们就会猜测铁路公司有60辆火车。

但是极大似然估计的结果从均方误差的角度来看并不理想。我们可以通过选择一个N^\hat{N}=ai使得估计量的均方误差达到最小,这里唯一需要做的是估计好a的值。

假设实际会有N辆火车。我们看到编号为i的火车后就猜测铁路公司有ai辆火车,那么平方误差为(aiN)2(a^*i^*-N)^2

假设我们观测了N次,且看到了所有的火车,那么均方差就是: MSE=1Ni=1N(aiN)2 MSE = \frac{1}{N}\sum_{i=1}^{N}(ai-N)^2 为了使均方差最小,我们对a进行求导: dMSEda=1Ni=1N2i(aiN)=0 \frac{dMSE}{da} = \frac{1}{N}\sum_{i=1}^{N}2i(ai-N)=0 解得: a=3N2N+1 a = \frac{3N}{2N+1} 粗略看来这个结果似乎没什么用,因为等式的右边有N。我们要想知道a就必须知道N,但是如果我们知道了N,上面的估计就全然不需要了。

考虑到N很大的时候a收敛于3/2,所以这里选择对N的估计量为N^\hat{N}=3i/2。

再从无偏估计的角度出发,首先计算平均误差: ME=1Ni=1N(aiN) ME = \frac{1}{N}\sum_{i=1}^{N}(ai-N) 令ME=0,得到 a=2NN1 a= \frac{2N}{N-1} 当N很大时,a收敛到2,所以这里我们又可以选择N^=2i\hat{N}=2i

到目前为止,我们已经构造了3个估计量:i、3i/2和2i,分别满足极大似然、最小均方误差和无偏性。

我们还有另外一种估计的方式:选择满足总体均值等于样本均值的N^\hat{N}作为N的估计量。假设我们观测到一辆火车的编号为i,样本均值恰好也为i,这时满足群体的均值N^\hat{N}等于样本均值的火车数的估计量为N^=2i1\hat{N}=2i-1

最后,我们从贝叶斯统计的角度来回答这个问题。我们计算 P(Hni)=P(iHn)P(Hn)P(i) P(H_n|i) = \frac{P(i|H_n)P(H_n)}{P(i)} 这里Hn是一个假设,假设总共有n辆火车。i表示我们的观测,即观测到的火车的编号。当i<n时, P(iHn)=1/nP(i|H_n)=1/n,其余情况为0,P(i)是归一化常数。

如果N的先验分布是1到200上的一个均匀分布,我们就遍历这200种假设,并计算每种假设的后验概率。读者可以从http://thinkstats.com/locomotive.py下载到实现所用的代码。最终的结果如图8-1所示。

enter image description here

图8-1 火车数后验分布

根据贝叶斯后验概率,我们得到90%的可信区间为[63, 189],这仍然是一个非常大的范围。仅仅观测到一辆火车的编号并不能为任意假设提供非常强的证据,虽然它将n<i的可能性排除在外了。

如果一开始我们设定不同的先验概率,后验概率就会显著不同,这能帮助你理解为什么其他估计量有多个。

本节我们针对同一个参数构造了多个不同的估计量,我们可以认为这些是从一些不够精确的先验出发得到的结果。假如有足够的先验信息,那么这些估计量将倾向于收敛到同一个值。总之,在我们的例子里,没有一个统计量能够满足我们期望的所有性质。

习题8-6

推广locomotive.py,使之能处理当观测到的火车数多于1辆时的情况,读者只需改动几行代码就可以完成。

看看你能否回答上述情况下的其他问题。维基百科页面http://wikipedia.org/wiki/German_tank_problem提供了更多的问题和讨论。