1、一个统计模型
有这样一个数据集,它按照时间顺序,收录了英国从1851年到1962年每年的矿难发生次数。如下图所示:
我们可以假设,矿难发生的概率服从一个Poisson过程,在某一年泊松过程的参数发生了变化,在时间轴的早些时候,矿难发生的概率较高,后来矿难发生的概率比较低。
我们将上述概念模型转化为统计模型:
以上模型参数定义如下:
- D_t: 第t年矿难发生的次数;
- r_t: 第t年Posson过程的参数;
- s: 泊松过程参数发生改变的那一年;
- e: 第s年之前,泊松过程的参数;
- l:第s年之后,泊松过程的参数;
- t_l,t_h: 年份t的下限和上限;
- r_e,r_l:e和l的先验分布
由于在模型中我们定义了D依赖于s,e,l,所以我们把D称作s,e,l的子变量,类似的,s,e,l称为D的父变量。
2、变量的两种类型
PyMC包中定义类两种随机变量类型,分别为stochastic和Deterministic。
模型中唯一的Deterministic变量是r,因为当我们知道r的父参数(s,l,e)后,我们可以准确地计算出r的值。
另一方面,s,D(在观察到数据之前)是stochastic变量,因为即使观察到他们的父变量,任然不能确定它们的值。
我们将模型写在一个名为 disaster_model.py 的python脚本中:
"""
导入numpy和pymc
"""
from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform
import numpy as np
"""
导入英国矿难数据集
"""
disasters_array = \
np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
"""
定义转折点s:
取值范围0-110
均匀离散分布
"""
switchpoint = DiscreteUniform('switchpoint', lower=0, upper=110, doc='Switchpoint[year]')
"""
定义e、l
指数分布
"""
early_mean = Exponential('early_mean', beta=1.)
late_mean = Exponential('late_mean', beta=1.)
"""
定义r
"""
@deterministic(plot=False)
def rate(s=switchpoint, e=early_mean, l=late_mean):
''' Concatenate Poisson means '''
out = np.empty(len(disasters_array))
out[:s] = e
out[s:] = l
return out
"""
定义矿难发生次数
服从泊松分布
"""
disasters = Poisson('disasters', mu=rate, value=disasters_array, observed=True)
3、父变量与子变量
我们已经使用PyMC创建了统计模型,PyMC中提供方法查看模型中参数之间的关系,试例代码如下:
from pymc.examples import disaster_model
disaster_model.switchpoint.parents #显示s的父参数
#输出{'lower': 0, 'upper': 110}
disaster_model.disasters.parents #显示disasters的父参数
#输出{'mu':
}
disaster_model.rate.children #显示rate的子参数
#输出{
.new_class 'disasters' at 0x000000000B791C18>}
4、变量的值
所有的PyMC变量都具有value属性,查看value值示例代码如下:
disaster_model.disasters.value
"""输出
array([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6, 3, 3, 5, 4, 5, 3, 1,
4, 4, 1, 5, 5, 3, 4, 2, 5, 2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3,
0, 0, 1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1, 0, 1, 0, 1, 0,
0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2, 3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2,
0, 0, 1, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
"""
disaster_model.switchpoint.value
#输出 array(40)
disaster_model.early_mean.value
#输出 array(1.1444157379406001)
disaster_model.late_mean.value
#输出 array(0.027985496189503425)
5、使用马尔科夫链蒙特卡洛(MCMC)拟合模型
PyMC提供MCMC方法拟合模型,使用方法如下:
from pymc.examples import disaster_model
from pymc import MCMC
M = MCMC(disaster_model)
M.sample(iter=10000, burn=1000, thin=10)
MCMC算法输出模型中每个变量的样本,获得样本方法如下:
M.trace('switchpoint')[:]
#输出array([43,43,44,...44,44])
from pymc.Matplot import plot
plot(M)
采样序列图可以判断MCMC是否收敛,如果采样序列分布近似于白噪声,那么可以判断MCMC已经收敛。