- 貝葉斯簡介
先驗：體現了對參數的了解，設定參數的分布
似然：引入觀測數據，反映在給定參數下，觀測數據的可信度
后驗：基於觀測數據，獲得參數的分布和可信度
貝葉斯定理：后驗正比於似然乘以先驗
- 拋硬幣案例
以拋硬幣為例
參數為硬幣朝上的概率A（表示拋一次硬幣時正面朝上的可能性），定義為P(A)，此為先驗（此處指的是概率P的概率）
似然為已知該參數A，觀測值B（N次拋硬幣正面朝上的次數）對應的概率，也就是已知硬幣朝上的概率，實驗中出現這個觀測值的概率
后驗則為通過一系列觀測值B，獲得參數A的分布
首先選擇似然分布，此處選擇了二項分布，表示N次拋硬幣實驗y次正面朝上的概率分布
（似然是指的觀測值的可能分布）
（此時二項分布表達式內既包含了參數A，也包含了觀測值B，也包含了試驗次數N）
實際上就是已知硬幣朝上的概率P（參數A），實驗了N次，求實驗中發生y次硬幣朝上（參數B）的概率
其次為參數A選擇先驗分布，其只與參數A直接相關
選用了beta分布，因為其輸出為0到1，與P一致
計算后驗，將似然去乘以先驗
由於只考慮了其正比的比例，所以可以把與A和B不相關的系數都去掉，得到的后驗分布表達式和先驗分布很相似
也是一個beta分布，此時就可以使先驗beta分布的系數和后驗beta分布的系數構成一個等式，后驗的分布就可以使用先驗的系數
- 注意點
貝葉斯分析的結果是后驗分布而不是某個值，描述了根據給定數據和模型得到不同數值的可能性
后驗最可能的值是根據后驗分布的形態決定的（也就是后驗分布的峰值）
后驗分布的離散程度與我們對參數的不確定性有關，分布越離散，不確定性越大
給定足夠多的數據時，兩個不同先驗的貝葉斯模型會趨近於相同的結果，這個過程也稱之為收斂
不同后驗收斂到相同分布取決於數據和模型
- 先驗的選擇
使用信息量較小的先驗，比如只使用先驗來限定其取值范圍等等，稱為正則化先驗
可以使用多個先驗來驗證先驗的正確性
- 最大后驗密度區間（HPD）
包含一定比例概率密度的最小區間，一般有95%HPD或者98%HPD
如果HPD區間為[2,5]，其含義是參數位於[2,5]的概率是0.95
在貝葉斯學派中，可以告訴我們參數取值的概率，這在頻率學中是不可能的
因為頻率學認為參數是一個固定值，置信區間只能包含或者不包含參數的真實值
def naive_hpd(post): sns.kdeplot(post) HPD = np.percentile(post, [2.5, 97.5]) print(HPD) plt.plot(HPD, [0, 0], label='HPD {:.2f} {:.2f}'.format(*HPD), linewidth=8, color='k') plt.legend(fontsize=16); plt.xlabel(r"$\theta$", fontsize=14) plt.gca().axes.get_yaxis().set_ticks([]) np.random.seed(1) post = stats.beta.rvs(5, 11, size=1000) #根據beta分布（5，11），生成數量為1000的array naive_hpd(post) plt.xlim(0, 1) plt.show()
實際上是找到2.5，97.5概率對應的值即可
- 后驗預測檢查
一旦獲得后驗分布，就可以根據后驗獲得未來的數據來進行預測
通過比較未來數據和已有數據，可以進行一致性檢查
- 拋硬幣代碼實現
theta_real = 0.35 #參數真值 trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150] # 實驗次數 data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48] # 每次實驗的正面朝上的次數 beta_params = [(1, 1), (0.5, 0.5), (20, 20)] # beta先驗分布的參數對 dist = stats.beta x = np.linspace(0, 1, 100) # 0到1之間均勻采樣100個數，用於計算后驗分布的概率密度 print(x) for idx, N in enumerate(trials): # 一次循環表示一次實驗過程 print(idx) print(N) if idx == 0: plt.subplot(4,3, 2) else: plt.subplot(4,3, idx+3) y = data[idx] #每一次實驗對應的實驗數據 print(y) for (a_prior, b_prior), c in zip(beta_params, ('b', 'r', 'g')): # 使用不同的先驗參數，繪制三張圖，分別對應藍色，紅色，綠色 p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y) # 將beta先驗分布參數更換為后驗分布的參數（此處就相當於引入了觀測值y及N），並求解x對應的概率密度，此處已經獲得了后驗分布的參數 print(p_theta_given_y) plt.plot(x, p_theta_given_y, c) plt.fill_between(x, 0, p_theta_given_y, color=c, alpha=0.6) plt.axvline(theta_real, ymax=0.3, color='k') # 用豎直黑線繪制出參數的真值 plt.plot(0, 0, label="{:d} experiments\n{:d} heads".format(N, y), alpha=0) plt.xlim(0,1) plt.ylim(0,12) plt.xlabel(r"$\theta$") plt.legend() plt.gca().axes.get_yaxis().set_visible(False) plt.tight_layout() plt.show()
觀測值B包含了實驗次數（experiments）以及朝上次數（heads）
將值代入后驗分布的參數中，並均勻采樣參數A的可能值導入后驗分布中，即得到了參數A每個可能值的概率密度分布
核心為
p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y)
此處就實現了使用觀測值來更改先驗，使其變成后驗，從而得到參數的可能值分布
（參數的后驗由先驗和觀測值共同組成）
- 自動的后驗
實際上通過分析方法獲得后驗分布在實際中是很難成立的
通過其他非分析解的方法，也可以對后驗進行計算
非馬爾可夫方法：
網格計算，二次近似，變分方法（用簡單的分布來近似后驗）
馬爾科夫方法：
馬爾科夫鏈蒙特卡洛（MCMC）
蒙特卡洛：一種隨機采樣的方法
馬爾科夫鏈：每個狀態轉移到其他狀態的概率只與當前狀態有關，此狀態鏈就成為馬爾科夫鏈
為了找到細節平衡的馬爾科夫鏈，可以使用Metropolis-Hasting算法采樣到一條鏈
該算法擁有自己的方式（即拒絕和接受的標准，采樣是隨機的，但是接受和拒絕的標准是人為的）（接收新的值的概率正比於新舊值之比）來確定采樣方向，從而得到一條采樣軌跡或者采樣鏈
采樣值就是后驗的近似，出現次數最多的值就是后驗中最可能的值
混合蒙特卡洛方法：避免采樣過慢且結果自相關，需要計算梯度
NUTS算法是其中一種HMC算法，不需要手動設置參數
- 使用網格算法硬幣案例
def posterior_grid(grid_points=100, heads=6, tosses=9): """ A grid implementation for the coin-flip problem """ # 參數構成的網格（0和1之前構成的網格） grid = np.linspace(0, 1, grid_points) print(grid) # 定義先驗，將5重復網格數次 prior = np.repeat(5, grid_points) print(prior) #prior = (grid <= 0.4).astype(int) # truncated #prior = abs(grid - 0.5) # "M" prior # 計算每個網格點的似然值,引入可觀測值和網格點 likelihood = stats.binom.pmf(heads, tosses, grid) print(likelihood) # 計算似然乘以先驗，似然相當於就是更新先驗的一個參數因子 unstd_posterior = likelihood * prior print(unstd_posterior) # 標准化后驗，使其值為1 posterior = unstd_posterior / unstd_posterior.sum() print(posterior) return grid, posterior points = 15 h, n = 1, 4 # 四次實驗只有一次正面朝上 grid, posterior = posterior_grid(points, h, n) plt.plot(grid, posterior, 'o-') plt.plot(0, 0, label='heads = {}\ntosses = {}'.format(h, n), alpha=0) plt.xlabel(r'$\theta$', fontsize=14) plt.legend(loc=0, fontsize=14) plt.show()
網格：使用0到1的均布數據
[0. 0.07142857 0.14285714 0.21428571 0.28571429 0.35714286
0.42857143 0.5 0.57142857 0.64285714 0.71428571 0.78571429
0.85714286 0.92857143 1. ]
先驗：定義全為5的均勻分布
[5 5 5 5 5 5 5 5 5 5 5 5 5 5 5]
似然：引入觀測值（此處是引入試驗次數和朝上次數）設定一個似然分布，然后導入網格，即可獲得每個網格值對應的概率分布
[0. 0.22875885 0.35985006 0.41576426 0.41649313 0.37952936
0.31986672 0.25 0.17992503 0.11713869 0.0666389 0.03092461
0.00999584 0.0013536 0. ]
（似然可以看做是更新先驗的一個因素）
后驗：先驗乘以似然來更新先驗獲得后驗
[0. 1.14379425 1.79925031 2.07882132 2.08246564 1.89764681
1.59933361 1.25 0.89962516 0.58569346 0.3331945 0.15462307
0.04997918 0.00676801 0. ]
將后驗歸一化
[0. 0.08239883 0.12961782 0.14975809 0.15002063 0.1367063
0.11521584 0.09004988 0.06480891 0.0421933 0.0240033 0.01113903
0.0036005 0.00048757 0. ]
考慮網格數為15，實驗4次，朝上1次
考慮網格數為1000，實驗200次，朝上85次
- 求圓的面積蒙特卡洛代碼實現
蒙特卡洛主要用於表征隨機性
在變成為2R的正方形內隨機撒N個點
在正方形內畫一個半徑為R的圓，計算在圓內點的個數
（使用關於x，y的不等式來判斷其是否在圓內）
得出pi的估計值
圓內點數除以總點數等於圓面積與正方形面積之比（pi*R*R/4*R*R）
(使用點的個數來代表面積)
N = 10000 x, y = np.random.uniform(-1, 1, size=(2, N)) # 隨機撒點，為中心為原點的正方形 inside = (x**2 + y**2) <= 1 # 計算位於圓內的點的個數 pi = inside.sum()*4/N # 換算出pi error = abs((pi - np.pi)/pi)* 100 outside = np.invert(inside) plt.plot(x[inside], y[inside], 'b.') plt.plot(x[outside], y[outside], 'r.') plt.plot(0, 0, label='$\hat \pi$ = {:4.3f}\nerror = {:4.3f}%'.format(pi, error), alpha=0) plt.axis('square') plt.legend(frameon=True, framealpha=0.9, fontsize=16); plt.show()
獲得pi值
- Metropolis-Hasting算法代碼實現
def metropolis(func, steps=10000): """A very simple Metropolis implementation""" samples = np.zeros(steps) # 指定一維的0矩陣，作為采樣鏈的初始值 old_x = func.mean() # 計算beta分布的平均值，作為采樣的起始位置 old_prob = func.pdf(old_x) # 計算采樣起始值的概率值，應該具有比較大的概率值，用於后續的比較 for i in range(steps): # 一次循環即采樣一次 new_x = old_x + np.random.normal(0, 1) # 以beta分布平均值作為起始值，將正態分布作為更新的偏差值 new_prob = func.pdf(new_x) # 計算采樣值對應的概率 acceptance = new_prob/old_prob # 計算接受概率，即新舊測量值概率之比 if acceptance >= np.random.random(): #隨機取值作為篩選條件，更新采樣值 samples[i] = new_x old_x = new_x old_prob = new_prob else: samples[i] = old_x return samples np.random.seed(345) func = stats.beta(0.4, 2) # 定義了一個指定參數的beta分布 samples = metropolis(func=func) # 對beta分布進行采樣，返回一條采樣鏈 print(samples) x = np.linspace(0.01, .99, 100) y = func.pdf(x) plt.xlim(0, 1) plt.plot(x, y, 'r-', lw=3, label='True distribution') # 真實的參數分布及其對應概率密度函數，也就是最先定義的beta分布 plt.hist(samples, bins=100, density=True, edgecolor='k', label='Estimated distribution') # 根據采樣鏈，獲得其概率密度分布 plt.xlabel('$x$', fontsize=14) plt.ylabel('$pdf(x)$', fontsize=14) plt.legend(fontsize=14) plt.show()
該算法的作用主要是
在給定一個分布（比如此處的beta分布）后，通過采樣獲得一條采樣鏈，該采樣鏈對應的概率密度分布與原分布基本一致
此處就實現了一個由分布到采樣數據之間的轉換過程
- 基於pymc3的代碼實現
np.random.seed(123) n_experiments = 4 theta_real = 0.35 # 生成數據使用的真值，是需要估計的值 data = stats.bernoulli.rvs(p=theta_real, size=n_experiments) # 按照二項分布及真值生成數據，作為了觀測數據 print(data) if __name__ == '__main__': with pm.Model() as our_first_model: # 指定參數的先驗（也就是硬幣朝上的概率） # 按照先驗theta來生成采樣鏈 theta = pm.Beta('theta', alpha=1, beta=1) # 定義似然的分布（也就是已知硬幣朝上的概率時，在實驗中硬幣朝上的概率分布，也就是觀測值的可能性），此處傳入了觀測值 # 似然不會進行采樣，因為其是一個觀測變量 # 似然表示的是觀測值的可能性，其基於參數和觀測值，作為先驗的更新因子來獲得后驗 y = pm.Bernoulli('y', p=theta, observed=data) #y = pm.Binomial('theta',n=n_experimentos, p=theta, observed=sum(datos)) start = pm.find_MAP() step = pm.Metropolis() trace = pm.sample(1000, step=step, start=start) # 從最大后驗值處進行采樣，此處獲得是參數theta的后驗采樣trace burnin = 100 # no burnin # 該參數用於處理老化，因為采樣的開始部分可能效果較差 chain = trace[burnin:] # trace就是獲得的后驗鏈，由之前的采樣鏈計算得到 pm.traceplot(chain) # 每一條采樣鏈的后驗概率密度及其采樣圖像 pm.forestplot(chain) # 每一條采樣鏈的森林圖，表示了其一定的HDI區間 pm.autocorrplot(chain) # 每一條采樣鏈的自相關圖像 pm.plot_posterior(chain) # 后驗的hpd曲線 plt.show() print(pm.summary(chain)) # 獲得后驗的詳細參數
指定參數先驗
指定似然分布
對先驗進行采樣
計算出后驗
每條采樣鏈的森林圖
后驗的所有信息
