Loading [MathJax]/jax/output/HTML-CSS/jax.js

Think Bayes Note 5: 火车头问题

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)=1000n=1P(Hn)P(D|Hn)=110001000n=601n0.002822

  得到一共有 60 个火车头的后验概率为:

P(H60|D)=P(H60)P(D|H60)P(D)0.005905

  使用相同方法可以得到任意 P(Hn|D),60n1000

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
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
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_distvalues 的分布,默认为均匀分布。使用方法为:

hypos = range(1, 1001)
train = Train(hypos)
train.update(60)
print("H_60 = {}".format(train.prob(60)) )
print("mean = {}".format(train.mean()))
train.plot()
hypos = range(1, 1001) train = Train(hypos) train.update(60) print("H_60 = {}".format(train.prob(60)) ) print("mean = {}".format(train.mean())) train.plot()
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
H_60 = 0.005905417875729856 mean = 333.41989326370776
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
prob value 60 0.005905 61 0.005809 62 0.005715 63 0.005624 64 0.005536
           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()))
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()))
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
upper_bound = 500, mean = 207.08 upper_bound = 1000, mean = 333.42 upper_bound = 2000, mean = 552.18
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()))
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()))
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
upper_bound = 500, mean = 151.85 upper_bound = 1000, mean = 164.31 upper_bound = 2000, mean = 171.34
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()
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()
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()))
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()))
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
upper_bound = 500, mean = 130.71 upper_bound = 1000, mean = 133.28 upper_bound = 2000, mean = 134.00
upper_bound = 500, mean = 130.71
upper_bound = 1000, mean = 133.28
upper_bound = 2000, mean = 134.00

可见不同上界带来的差异更小了。

  使用 train.make_cdf() 可以生成 train 对应的累积分布函数类 Cdf。通过 Cdfpercentile() 方法可以得到 90% 的置信区间:

cdf = train.make_cdf()
print("({}, {})".format(cdf.percentile(5), cdf.percentile(95)))
cdf = train.make_cdf() print("({}, {})".format(cdf.percentile(5), cdf.percentile(95)))
cdf = train.make_cdf()
print("({}, {})".format(cdf.percentile(5), cdf.percentile(95)))

输出为:

(91, 243)
(91, 243)
(91, 243)