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已經收斂。