贝叶斯思维——chapter3(估计)

贝叶斯思维——chapter3(估计)

---

写在前面：
库文件链接： thinkbayes.py

---

3.1 骰子问题

假设我有５颗骰子，分别为4面，6面，8面，12面，20面。
现在我从里面随机选取一个，投掷以后得到6，问我最有可能取出的是哪个骰子？
我们通过一下3个步骤来解决这个问题：

1. 选择假设的表示方法
2. 选择数据的表示方法
3. 编写似然度函数

现在我们用[4,6,8,12,20]来表示假设：

suite = Dice([4,6,8,12,20])

似然函数：

def Likelihood(self,data,hypo):#data为得到的投掷结果，hypo是我们的假设
    if hypo < data:
    ＃意味着骰子的投掷值要大于骰子的面数，显然这不可能
        return 0
    else
    ＃可能情况下得到投掷值的概率自然就是1/hypo
        return 1.0/hypo

随着投掷数据的变动，我们可以进一步得到每一个假设的后验概率分布，
先验实际上也能看成没有任何事件发生情况下的假设的可能性。

$$贝叶斯公式： P(hopy|data) = \frac{P(hopy)P(data|hopy)}{P(data)}\\ 后验概率和为1:\sum_{i=1}^nP(hopy_i|data)=1$$

class Dice(suite):

    def __init__(self,hypos):#对每个假设赋值先验概率,等可能性
        for hypo in hypos:
            self.Set(hypo,1)
        ＃数据归一化
        self.Normalize()

    def Update(self,data):
        for hypo in self.Values():
            #我们得到了hypo假设下，数据data发生的概率，也就是我们说的该假设下的似然度（至少数值上是相等的）
            like = self.Likelihood(data,hypo)
            #事实上,我们在不断得更新假设的概率，这里我们没有取计算P(data)的概率。事实上，在确定假设之前P(data)是未知的，我么直接计算的P(hopy)P(data|hopy),因为我们在接下来进行了归一化，避免了P(data)的计算
            self.Mult(hypo,like)
        self.Normalize()

所以此时，我们投掷得到6，就相当于：

suite.Update(6)

得到：

$$\begin{array}{c|rr} hypo& \text{先验概率} & \text{data=6}\\ \hline 4 & 0.2 & 0.0\\ 6 & 0.2 & 0.39\\ 8 & 0.2 & 0.29 \\ 12 & 0.2 & 0.196 \\ 20 & 0.2 & 0.117 \end{array}$$
现在我们可以得到问题的解为6。

3.2 火车头的问题

铁路上以1到N命名火车，有一天你看到一个标号为60的火车头，估计铁路上有多少火车头？
解：可以肯定的是铁路上至少有60个火车头，但这个数字到底是多少？
我们可以把问题分成两个步骤：

1. 在得到数据之前，我们对N的人事认识是什么？（先验概率）
2. 已知一个N的任意值后，得到数据(“标志为60号的火车头”)的似然度？（似然度）
   我们假设，N是从1到1000的等概率的任何值。

hypos = range(1,1000)#假设

其实到这里，应该已经看出来了，我们进入到了跟投掷骰子一样的情况里面。
这里我们可以得到的是：
可能性最大化的解是60
不过，这还不是我们的目标，另一个可选的方案是计算后验概率的平均值分布：

听得莫名其妙，简单点说，就是期望。

$$E(火车数量)＝\sum_{i=1}^nhopy_{火车数量_i}P(火车数量_i)$$
可以得到在取假设1——1000的情况下，我们能得到期望为333

3.3 怎样看待先验概率

如果细看上面的火车头问题的话，会发现我们的假设是非常武断的

1. 我们的假设区间的先验假设是平均分布的
2. 选择1000作为上界
   但我们选择1000作为上界时，我们的期望值是333;
   但当我们选择500作为上界时，我们的期望值是207;
   2000为上界是，期望值是552。
   结论：猜测结果对上界敏感

有两种方法可以进一步减小我们误差：

1. 获取更多的数据。（更多的数据，基于不同的先验分布，后验分布会趋于收敛）
2. 更多的背景资料（先验概率）

假设，我们还看到了火车头30和90,我们就能向骰子一样更新分布。
采用数据后，期望值是：

$$\begin{array}{c|r} 上限值& \text{期望值}\\ \hline 500 & 152 \\ 1000 & 164 \\ 2000 & 171 \\ \end{array}$$
期望值的差异变小了。

3.4 其他先验概率

如果没有更多的数据，我们就需要是通过收集背景资料来优化先验。
事实上，公司规模的分布往往遵循幂律，这项规律表明，如果少于10个火车头的公司有1000家，100个火车头的公司可能有100家，1000个火车头的公司有10家,10000个火车头的公司可能只有1家。
数学上，幂律表示公司数量与公司规模成反比，或：

$$PMF(x)\propto\left(\frac{1}{x}\right)^\alpha$$
其中 $PMF(x)$是x的概率密度函数， $\alpha$通常是一个接近于1的参数。
我们可以构造一个服从幂律分布的先验：

class train():
    def __init__(self, hypos, alpha = 1.0):
        for hypo in hypos:
            self.Set(hypo, hypo**(-alpha))
        self.Normalize()

再次说明，上限是任意的，但对于幂律分布的先验概率，后验概率对这一选择的敏感度较小。
如果基于这种先验概率，在观察到列车30,60,90时，期望分别是：

$$\begin{array}{c|r} 上限值& \text{期望值}\\ \hline 500 & 131 \\ 1000 & 133 \\ 2000 & 134 \\ \end{array}$$
事实上，考虑一个任意大的上界，期望都将收敛到134。
所以基于幂律分布的先验概率是比较现实的，因为它是基于公司规模的一些情况，并且在实际中表现更好。

3.5 置信区间

一旦计算出了后验分布，通过单点估计或区间对后验分布进行总结通常是有用的。
点估计：
1. 平均数
2. 中位数
3. 最大似然值

区间，我们通常需要计算两个值，使得未知量有90%的可能落在两个值之间。这些值定义了一个置信区间。
计算置信区间的简单做法是在后验概率分布中累加其中的概率，并记录对应与概率5%和95%的值。

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)。如此大的区间确切表明，（尽管期望值已经收敛了）我们任然相当不确定究竟有多少个火车头的存在。