德国坦克问题

在统计学理论的估计中，用不放回抽样来估计离散型均匀分布最大值问题中著名的**德国坦克问题**（），它因在第二次世界大战中用于估计德国坦克数量而得名。

这些分析说明了频率推断和贝叶斯推断之间的不同。

基于“单个”样本估计的样本总数各有不同，而在“多个”样本的基础上估计则是现实生活中一个很有意义的估计问题，它的答案很简单，但并不那么明显。

## 范例

假设某个情报人员已经发现了辆坦克，其序列号分别为2、6、7、14，观测到的最大的序列号为。坦克未知的总数设为*N*。

按“频率论”方法的估算公式，计算出坦克总数为：
:N \approx m + \frac{m}{k} - 1 = 16.5

而“贝叶斯”分析会给出一个估计坦克数量的概率质量函数，如下所示：
:\Pr(N=n) = \begin{cases}
   0 &\text{if } n 

从中，我们可以估计坦克数量：

:\begin{align}
       N &\approx \mu \pm \sigma = 19.5 \pm 10 \\
     \mu &= (m - 1)\frac{k - 1}{k - 2} \\
  \sigma &= \sqrt{\frac{(k-1)(m-1)(m-k+1)}{(k-3)(k-2)^2}}
\end{align}

这个分布的偏度为正，因为其中有一个事实，即至少有14辆坦克。

## 历史上的问题

在战争的过程中，西方盟国一直在努力确定德国的生产速度，并以两种主要方式得到这个数据：常规情报收集和统计估计。大部分情况下，统计估计的可信度很大程度上决定于常规情报收集。而有时常规情报收集会和统计估计一起结合使用，如在D日前夕估计豹式坦克的生产速度。

盟军指挥机构已经注意到了出现在意大利的五号坦克（豹式坦克），速度快，装备有75 mm/L70长管戰車砲，是很不寻常的重型坦克，但在法国北部编号就很小，就和虎I坦克在突尼斯那样。由于谢尔曼坦克在对阵三号坦克和四号坦克时表现良好，美国陆军对它非常自信，以至于迫不及待地希望与他们相遇。就在D日的不久前，有消息称德国手中有大量的五号坦克。

为了确定这种说法的真实性，盟军要尝试估计正在生产坦克的数量。为此，他们要使用被缴获或被摧毁的坦克的序列号。他们主要使用变速箱的号码，因为其为两个完整的序列。底盘和发动机号码也有使用，但它们的编号规则比较复杂。各种其他的零件用来对分析结果进行交叉检查。对轮胎也做了类似的分析，其上观察到了顺序编号（即1、2、3、......、*N*)。

对坦克车轮的分析产生了对使用中的车轮模具数量的估计。在与英国车轮制造商讨论过后，他们估计了这么多的模具可以生产多少车轮，进而是每个月可生产的坦克数量。对两辆坦克（每辆32个车轮，总计64个车轮）车轮的分析的结果是1944年2月的生产数量估计在270左右，大大超出此前预期。

德国战后公布的记录显示，1944年2月一个月的生产量是276辆。统计方法结果的精确度是常规情报收集方法所远远不能达到的，而“德国坦克问题”这个词也成为了这种统计分析问题的标志。

估计产量并不是这种序列号分析的唯一应用。它也用于探查更多德国生产的信息，包括工厂数目、工厂的相对重要性、供应链长度（基于生产和使用之间的滞后程度）、生产工艺的改变、及对诸如橡胶等资源的使用。

### 具体数据
根据常规盟军情报的估计，德国在1940年6月和1942年9月之间，每月大约能生产1,400辆坦克。将缴获坦克的序列号代入下文的公式，可计算出每月246辆。战后，从阿尔伯特·斯佩尔所管辖的部门缴获的德国生产记录显示，实际数目是245辆。

某些特定月份的估计如下：

### 类似分析

在二战期间，类似的序列号分析也在其他军事装备上应用过，其中最成功的属V-2火箭。

二战期间，德国情报部门分析了苏联军事装备工厂，朝鲜战争期间苏联装备工厂也被分析过。苏联在二战期间也估计过德国坦克的生产数。

20世纪80年代，一些美国人被获准进入以色列的梅卡瓦坦克的生产线。生产编号是分类过的，但水箱上有序列号，可以用来估计生产量。

该公式在非军事中也有使用，如估计Commodore 64计算机的总数，其结果（1.25亿）与官方数字相当匹配。

### 应对措施

为了防止序列号分析，可以不使用序列号，或者减少可用的辅助信息。此外可以使用抗密码分析的序列号，最有效的方法是从比生产物品数量大得多的列表中随机选择一个号码，之后就不再使用（一次性密码本），或产生随机数后检查是否在已分配编号列表上；编号可能会发生碰撞，除非编号可能的位数数目在生产物品数字位数的两倍以上（其中序列号可以以任何数为底）；参见生日问题。对此可使用密码学安全伪随机数生成器。要想从序列号中反推生产订单，对所有这些方法都需要一个查找表（或破解加密器），其使序列号变得复杂：如对于加密器，不能一次反推一串序列号，每个号码都必须单独查找，或者是一张生成列表。

另外，还可以对连续的序列号进行加密，同时方便解码，但这会带来潜在的已知明文攻击：即使从任意数开始，明文都会具有特定模式（即数字是有序的）。一个例子是在肯·福勒特的小说《代码为零》中，木星-C火箭的序列号加密规则为：



此处的代码字为Huntsville（省略重复字母），用以获得一张10字的键值表。因此，13号火箭编为“HN”，而24号火箭则是“UT”。

## 频率论分析

### 最小方差无偏估计
对于点估计（估算出单个总体（\hat{N}）值），最小方差无偏估计（MVUE，或UMVU估计）由下式给出：

:\hat{N} = m\left(1 + k^{-1}\right) - 1

其中*m*是所观察到的最大序号（样本最大值），而*k*是观察到的坦克数目（样本容量）。注意，一旦观察到一个序列号，它就不再在样本池中，也不会被再次观察到。

其方差为

: \operatorname{var}(\hat{N}) = \frac{1}{k}\frac{(N-k)(N+1)}{(k+2)} \approx \frac{N^2}{k^2} \text{ for small samples } k \ll N

因而标准差约为*N*/*k*，即样本间距的（总体）平均大小；注意与前文中的*m*/*k*相比。

#### 直观
公式可以直观地理解为样本最大值加上样本中观测值之间的平均间距，由于是最大似然估计，样本最大值被用作初始估计值，再加上间距以补偿样本最大值的负偏差，以此作为总体最大值的一个估计，因而可写成：

:\hat{N} = m + \frac{m - k}{k}= m + mk^{-1} - 1 = m\left(1 + k^{-1}\right) - 1

可以想象样品在整个区间内均匀分布，而更多的样本就在0至*N* + 1的区间之外。
如果在0和编号最小样本（样本最小值）之间选一个初始间距，那么样本间的平均间距是(m - k)/k；有-k是因为样本本身在计算样本间距时不应算入。

这一理念确立并推广了最大间距估计的方法。

#### 推导
样本最大值等于*m*的概率为\tbinom{m - 1}{k - 1}\big/\tbinom Nk，其中\tbinom \cdot\cdot是二项式系数。

样本最大值的期望值为
:\begin{align}
            \mu &= \sum_{m=k}^N m\frac{\tbinom{m - 1}{k - 1}}{\tbinom Nk} = \frac{k(N + 1)}{k + 1} \\
  \Rightarrow N &= \mu\left(1 + k^{-1}\right) - 1
\end{align}

因而
:\begin{align}
  \mu\left(1 + k^{-1}\right) - 1 &= E\left[m\left(1 + k^{-1}\right) - 1\right] \\
             \Rightarrow \hat{N} &= m\left(1 + k^{-1}\right) - 1
\end{align}

为*N*的无偏估计。

为了表明这是UMVU估计：
* 首先表明样本最大值是总体最大值的充分统计，使用类似于在充分性：均匀分布中具体阐述的方法（但对于德国坦克问题，我们必须排除在样本中出现两次序列号的结果）;
* 接着表明这是一个完整统计。 
* 那么根据莱曼–薛费定理，如上纠正过偏差的样本最大值是无偏的，也就是UMVU估计。

### 置信区间
此外，或者补充来说，对于“点”估计可以进行“区间”估计，如置信区间。基于*k*个样本落在覆盖全范围*p*（0 ≤ *p* ≤ 1）的区间的概率为*p**k*，可以很容易地计算出这些结果（在此段中假设取样后“放回”以简化计算，如果取样后不放回则会夸大概率，区间也将过于保守）。

因此，样本最大值的分位数的抽样分布曲线是*x*1/*k*在0到1上的图像：样本最大值*m*的第*p*到第*q*个分位数是区间[*p*1/*k**N*, *q*1/*k**N*]。反转这个结果就产生总体最大值的相应置信区间[*m*/*q*1/*k*, *m*/*p*1/*k*]。下界非常接近*m*，所以取*p* = 5%到100%得到的非对称置信区间更有意义；*k* = 5时会产生\scriptstyle 0.05^{1/5} \;\approx\; 0.55，所以区间为[*m*, 1.82*m*]。

例如，对*k* = 5，取对称的95%区间，即令*p* = 2.5%，*q* = 97.5%，可得\scriptstyle 0.025^{1/5} \;\approx\; 0.48,\, \scriptstyle 0.975^{1/5} \;\approx\; 0.995，所以置信区间大约为\scriptstyle\left[1.005m,\, 2.08m\right]。

更一般地，（下偏）95％置信区间为\scriptstyle\left[m,\, m/0.05^{1/k}\right] \;=\; \left[m,\, m \cdot 20^{1/k}\right]。对于一系列的*k*，附上UMVU的点估计值（加1以便于阅读）以便参考，可得：

从中可发现：
* 对于较小的样本容量，置信区间很宽，反映了其估计值具有很大的不确定性。 
* 范围迅速缩小，反映了“所有”样品远低于最大值的概率呈指数衰减。
* 置信区间呈正偏态，因为*N*永远不可能低于样本最大值，但却可能比它任意大。

注意，*m*/*k*不能天真地用作（或者更确切地说，是(*m* + *m*/*k* − 1)/*k*）标准差*SE*的估计值，因为标准差的估计值基于“总体”最大值（一个参数），在这里用估计值来估计误差会造成循环论证。

在某些领域，尤其是在未来学中，以这种方式基于单个样本估计置信区间——考虑到它是一个随机取样的分位数（根据平庸原理）——称为哥白尼原则。这一点特别适用于基于目前的年龄来估计寿命，尤其是末日论将其用于估计人类的预期存活时间。

## 贝叶斯分析
在德国坦克问题中，贝叶斯方法要考虑当观察到的坦克总数\scriptstyle K等于数\scriptstyle k、序列号最大值\scriptstyle M等于数\scriptstyle m时敌方坦克总数\scriptstyle N等于数\scriptstyle n的可信度\scriptstyle (N=n\mid M=m, K=k)

为了简便起见，以下将\scriptstyle (N=n\mid M=m,K=k)写成\scriptstyle (n\mid m,k)。

条件概率的法则会给出
:(n\mid m,k) = (m\mid n,k)\frac {(n\mid k)}{(m\mid k)}

表达式\scriptstyle (m\mid n,k)=(M=m\mid N=n,K=k)是当敌方坦克总数等于已知的\scriptstyle n、且观察到了\scriptstyle k辆敌方坦克时，所观察到的序列号最大值等于\scriptstyle m的条件概率。其为
:
  (m\mid n,k) =
  \begin{cases}
    \frac{\binom{m - 1}{k - 1}}{\binom{n}{k}} &\text{if } k \le m \le n\\
    0                                         &\text{otherwise}
  \end{cases}


其中二项式系数\scriptstyle \binom n k是从总数为\scriptstyle n的总体中取\scriptstyle k个样本的方法数。

表达式\scriptstyle (m\mid k)=(M=m\mid K=k)是在真正观测前，假设观察了*k*辆坦克后，最大序号等于*m*的概率。通过合并所有可能的\scriptstyle n，\scriptstyle (m\mid k)可以改写为其它形式。
:\begin{align}
  (m\mid k)
    &=(m\mid k)\cdot 1 \\
    &=(m\mid k){\sum_{n=0}^\infty(n\mid m,k)} \\
    &=(m\mid k){\sum_{n=0}^\infty(m\mid n,k)\frac {(n\mid k)}{(m\mid k)}} \\
    &=\sum_{n=0}^\infty(m\mid n,k)(n\mid k)
\end{align}

表达式\scriptstyle (n\mid k)=(N=n\mid K=k)为当已观察到*k*辆坦克、但序列号尚未统计前坦克总数等于*n*的可信度。假定它是某种离散均匀分布
:
  (n\mid k) = 
  \begin{cases}
    \frac 1{\Omega - k} &\text{if } k \le n 

上限\Omega必须是有限的，因为该函数
:
  f(n)=\lim_{\Omega\rarr\infty}
  \begin{cases}
    \frac 1{\Omega - k} &\text{if } k \le n 

的结果为:f(n) = 0，而这不是一个概率质量函数。

因而
:
  (n\mid m,k) = 
  \begin{cases}
    \frac{(m|n,k)}{\sum_{n=m}^{\Omega - 1} (m|n,k)} &\text{if } m \le n 

如果\scriptstyle \sum_{n=m}^\infty(m|n,k)，那么不受欢迎的变量\scriptstyle \Omega就能从表达式中消失。

:
  (n\mid m,k) =
  \begin{cases}
    0                                          &\text{if } n 

当*k* ≥ 1时，敌方坦克数量分布的众数为*m*。

当*k* ≥ 2时，敌方坦克数量“等于”n的可信度为
:
  (N=n|M=m\ge k,K=k\ge 2) =
  \begin{cases}
    0                                                      &\text{if } n 

而敌方坦克数量\scriptstyle N“大于”\scriptstyle n的可信度为
:
  (N > n|M = m \ge k , K = k \ge 2) = 
  \begin{cases}
    1                                              &\text{if } n 

当*k* ≥ 3时，N的均值有限：
:\frac{(m - 1)(k - 1)}{k - 2}

当*k* ≥ 4时，\scriptstyle N的标准差有限：
:\sqrt{\frac{(m - 1)(k - 1)(m + 1 - k)}{(k - 2)^2(k - 3)}}

这些公式将在下面推导。

### 求和公式
下面的等式用二项式系数表示法来简化德国坦克问题中的级数。
:\sum_{n=m}^\infty \frac 1 {\binom n k} = \frac k{k - 1}\frac 1 {\binom{m - 1}{k - 1}}

这个求和公式有点类似于积分公式
:\int_{n=m}^\infty \frac {dn}{n^k} = \frac 1{k - 1}\frac 1{m^{k - 1}}

这些公式在*k* > 1时适用。

### 一辆坦克
从*n*辆坦克的总体中随机观察一辆坦克，当*m* ≤ *n*时，其序列号为*m*的概率为1/*n*，而当*m* > *n*时概率是零。用艾弗森括号表示法可写成
:(M=m|N=n,K=1) = (m|n) = \frac{[m \le n]}{n}

这是\scriptstyle m的条件概率质量分布函数。

当*m*为定值时，这是一个*n*的似然函数。
:\mathcal{L}(n) = \frac{[n\ge m]}{n}

最大似然估计的坦克总数为*N*0 = *m*。

总概率为无穷大，因为尾部为一个调和级数数列。

:\sum_n \mathcal{L}(n) = \sum_{n=m}^\infty \frac{1}{n} = \infty

但
:\begin{align}
  \sum_n \mathcal{L}(n)[n 

其中H_n为调和数。

可信度质量分布函数依赖于先前的限制\scriptstyle \Omega：
:\begin{align}
       &(N=n\mid M=m,K=1) \\
  = {} &(n|m) = \frac{[m\le n]}{n} \frac{[n

\scriptstyle N的均值为
:\begin{align}
  \sum_n n\cdot(n|m) &= \sum_{n=m}^{\Omega - 1} \frac{1}{H_{\Omega - 1} - H_{m - 1}} \\
                     &= \frac{\Omega - m}{H_{\Omega - 1} - H_{m - 1}} \\
                     &\approx \frac{\Omega - m}{\log\left(\frac{\Omega - 1}{m - 1}\right)}
\end{align}

### 两辆坦克
如果观察到了两辆坦克，而不是一辆，那么所观察到的两个序列号中较大值为*m*的概率为

:(M=m\mid N=n,K=2) = (m|n) = [m \le n]\frac{m - 1}{\binom{n}{2}}

当*m*为定值时，这是一个*n*的似然函数
:\mathcal{L}(n) = [n \ge m]\frac{m - 1}{\binom{n}{2}}

总概率为

:\begin{align}
  \sum_{n}\mathcal{L}(n) &= \frac{m - 1}{1} \sum_{n=m}^\infty \frac{1}{\binom n 2} \\
                         &= \frac{m - 1}{1} \cdot \frac{2}{2 - 1} \cdot \frac{1}{\binom {m - 1}{2 - 1}} \\
                         &= 2
\end{align}

可信度质量分布函数为

:\begin{align}
       &(N=n\mid M=m,K=2) \\
  = {} &(n\mid m) \\
  = {} &\frac{\mathcal{L}(n)}{\sum_n \mathcal{L}(n)} \\
  = {} &[n \ge m]\frac{m - 1}{n(n - 1)}
\end{align}

中位数\scriptstyle \tilde{N}满足
:\sum_n [n \ge \tilde{N}](n|m) = \frac{1}{2}

所以
:\frac{m - 1}{\tilde N - 1} = \frac{1}{2}

因而中位数为
:\tilde{N} = 2m - 1

但*N*的均值为无穷大
:\mu = \sum_n n \cdot (n|m) = \frac{m - 1}1\sum_{n=m}^\infty \frac{1}{n - 1} = \infty

### 多辆坦克

#### 可信度质量分布函数
在序列号{1,...,*n*}中观测到的最大值为*k*的条件概率*m*为
:\begin{align}
       &(M=m|N=n,K=k\ge 2) \\
  = {} &(m\mid n,k) \\
  = {} &[m\le n]\frac{\binom{m - 1}{k - 1}}{\binom{n}{k}}
\end{align}

*n*的似然函数表达式与此相同
:\mathcal{L}(n) = [n \ge m]\frac{\binom{m - 1}{k - 1}}{\binom{n}{k}}

*k* ≥ 2的总概率为一个有限值：
:\begin{align}
  \sum_n \mathcal{L}(n)
    &= \frac{\binom{m - 1}{k - 1}}{1} \sum_{n=m}^\infty {1 \over \binom n k} \\
    &= \frac{\binom{m - 1}{k - 1}}{1} \cdot \frac{k}{k-1} \cdot \frac{1}{\binom{m - 1}{k - 1}} \\
    &= \frac k{k - 1}
\end{align}

可信度质量分布函数为
:\begin{align}
       &(N=n|M=m,K=k \ge 2) = (n|m,k) \\
  = {} &\frac{\mathcal{L}(n)}{\sum_n \mathcal{L}(n)} \\
  = {} &[n\ge m]\frac{k-1}{k} \frac{\binom{m - 1}{k - 1}}{\binom n k} \\
  = {} &[n\ge m]\frac{m-1}{n} \frac{\binom{m - 2}{k - 2}}{\binom{n - 1}{k - 1}} \\
  = {} &[n\ge m]\frac{m-1}{n} \frac{m - 2}{n - 1} \frac{k - 1}{k - 2} \frac{\binom{m - 3}{k - 3}}{\binom{n-2}{k-2}}
\end{align}

互补累积分布函数是*N* > *x*的可信度
:\begin{align}
       &(N>x\mid M=m,K=k) \\
  = {} &\begin{cases}
                                    1 &\text{if }x   

累积分布函数是*N* ≤ *x*的可信度
:\begin{align}
       &(N\le x|M=m,K=k) \\
  = {} &1 - (N>x\mid M=m,K=k) \\
  = {} &[x \ge m]\left(1 - \frac{\binom{m - 1}{k - 1}}{\binom{x}{k - 1}}\right)
\end{align}

#### 数量级
敌方坦克数目的数量级为
:\begin{align}
  \mu &= \sum_n n\cdot(N=n|M=m,K=k) \\&
       = \sum_n n [n\ge m]\frac {m-1}n \frac {\binom{m-2}{k-2}}{\binom{n-1}{k-1}}  \\&
       = \frac{m-1}1 \frac{\binom{m-2}{k-2}}1\sum_{n=m}^\infty \frac 1{\binom{n-1}{k-1}}\\&
       = \frac{m-1}1 \frac{\binom{m-2}{k-2}}1 \cdot \frac{k-1}{k-2}\frac {1}{\binom{m-2}{k-2}}\\&
       = \frac{m-1}1 \frac{k-1}{k-2}
\end{align}

#### 统计不确定度
统计的不确定度用标准差*σ*表示，其满足等式
:\sigma^2 + \mu^2 = \sum_n n^2 \cdot (N=n\mid M=m,K=k)

所以
:\begin{align}
  \sigma^2+\mu^2-\mu &= \sum_n n(n-1)\cdot(N=n\mid M=m,K=k)\\&
                      = \sum_{n=m}^\infty n(n-1)\frac{m-1}n \frac{m-2}{n-1} \frac{k-1}{k-2} \frac{\binom{m-3}{k-3}}{\binom{n-2}{k-2}}\\&
                      = \frac{m-1}1 \frac{m-2}1 \frac{k-1}{k-2} \cdot \frac{\binom{m-3}{k-3}}1 \sum_{n=m}^\infty \frac 1{\binom{n-2}{k-2}}\\
& = \frac{m-1}1 \frac{m-2}1 \frac{k-1}{k-2} \frac{\binom{m-3}{k-3}}1 \frac{k-2}{k-3} \frac 1{\binom{m-3}{k-3}}\\
& = \frac{m-1}1 \frac{m-2}1 \frac{k-1}{k-3}\\&
\end{align}

及
:\begin{align}
  \sigma &= \sqrt{\frac{m-1}1 \frac{m-2}1 \frac{k-1}{k-3}+\mu-\mu^2} \\&
          = \sqrt{\frac{(k-1)(m-1)(m-k+1)}{(k-3)(k-2)^2}} \\&
\end{align}

方差均值比则为
:\frac{\sigma^2}\mu = \frac{m - k + 1}{(k - 3)(k - 2)}

## 参见
* 标记重捕法，另一种估计样本总数的方法
* 最大间距估计，对“假设均匀分布”直觉的推广
* 哥白尼原则，类似的预测，基于单个样本（目前的年龄）计算预期寿命。
** 一个应用：末日论，估计人类的预期生存时间

### 其他对估计的讨论
* 最大似然估计
* 估计偏差#离散均匀分布的最大值
* 似然函数#例2

## 参考
;注释

;引用

;文献

* 
* 
* 
*