Think Bayes Note 5: 火车头问题
1. 问题描述
铁路上有 N 个火车头，依次编号为 1 到 N，有一天你看到了一个编号为 60 的火车头，估计铁路上一共有多少个火车头。
2. 推导求解
在看到任何火车头之前，我们只知道铁路上有 N 个火车头，为了能够进一步分析，先假设火车头总数 N 不会超过 1000，且在 1 到 1000 的范围内均匀分布。假设 $H_n$ 为铁路上有 $n$ 个火车头，则有：
\begin{equation}
P(D|H_n) = \frac{1}{1000}, \; n = 1, 2, \cdots, 1000
\end{equation}
记 D 为看到编号为 60 的火车头，D 发生之后，可以知道火车头的数量不可能少于 60，即有：
\begin{equation}
P(D|H_n) = 0, \; n < 60
\end{equation}
\begin{equation}
P(H_n|D) = \frac{P(H_n)P(D|H_n)}{P(D)} = 0, \; n < 60
\end{equation}
如果一共有 60 个火车头，即 $n = 60$ 时，看到任意编号的火车头的概率都为 $\frac{1}{60}$，则看到编号为 60 的火车头的概率也为 $\frac{1}{60}$，即：
\begin{equation}
P(D|H_{60}) = \frac{1}{60}
\end{equation}
计算 $P(D)$ 为：
\begin{equation}
P(D) = \sum_{n = 1}^{1000}P(H_n)P(D|H_n) = \frac{1}{1000} \sum_{n=60}^{1000}\frac{1}{n} \approx 0.002822
\end{equation}
得到一共有 60 个火车头的后验概率为：
\begin{equation}
P(H_{60}|D) = \frac{P(H_{60})P(D|H_{60})}{P(D)} \approx 0.005905
\end{equation}
使用相同方法可以得到任意 $P(H_n|D), \; 60 \leq n \leq 1000$。
3. 代码求解
定义 Train类如下：
class Train(Suite):
def __init__(self, values=None, prob_dist=lambda hypos: [1.0 / len(hypos) for _ in hypos]):
super().__init__(values, value_desc='Train number')
probs = prob_dist(values)
for hypo, prob in zip(self.values(), probs):
self.set(hypo, prob)
self.normalize()
def likelihood(self, data, hypo):
if hypo < data:
return 0
else:
return 1.0 / hypo
其中 __init__() 的参数 values 为假设火车头总数的列表，prob_dist 为 values 的分布，默认为均匀分布。使用方法为：
hypos = range(1, 1001)
train = Train(hypos)
train.update(60)
print("H_60 = {}".format(train.prob(60)) )
print("mean = {}".format(train.mean()))
train.plot()
输出为：
H_60 = 0.005905417875729856 mean = 333.41989326370776
图像为：
可见当看见编号为 60 的火车头后，火车数量少于 60 的概率为 0。使用 train.top(5) 得到最高概率的 5 个结果如下：
prob value 60 0.005905 61 0.005809 62 0.005715 63 0.005624 64 0.005536
如果想要以最大概率得到完全正确的答案，应当猜测一共有 60 个火车头。
前面通过 train.mean() 得到均值为 333。如果想要使猜测的均方误差最小，则应当以均值作为猜测。
4. 新的问题
前面假设火车的总数不超过 1000 且具有均匀分布，这些假设比较随意。分别假设火车总数不超过 500 、1000 和 2000，使用如下代码计算均值：
for upper_bound in [500, 1000, 2000]:
train = Train(range(1, upper_bound + 1))
train.update(60)
print("upper_bound = {}, mean = {:.2f}".format(upper_bound, train.mean()))
输出为：
upper_bound = 500, mean = 207.08 upper_bound = 1000, mean = 333.42 upper_bound = 2000, mean = 552.18
几个均值的差距很大，对假设的火车总数敏感。如果不同的人对火车总数做出的不同的假设，就会得到差异很大的结果。
为了得到更可靠的预测，可以：
- 收集更多数据
- 收集更多背景知识
5. 收集更多数据
假设我们又看到了编号为 30 和 90 的火车头，计算此时的后验概率和均值：
for upper_bound in [500, 1000, 2000]:
train = Train(range(1, upper_bound + 1))
for data in [60, 30, 90]:
train.update(data)
print("upper_bound = {}, mean = {:.2f}".format(upper_bound, train.mean()))
输出为：
upper_bound = 500, mean = 151.85 upper_bound = 1000, mean = 164.31 upper_bound = 2000, mean = 171.34
此时几个均值之间的差异就小多了。
6. 收集更多背景知识
假设我们通过调研相关背景资料，发现一个行业内的公司的规模往往遵循幂律，即公司数量与公司规模成反比。例如在铁路运输行业，有 1000 家公司拥有 10 个火车头，有 100 家公司拥有 100 个火车头，有 10 家公司拥有 1000 个火车头，有 1 家公司拥有 10000 个火车头，即：
\begin{equation}
PMF(x) \propto \frac{1}{x^\alpha}
\end{equation}
其中 $PMF(x)$ 是 $x$ 的概率质量函数，$\alpha$ 是一个接近于 $1$ 的参数。
使用一下代码构造一个生成幂律分布的函数，并以此作为 Train 类的分布。
def zipf_dist(hypos, alpha=1.0):
return [hypo ** (-alpha) for hypo in hypos]
train = Train(range(1, 1001), prob_dist=zipf_dist)
train.update(60)
train.plot()
图像为：
分别假设火车总数不超过 500 、1000 和 2000，再次计算均值：
for upper_bound in [500, 1000, 2000]:
train = Train(range(1, upper_bound + 1), prob_dist=exp_dist)
for data in [30, 60, 90]:
train.update(data)
print("upper_bound = {}, mean = {:.2f}".format(upper_bound, train.mean()))
输出为：
upper_bound = 500, mean = 130.71 upper_bound = 1000, mean = 133.28 upper_bound = 2000, mean = 134.00
可见不同上界带来的差异更小了。
使用 train.make_cdf() 可以生成 train 对应的累积分布函数类 Cdf。通过 Cdf 的 percentile() 方法可以得到 90% 的置信区间：
cdf = train.make_cdf()
print("({}, {})".format(cdf.percentile(5), cdf.percentile(95)))
输出为：
(91, 243)