Think Bayes Note 5: 火车头问题
Contents [show]
1. 问题描述
铁路上有 N 个火车头,依次编号为 1 到 N,有一天你看到了一个编号为 60 的火车头,估计铁路上一共有多少个火车头。
2. 推导求解
在看到任何火车头之前,我们只知道铁路上有 N 个火车头,为了能够进一步分析,先假设火车头总数 N 不会超过 1000,且在 1 到 1000 的范围内均匀分布。假设 Hn 为铁路上有 n 个火车头,则有:
P(D|Hn)=11000,n=1,2,⋯,1000
记 D 为看到编号为 60 的火车头,D 发生之后,可以知道火车头的数量不可能少于 60,即有:
P(D|Hn)=0,n<60
P(Hn|D)=P(Hn)P(D|Hn)P(D)=0,n<60
如果一共有 60 个火车头,即 n=60 时,看到任意编号的火车头的概率都为 160,则看到编号为 60 的火车头的概率也为 160,即:
P(D|H60)=160
计算 P(D) 为:
P(D)=1000∑n=1P(Hn)P(D|Hn)=110001000∑n=601n≈0.002822
得到一共有 60 个火车头的后验概率为:
P(H60|D)=P(H60)P(D|H60)P(D)≈0.005905
使用相同方法可以得到任意 P(Hn|D),60≤n≤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 个火车头,即:
PMF(x)∝1xα
其中 PMF(x) 是 x 的概率质量函数,α 是一个接近于 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)