数学概率统计29 pymc入门

贝叶斯思维¶

和更传统的统计推断不同，贝叶斯推断会保留不确定性，在贝叶斯派的世界观中，概率是被解释为我们对一件事情发生的相信程度或者说信心 (飞机事故，总统选举)。

需注意的是，我们每个人都可以给事件赋概率值，而不是存在某个唯一的概率值，因为不同的人拥有不同的信息，因此他们对同一事件发生的信心也可以有不同的值，但这些不同并不说明其他人是错误的。

飞机事故：综合某航空公司过去十年出现飞机事故的次数，来推断该公司出现飞机事故的概率;
总统选举：根据候选人的声望、民意等等，推断他能当上总统的概率。

概率编程入门实例——贝叶斯点估计¶

我们来以《贝叶斯统计》 (韦来生著) 这本书上的一个习题作为我们入门的实例，因为学习 PyMC3 必然是要落地的，通过书上的问题能帮我们迅速拉近和PyMC3 的距离。

观察问题¶

设 $X_1, \cdots, X_5$ 为从同一批生产的电子零件中抽取的样本, 进行检验以确定其平均寿命. 假定 $X_i \sim \operatorname{Exp}(1 / \theta)$, 由过去的先验信息可知 $\theta$ 的先验分布为逆伽玛分布 $\Gamma^{-1}(10,100)$. 令 5 个具体观测值为 $5,12,14,10,12$, 求 $\theta$ 的广义最大似然估计和后验均值估计, 以及它们各自的后验方差.

这里有五个样本 $X_1$ 到 $x_5$ ，而且已经有了观测值然后这些样本服从指数分布 $\operatorname{Exp}(1 / \theta)$ ，而 $\theta$ 的先验分布是逆伽马分布 $\Gamma^{-1}(10,100)$ ，然后我们要求什么呢，我们要求参数 $\theta$ 的广义最大似然估计，我们现在使用编程来解决这个问题。

确定观测数据¶

导入第三方库，以及确定观测数据:

import pymc3 as pm
X = [5, 12, 14, 10, 12]

定义模型与变量¶

接着第三步，定义模型，我们使用 pm.Model(）方法，以及 with 上下文管理器将该代码块中的变量都添加到 my_first_model 模型，注意，这个 with pm.Model() as my_first_model 不同于其他的有关 with 的使用方法，my_first_model 这个模型不会在退出with 时自动关闭，而是会一直保存，作为连接整个代码程序的枢纽:

with pm.Model() as my_first_model:

    theta = pm.InverseGamma('theta', alpha=10, beta=100)
    X_obs = pm.Exponential('X_obs', lam=1 / theta, observed=X)

确定了 $\theta$ 之后，我们要来确定需要拟合的变量，即 X_obs，它服从指数分布 Exp(1 / \theta)，因此传入的参数是 1 / theta，而参数 observed=X 是告诉模型这个变量的值是已经被观测到了的，不会被拟合算法改变，而观测值就是 X，在这X是一开始我们定义的列表。

定义了模型之后，我们就要开始对模型进行统计推断了，这里选择 MAP 作为推断方法，那什么是 MAP 呢，MAP（maximum a posteriori probability estimate）中文名称叫极大似然估计，其实它就是《贝叶斯统计》这本书上所讲的广义最大似然估计，即后验众数估计，PyMC3 也提供了对应的函数 find_MAP()，不过它的计算方式和我们书上说的不同，因为书上是通过分析计算的方式推断出参数的估计值，但是这种方式一般只适合于标准分布，而实际上我们碰到的很多分布都是非标准分布，很难或无法通过分析计算进行推断，这个时候就需要通过 find_MAP() 所采用的数值优化的方法得出参数的估计值。

map_estimate = pm.find_MAP(model=my_first_model, start={'theta':10})

In [2]:

import pymc3 as pm
X = [5, 12, 14, 10, 12]

WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

In [3]:

with pm.Model() as my_first_model:
    theta = pm.InverseGamma('theta', alpha=10, beta=100)
    X_obs = pm.Exponential('X_obs', lam=1 / theta, observed=X)

In [4]:

map_estimate = pm.find_MAP(model=my_first_model, start={'theta':10})

100.00% [5/5 00:00<00:00 logp = -18.891, ||grad|| = 0.7]

In [5]:

map_estimate

Out[5]:

{'theta_log__': array(2.25784936), 'theta': array(9.56250158)}