pymc_實現貝葉斯統計模型和馬爾科夫鏈蒙塔卡洛

 python機器學習-乳腺癌細胞挖掘（博主親自錄制視頻）

https://study.163.com/course/introduction.htm?courseId=1005269003&utm_campaign=commission&utm_source=cp-400000000398149&utm_medium=share

項目合作：QQ231469242

 

http://hao.jobbole.com/pymc/

 

PyMC是一個實現貝葉斯統計模型和馬爾科夫鏈蒙塔卡洛采樣工具擬合算法的Python庫。PyMC的靈活性及可擴展性使得它能夠適用於解決各種問題。除了包含核心采樣功能，PyMC還包含了統計輸出、繪圖、擬合優度檢驗和收斂性診斷等方法。

 

特性

PyMC使得貝葉斯分析盡可能更加容易。以下是一些PyMC庫的特性：

• 用馬爾科夫鏈蒙特卡洛算法和其他算法來擬合貝葉斯統計分析模型。
• 包含了大范圍的常用統計分布。
• 盡可能地使用了NumPy的一些功能。
• 包括一個高斯建模過程的模塊。
• 采樣循環可以被暫停和手動調整，或者保存和重新啟動。
• 創建包括表格和圖表的摘要說明。
• 算法跟蹤記錄可以保存為純文本，pickles，SQLite或MySQL數據庫文檔或HDF5文檔。
• 提供了一些收斂性診斷方法。
• 可擴展性：引入自定義的步驟方法和非常規的概率分布。
• MCMC循環可以嵌入在較大的程序中，結果可以使用Python進行分析。

 

 

使用

首先，在文件中定義你的模型，並命名為mymodel.py

# Import relevant modules import pymc import numpy as np # Some data n = 5*np.ones(4,dtype=int) x = np.array([-.86,-.3,-.05,.73]) # Priors on unknown parameters alpha = pymc.Normal('alpha',mu=0,tau=.01) beta = pymc.Normal('beta',mu=0,tau=.01) # Arbitrary deterministic function of parameters @pymc.deterministic def theta(a=alpha, b=beta): """theta = logit^{-1}(a+b)""" return pymc.invlogit(a+b*x) # Binomial likelihood for data d = pymc.Binomial('d', n=n, p=theta, value=np.array([0.,1.,3.,5.]),\ observed=Tr

# -*- coding: utf-8 -*-
"""
Created on Mon Jul 24 10:56:07 2017

@author: toby
"""

# Import relevant modules
import pymc
import numpy as np

# Some data
n = 5*np.ones(4,dtype=int)
x = np.array([-.86,-.3,-.05,.73])

# Priors on unknown parameters
alpha = pymc.Normal('alpha',mu=0,tau=.01)
beta = pymc.Normal('beta',mu=0,tau=.01)

# Arbitrary deterministic function of parameters
@pymc.deterministic
def theta(a=alpha, b=beta):
    """theta = logit^{-1}(a+b)"""
    return pymc.invlogit(a+b*x)

# Binomial likelihood for data
d = pymc.Binomial('d', n=n, p=theta, value=np.array([0.,1.,3.,5.]),\
                  observed=True)

 

測試腳本

# -*- coding: utf-8 -*-
"""
Created on Mon Jul 24 11:21:23 2017

@author: toby
"""

import pymc
import mymodel
 
S = pymc.MCMC(mymodel, db='pickle')
S.sample(iter=10000, burn=5000, thin=2)
pymc.Matplot.plot(S)

 

import pymcimport mymodelS = pymc.MCMC(mymodel, db='pickle')S.sample(iter=10000, burn=5000, thin=2)pymc.Matplot.plot

這個例子會產生10000個后驗樣本。這個樣本會存儲在Python序列化數據庫中。

 

教程示例

教程會指導用戶完成常見的PyMC應用。

如何用MCMC來擬合模型

PyMC提供了一些可以擬合概率模型的方法。最主要的擬合模型方法是MCMC，即馬爾科夫蒙特卡洛算法。生成一個MCMC對象來處理我們的模型，導入disaster_model.py並將其作為MCMC的參數。

 

調用MCMC中的sample()方法（或者交互采樣函數isample())來運行采樣器

# -*- coding: utf-8 -*-
"""
Created on Mon Jul 24 11:26:27 2017

@author: toby
"""

from pymc.examples import disaster_model
from pymc import MCMC
M = MCMC(disaster_model)
M.sample(iter=10000, burn=1000, thin=10)

等待幾秒鍾后，便可以看到采樣過程執行完成，模型已經完成擬合。

 

 

http://blog.csdn.net/dmsgames/article/details/52525636

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)

來自CODE的代碼片

snippet_file_0.txt

3、父變量與子變量

我們已經使用PyMC創建了統計模型，PyMC中提供方法查看模型中參數之間的關系，試例代碼如下：

from pymc.examples import disaster_model disaster_model.switchpoint.parents #顯示s的父參數 #輸出{'lower': 0, 'upper': 110} disaster_model.disasters.parents #顯示disasters的父參數 #輸出{'mu': <pymc.PyMCObjects.Deterministic 'rate' at 0x000000000B791BE0>} disaster_model.rate.children #顯示rate的子參數 #輸出{<pymc.distributions.new_dist_class.<locals>.new_class 'disasters' at 0x000000000B791C18>}

來自CODE的代碼片

snippet_file_0.txt

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)

來自CODE的代碼片

snippet_file_0.txt

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)

來自CODE的代碼片

snippet_file_0.txt

MCMC算法輸出模型中每個變量的樣本，獲得樣本方法如下：

M.trace('switchpoint')[:] #輸出array([43,43,44,...44,44])

來自CODE的代碼片

snippet_file_0.txt

畫出每個變量的采樣序列圖、后驗邊緣分布直方圖、自相關性圖，代碼如下：

from pymc.Matplot import plot plot(M)

來自CODE的代碼片

snippet_file_0.txt

采樣序列圖可以判斷MCMC是否收斂，如果采樣序列分布近似於白噪聲，那么可以判斷MCMC已經收斂。

 

 https://study.163.com/provider/400000000398149/index.htm?share=2&shareId=400000000398149（ 歡迎關注博主主頁，學習python視頻資源，還有大量免費python經典文章）