贝叶斯思维（实例2）——估计

火车头问题

铁路上以1到N命名火车头。有一天你看到一个标号60的火车头，请估计铁路上有多少火车头？

应用贝叶斯进行推理，可以将这个问题分成两步：
1. 在得到数据之前，我们对N的认识是什么？
2. 已知一个N的任意值后，得到数据（“标志为60号的火车头”）的似然度？

第一个问题的答案就是问题的前置概率，第二个问题是似然度。

在选择前置概率上，我们还没有太多的基本信息，先假设N可以是从1到1000等概率的任何值。

hypos = xrange(1, 1001)

接着需要的是似然函数。先假设存在一个有N个火车头的车队，我们看到60号火车头的概率是多少？假设只有一个列车运营公司，看到任一个火车头有同等可能，那么看到任何特定火车头的机会为1/N。似然度函数如下：

import logging#定义字典对象class _DictWrapper(object):    #具体代码见《贝叶斯思维（实例1）》    def InitSequence(self, values):        for value in values:            self.Set(value, 1)    def __len__(self):        return len(self.d)#概率质量函数对象class Pmf(_DictWrapper):     #具体代码见《贝叶斯思维（实例1）》def Probability(o):    return o / (o + 1)class Suite(Pmf):    def Update(self, data):        for hypo in self.Values():            like = self.Likelihood(data, hypo)            self.Mult(hypo, like)        return self.Normalize()class Train(Suite):    def Likelihood(self, data, hypo):        if hypo<data:            return 0        else:            return 1.0/hypo

Update如下：

suite = Train(hypos)suite.Update(60)

结果为0.003。

如果非要猜，最优可能的值是60。不过，这不是我们的目标。另一个可选的方法是计算后验概率的平均分布：

def Mean(suite):    total  = 0    for hypo, prob in suite.Items():        total += hypo*prob    return totalprint Mean(suite)

后验的平均值是333，所以如果想最大限度地减少错误，这也许是个好的猜测结果。

为了进一步解决火车头问题，我们必须进行一些假设。
但仅仅依赖一次观察的小量数据，情况可能如此敏感。回忆一下，从1到1000的均匀分布的先验概率，后验概率的平均值是333。同样上界为500，得到的后验平均值为207，一个2000的上界后验平均值为552。

所以结果很糟。有两种方法继续进行分析：
- 获取更多的数据
- 更多的背景信息

有了更多的数据后，基于不同的先验概率，后验分布趋于收敛。
例如，假设除了列车60我们也看到列车30和90,我们可以这样更新分布：

for data in[60, 30, 90]:    suite.Update (data)

采样这些数据时，后验概率的均值时：

上限 后验均值 500 152 1000 164 2000 171

这样差异就较小了。

如果没有更多的数据，另一个方法是通过收集背景资料优化先验。
即使没有深入了解铁路产业的一些具体情况，我们也可以猜测，在大多数领域，有大量小公司，一些中型公司，一个到两个非常大型的公司。事实上，公司规模的分布往往遵循幂律。这规律表明，如果少于10个火车头的公司有1000家，100个火车头的公司可能有100家，1000个火车头的公司有10家，10000个火车头的公司可能仅有1家。

在数学上，幂律表示公司的数量与公司规模成反比，或

PMF(x)∝(1x)α

其中PMF(x)是x的概率质量函数，α是一个通常接近于1的参数。
我们可以构造一个服从幂律分布的先验如下：

class Dice(Suite):    def Likelihood(self, data, hypo):        if hypo < data:            return 0        else:            return 1.0/hypoclass Train (Dice):    def __init__ (self, hypos, alpha = 1.0):        Pmf.__init__ (self)        for hypo in hypos:            self.Set (hypo, hypos** (-alpha))        self.Normalize()

而下面是生成前置概率的代码：

hypos = range(1, 1001)suite = Train(hypos)

如果基于这种先验概率，在观察到列车30、60和90时，后验概率分别是：

上限 后验均值 500 131 1000 133 2000 134

事实上，考虑一个任意大的上界，平均值都收敛于134。所以基于幂律分布的先验概率是比较现实的，因为它基于公司规模的一般情况，并且在实际中表现的更好。

置信区间
对于点估计，通常使用平均数、中位数或最大似然值。
对于区间，通常给出两个计算值，使得未知量有90%的可能落入这两个值之间。这些值定义了一个置信区间。

计算置信区间的一个简单方法是在后验概率分布中累积其中的概率，并记录对应于概率5%和95%的值。thinkbayes提供了一个函数计算百分位数：

def Percentile(pmf, percentage):    p = percentage / 100.0    total = 0    for val, prob in pmf.Items():        total += prob        if total >= p:            return val

下面是一个其应用代码：

interval = Percentile(suite, 5), Percentile(suite, 95)print interval

在前面的示例（看到了三个火车，且呈幂律分布的先验概率的火车头问题）中90%置信区间为(91, 243)。如此大的范围其实确切的表明，我们仍然相当不确定究竟有多少火车头存在。

贝叶斯当中，有两种途径选择先验分布。一些人建议选择最能代表问题相关背景资料的先验概率，这种情况下，先验被认为是“信息”。另一种方法是所谓的“无信息参考的先验”，其目的是为了让数据来说话，越没约束越好。

参考文献：
贝叶斯思维. Allen B.Downey. 许杨毅 译