參考：https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers

1、模型假設

貝葉斯線性回歸模型，假設參數具有正態分布（Normal）先驗，要預測的具有正態分布的觀測值?Y?，其期望?μ?是兩個預測變量的線性組合，?X1?和?X2?：

Y～N(μ,σ2)

μ=α+β1X1+β2X2Y

其中?α?是截距，?βi是變量?Xi的系數；?σ?代表觀察誤差。

模型中的未知變量需要賦予先驗分布，假設服從零均值的正態先驗。其中，系數?α?、?βi的方差為100，選擇半正態分布作為觀測誤差?σ的先驗分布。

α～N(0,100)

βi～N(0,100)

σ～|N(0,1)|

2、模擬數據生成

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3Dnp.random.seed(123)alpha=1 sigma=1 beta =[1, 2.5]N=100X1=np.random.randn(N) X2=np.random.randn(N)Y=alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(N)*sigma%matplotlib inline fig1,ax1 = plt.subplots(1, 2, figsize=(10,4)); ax1[0].scatter(X1, Y);ax1[0].set_xlabel('X1');ax1[0].set_ylabel('Y'); ax1[1].scatter(X2, Y);ax1[1].set_xlabel('X2');ax1[1].set_ylabel('Y');fig2 = plt.figure(2); ax2 = Axes3D(fig2); ax2.scatter(X1,X2,Y); ax2.set_xlabel('X1'); ax2.set_ylabel('X2'); ax2.set_zlabel('Y');

2、利用pymc3定義模型

import pymc3 as pm with pm.Model() as basic_model:alpha=pm.Normal('alpha',mu=0,sd=10)beta=pm.Normal('beta',mu=0,sd=10,shape=2)sigma=pm.HalfNormal('sigma',sd=1)mu=alpha+beta[0]*X1+beta[1]*X2Y_obs=pm.Normal('Y_obs',mu=mu,sd=sigma,observed=Y)

定義了三個具有正態分布先驗的隨機性隨機變量（stochastic?random veriables）。

定義了一個確定性隨機變量（deterministic?random veriable）,確定性指這個變量的值完全由右端值確定。

最后定義了?Y?的觀測值，這是一個特殊的觀測隨機變量（observed stochastic），表示模型數據的可能性。通過?observed?參數來告訴這個變量其值是已經被觀測到了的（就是?Y?變量），不會被擬合算法改變。可以用?numpy.ndarry?或?pandas.DataFrame?對象作為數據參數傳入。

3、模擬擬合

3.1 MAP方法

使用最優化方法找到點估計，快速簡單。PyMC3提供了?find_MAP?函數，返回參數的一個估計值（point）。

默認情況下，find_MAP?使用Broyden–Fletcher–Goldfarb–Shanno (BFGS) 算法進行最優化運算，找到對數后驗分布的最大值。這里也可以指定?scipy.optimize?模塊中的其他最優化算法完成尋優。

map_estimate = pm.find_MAP(model=basic_model) from scipy import optimize map_estimate2 = pm.find_MAP(model=basic_model,fmin=optimize.fmin_powell) print(map_estimate) print(map_estimate2) Optimization terminated successfully.Current function value: 148.980478Iterations: 6Function evaluations: 269 {'sigma': array(0.96298715), 'sigma_log__': array(-0.03771521), 'alpha': array(0.90661753), 'beta': array([0.94850377, 2.52246124])} {'sigma': array(0.96298742), 'sigma_log__': array(-0.03771493), 'alpha': array(0.90661033), 'beta': array([0.94849417, 2.5224695 ])}

值得注意的是，MAP估計不總是有效的，特別是在極端模型情況下。在高維后驗中，可能出現某一出概率密度異常大，但整體概率很小的情況，在層次化模型中較為常見。

大多數尋找MAP極大值算法都找到局部最優解，那么這種算法這在差異非常大的多模態后驗中會失效。

3.2 MCMC方法

雖然極大后驗估計是一個簡單快速的估計未知參數的辦法，但是沒有對不確定性進行估計是其缺陷。相反，比如MCMC這類基于采樣的方法可以獲得后驗分布接近真實的采樣值。

為了使用MCMC采樣以獲得后驗分布的樣本，PyMC3需要指定和特定MCMC算法相關的步進方法（采樣方法）,如Metropolis, Slice sampling, or the No-U-Turn Sampler (NUTS)。PyMC3中的?step_methods?子模塊提供了如下采樣器：?NUTS,?Metropolis,?Slice,?HamiltonianMC, and?BinaryMetropolis?？梢杂蒔yMC3自動指定也可以手動指定。自動指定是根據模型中的變量類型決定的：?
* 二值變量：指定?BinaryMetropolis?
* 離散變量：指定?Metropolis?
* 連續變量：指定?NUTS

手動指定優先，可覆蓋自動指定。

PyMC3有許多基本采樣算法，如自適應切片采樣、自適應Metropolis-Hastings采樣，但最厲害是的No-U-Turn采樣算法（NUTS），特別是在具有大量連續型參數的模型中。NUTS基于對數概率密度的梯度，利用了高概率密度的區域信息，以獲得比傳統采樣算法更快的收斂速度。PyMC3使用Theano的自動微分運算來進行后驗分布的梯度計算。通過自整定步驟來自適應的設置Hamiltonian Monte Carlo(HMC)算法的可調參數。NUTS不可用于不可微分變量（離散變量）的采樣，但是對于具有不可微分變量的模型中的可微分變量是可以使用的。

NUTS需要一個縮放矩陣作為參數，這個矩陣提供了分布的大致形狀，使NUTS不至于在某些方向上的步長過大而在另外一些方向上的步長過小。設置有效的縮放矩陣有助于獲得高效率的采樣，特別是對于那些有很多未知的（未被觀察的）隨機型隨機變量的模型或者具有高度非正態后驗的模型。不好的縮放矩陣將導致采樣速度很慢甚至完全停止。此外，在多數時候，采樣的初始位置對于高效的采樣也是很關鍵的。

幸運的是，PyMC3使用自動變分推理ADVI (auto-diff variational inference)來初始化NUTS算法，并在?step?參數沒有被指定的情況下會自動指定一個合適的迭代方法（step，采樣器）。

with basic_model:#trace = pm.sample(2000)step = pm.NUTS()#step = pm.Metropolis()trace = pm.sample(2000,step=step) Multiprocess sampling (4 chains in 4 jobs) NUTS: [sigma, beta, alpha] Sampling 4 chains: 100%|██████████| 10000/10000 [00:03<00:00, 2997.90draws/s] print("alpha",trace['alpha'].shape) print(trace['alpha'][0:5],"\n") print("beta",trace['beta'].shape) print(trace['beta'],"\n") print("sigma",trace['sigma'].shape) print(trace['sigma']) alpha (8000,) [1.04278775 0.88323412 0.88755258 0.87381961 0.89942359] beta (8000, 2) [[1.01046211 2.55993421][1.00954145 2.54110235][0.91388241 2.45626273]...[1.16787685 2.60887027][0.7255497 2.42006215][1.00161371 2.75913019]] sigma (8000,) [0.91877597 1.17095769 0.87067159 ... 1.02761247 0.95970655 1.05418962]

可以通過?step?參數指定特定的采樣器（迭代器）來替換默認的迭代器NUTS。這里給?sigma?變量指定?Slice?采樣器（step），那么其他兩個變量（?alpha?和?beta）的采樣器會自動指定（NUTS）。

with basic_model:# 用MAP獲得初始點start = pm.find_MAP(fmin=optimize.fmin_powell)# 實例化采樣器step = pm.Slice(vars=[sigma])# 對后驗分布進行5000次采樣trace = pm.sample(5000, step=step,start=start) /home/fjs/.local/lib/python3.5/site-packages/pymc3/tuning/starting.py:61: UserWarning: find_MAP should not be used to initialize the NUTS sampler, simply call pymc3.sample() and it will automatically initialize NUTS in a better way.warnings.warn('find_MAP should not be used to initialize the NUTS sampler, simply call pymc3.sample() and it will automatically initialize NUTS in a better way.') /home/fjs/.local/lib/python3.5/site-packages/pymc3/tuning/starting.py:102: UserWarning: In future versions, set the optimization algorithm with a string. For example, use `method="L-BFGS-B"` instead of `fmin=sp.optimize.fmin_l_bfgs_b"`.warnings.warn('In future versions, set the optimization algorithm with a string. ' logp = -149: 5%|▌ | 269/5000 [00:00<00:01, 3256.65it/s] Optimization terminated successfully.Current function value: 148.980478Iterations: 6Function evaluations: 269 Multiprocess sampling (4 chains in 4 jobs) CompoundStep >Slice: [sigma] >NUTS: [beta, alpha] Sampling 4 chains: 100%|██████████| 22000/22000 [00:08<00:00, 2581.19draws/s]

4、后驗分布分析

PyMC3提供了?traceplot?函數來繪制后驗采樣的趨勢圖。

pm.traceplot(trace);

左側是每個隨機變量的邊際后驗的直方圖，使用核密度估計進行了平滑處理。右側是馬爾可夫鏈采樣值按順序繪制。對于向量參數?beta會有兩條后驗分布直方圖和后驗采樣值。

同時也提供蚊子形式的后驗統計總結，使用summary函數獲得。

pm.summary(trace)

前面設置的參數如下

alpha=1 beta =[1, 2.5] sigma=1

這里通過后驗估計可以看出平均估計值：

mean(alpha)=[0.906] mean(beta)=[0.948,2.522] mean(sigma)=[0.991] 《新程序員》：云原生和全面數字化實踐50位技術專家共同創作，文字、視頻、音頻交互閱讀

總結

以上是生活随笔為你收集整理的概率编程库Pymc3案例之线性回归的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。