使用python進行貝葉斯統(tǒng)計分析|附代碼數(shù)據(jù)
原文鏈接：http://tecdat.cn/?p=7637
最近我們被客戶要求撰寫關(guān)于貝葉斯統(tǒng)計的研究報告，包括一些圖形和統(tǒng)計輸出。
本文講解了使用PyMC3進行基本的貝葉斯統(tǒng)計分析過程.?（?點擊文末“閱讀原文”獲取完整代碼數(shù)據(jù)********?）。
# Importsimport pymc3 as pm # python的概率編程包import numpy.random as npr # numpy是用來做科學計算的import numpy as npimport matplotlib.pyplot as plt # matplotlib是用來畫圖的import matplotlib as mplfrom collections import Counter # ? import seaborn as sns # ? # import missingno as msno # 用來應對缺失的數(shù)據(jù)# Set plotting style# plt.style.use('fivethirtyeight')sns.set_style('white')sns.set_context('poster')%load_ext autoreload%autoreload 2%matplotlib inline%config InlineBackend.figure_format = 'retina'import warningswarnings.filterwarnings('ignore')
**
拓端
，贊37
使用python進行貝葉斯統(tǒng)計分析?
??
貝葉斯公式?
貝葉斯主義者的思維方式?
根據(jù)證據(jù)不斷更新信念
pymc3?
常見的統(tǒng)計分析問題?
參數(shù)估計: "真實值是否等于X"
比較兩組實驗數(shù)據(jù): "實驗組是否與對照組不同? "
問題1: 參數(shù)估計?
"真實值是否等于X?"
或者說
"給定數(shù)據(jù)，對于感興趣的參數(shù)，可能值的概率分布是多少？"
例 1: 拋硬幣問題?
我把我的硬幣拋了?_n_次，正面是?_h_次。這枚硬幣是有偏的嗎？
參數(shù)估計問題parameterized problem?
先驗假設?
對參數(shù)預先的假設分布:??p～Uniform(0,1)
likelihood function(似然函數(shù), 翻譯這詞還不如英文原文呢):?data～Bernoulli(p)
# 產(chǎn)生所需要的數(shù)據(jù)from random import shuffletotal = 30n_heads = 11n_tails = total - n_headstosses = [1] * n_heads + [0] * n_tailsshuffle(tosses)
數(shù)據(jù)?
fig = plot_coins()plt.show()
點擊標題查閱往期內(nèi)容
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
左右滑動查看更多
01
02
03
04
MCMC Inference Button (TM)?
100%|██████████| 2500/2500 [00:00<00:00, 3382.23it/s]
結(jié)果?
pm.traceplot(coin_trace)plt.show()
In?[10]:
plt.show()
95% highest posterior density (HPD, 大概類似于置信區(qū)間)?包含了?region of practical equivalence (ROPE, 實際等同區(qū)間).
例 2: 藥品活性問題?
我有一個新開發(fā)的分子X; X在阻止流感病毒復制方面有多好？
實驗?
測試X的濃度范圍, 測量流感活性
計算?IC50: 能夠抑制病毒復制活性50%的X濃度.
data?
import pandas as pdchem_df = pd.DataFrame(chem_data)chem_df.columns = ['concentration', 'activity']chem_df['concentration_log'] = chem_df['concentration'].apply(lambda x:np.log10(x))# df.set_index('concentration', inplace=True)
參數(shù)化問題parameterized problem?
給定數(shù)據(jù), 求出化學物質(zhì)的IC50值是多少, 并且求出置信區(qū)間( 原文中the uncertainty surrounding it, 后面看類似置信區(qū)間的含義)?
先驗知識?
由藥學知識已知測量函數(shù)(measurement function):??m=β1+ex?IC50
測量函數(shù)中的參數(shù)估計, 來自先驗知識:?β～HalfNormal(1002)
關(guān)于感興趣參數(shù)的先驗知識:?log(IC50)～ImproperFlat
likelihood function:?data～N(m,1)
數(shù)據(jù)?
In?[13]:
fig = plot_chemical_data(log=True)plt.show()
MCMC Inference Button (TM)?
In?[16]:
pm.traceplot(ic50_trace[2000:], varnames=['IC50_log10', 'IC50']) ?# live: sample from step 2000 onwards.plt.show()
結(jié)果?
In?[17]:
pm.plot_posterior(ic50_trace[4000:], varnames=['IC50'], ? ? ? ? ? ? ? ? ?color='#87ceeb', point_estimate='mean')plt.show()
該化學物質(zhì)的 IC50?大約在[2 mM, 2.4 mM] (95% HPD). 這不是個好的藥物候選者. 在這個問提上不確定性影響不大, 看看單位數(shù)量級就知道IC50在毫摩的物質(zhì)沒什么用...
第二類問題: 實驗組之間的比較?
"實驗組和對照組之間是否有差別? "
例 1: 藥品對IQ的影響問題?
藥品治療是否影響(提高)IQ分數(shù)?
def ECDF(data): ? ?x = np.sort(data) ? ?y = np.cumsum(x) / np.sum(x) ? ? ? ?return x, ydef plot_drug(): ? ?fig = plt.figure() ? ?ax = fig.add_subplot(1,1,1) ? ?x_drug, y_drug = ECDF(drug) ? ?ax.plot(x_drug, y_drug, label='drug, n={0}'.format(len(drug))) ? ?x_placebo, y_placebo = ECDF(placebo) ? ?ax.plot(x_placebo, y_placebo, label='placebo, n={0}'.format(len(placebo))) ? ?ax.legend() ? ?ax.set_xlabel('IQ Score') ? ?ax.set_ylabel('Cumulative Frequency') ? ?ax.hlines(0.5, ax.get_xlim()[0], ax.get_xlim()[1], linestyle='--') ? ? ? ?return fig
In?[19]:
# Eric Ma自己很好奇, 從頻率主義的觀點, 差別是否已經(jīng)是具有"具有統(tǒng)計學意義"from scipy.stats import ttest_indttest_ind(drug, placebo) # (非配對) t檢驗. P=0.025, 已經(jīng)<0.05了
Out[19]:
Ttest_indResult(statistic=2.2806701634329549, pvalue=0.025011500508647616)
實驗?
參與者被隨機分為兩組:
給藥組?vs.?安慰劑組測量參與者的IQ分數(shù)
先驗知識?
被測數(shù)據(jù)符合t分布:??data～StudentsT(μ,σ,ν)
以下為t分布的幾個參數(shù):
均值符合正態(tài)分布:??μ～N(0,1002)
自由度(degrees of freedom)符合指數(shù)分布:??ν～Exp(30)
方差是positively-distributed:??σ～HalfCauchy(1002)
數(shù)據(jù)?
In?[20]:
fig = plot_drug()plt.show()
代碼?
In?[21]:
y_vals = np.concatenate([drug, placebo])labels = ['drug'] * len(drug) + ['placebo'] * len(placebo)data = pd.DataFrame([y_vals, labels]).Tdata.columns = ['IQ', 'treatment']
MCMC Inference Button (TM)?
結(jié)果?
In?[24]:
pm.traceplot(kruschke_trace[2000:], ? ? ? ? ? ? varnames=['mu_drug', 'mu_placebo'])plt.show()
In?[25]:
pm.plot_posterior(kruschke_trace[2000:], color='#87ceeb', ? ? ? ? ? ?varnames=['mu_drug', 'mu_placebo', 'diff_means'])plt.show()
IQ均值的差距為: [0.5, 4.6]
頻率主義的 p-value:?0.02?(!!!!!!!!)
注: IQ的差異在10以上才有點意義. p-value=0.02說明組間有差異, 但沒說差異有多大. 這個故事說的是雖然有差異, 但是差異太小了, 也沒啥意思.
In?[27]:
ax = adjust_forestplot_for_slides(ax)plt.show()
森林圖：在同一軸上的95％HPD（細線），IQR（粗線）和后驗分布的中位數(shù)（點），使我們能夠直接比較治療組和對照組。
In?[29]:
ax = pm.plot_posterior(kruschke_trace[2000:], ? ? ? ? ? ? ? ? ? ? ? varnames=['effect_size'], ? ? ? ? ? ? ? ? ? ? ? color='#87ceeb')overlay_effect_size(ax)
效果大?。–ohen's d, 效果微小, 效果中等, 效果很大）可以從微小到很大（95％HPD [0.0，0.77]）。
?這種藥很可能是無關(guān)緊要的。
沒有生物學意義的證據(jù)。
例 2: 手機消毒問題??
比較兩種常用的消毒方法, 和我的fancy方法, 哪種消毒方法更好
實驗設計?
將手機隨機分到6組: 4 "fancy" 方法 + 2 "control" 方法.
處理前后對手機表面進行拭子菌培養(yǎng)
count?菌落數(shù)量, 比較處理前后的菌落計數(shù)
Out[30]:
sample_id ? ? ? ? ? ? ? ? int32treatment ? ? ? ? ? ? ? ? int32colonies_pre ? ? ? ? ? ? ?int32colonies_post ? ? ? ? ? ? int32morphologies_pre ? ? ? ? ?int32morphologies_post ? ? ? ? int32year ? ? ? ? ? ? ? ? ? ?float32month ? ? ? ? ? ? ? ? ? float32day ? ? ? ? ? ? ? ? ? ? float32perc_reduction morph ? ?float32site ? ? ? ? ? ? ? ? ? ? ?int32phone ID ? ? ? ? ? ? ? ?float32no case ? ? ? ? ? ? ? ? float32frac_change_colonies ? ?float64dtype: object
數(shù)據(jù)?
In?[32]:
fig = plot_colonies_data()plt.show()
先驗知識?
菌落計數(shù)符合泊松Poisson分布. 因此...
菌落計數(shù)符合泊松分布:??dataij～Poisson(μij),j∈[pre,post],i∈[1,2,3...]
泊松分布的參數(shù)是離散均勻分布:??μij～DiscreteUniform(0,104),j∈[pre,post],i∈[1,2,3...]
滅菌效力通過百分比變化測量，定義如下:??mupre?mupostmupre
MCMC Inference Button (TM)?
In?[34]:
with poisson_estimation: ? ?poisson_trace = pm.sample(200000)
Assigned Metropolis to pre_musAssigned Metropolis to post_mus100%|██████████| 200500/200500 [01:15<00:00, 2671.98it/s]
In?[35]:
pm.traceplot(poisson_trace[50000:], varnames=['pre_mus', 'post_mus'])plt.show()
結(jié)果?
In?[39]:
pm.forestplot(poisson_trace[50000:], varnames=['perc_change'], ? ? ? ? ? ? ?ylabels=treatment_order) #, xrange=[0, 110])plt.xlabel('Percentage Reduction')ax = plt.gca()ax = adjust_forestplot_for_slides(ax)
點擊文末?“閱讀原文”
獲取全文完整資料。
本文選自《使用python進行貝葉斯統(tǒng)計分析》。
點擊標題查閱往期內(nèi)容
數(shù)據(jù)分享|R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python用MCMC馬爾科夫鏈蒙特卡洛、拒絕抽樣和Metropolis-Hastings采樣算法
R語言貝葉斯METROPOLIS-HASTINGS GIBBS 吉布斯采樣器估計變點指數(shù)分布分析泊松過程車站等待時間
R語言馬爾可夫MCMC中的METROPOLIS HASTINGS，MH算法抽樣（采樣）法可視化實例
python貝葉斯隨機過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實現(xiàn)
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機波動率模型、序列蒙特卡羅SMC、M H采樣分析時間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機波動率SV模型、粒子濾波、Metropolis Hasting采樣時間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進球數(shù)
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病
R語言中貝葉斯網(wǎng)絡（BN）、動態(tài)貝葉斯網(wǎng)絡、線性模型分析錯頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負擔能力數(shù)據(jù)集
R語言實現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學術(shù)能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計
