---

貝葉斯簡介

貝葉斯學派和頻率學派的核心區別：用信念度而不是長期頻率來衡量某件事發生的概率，這個信念度可能會隨著新信息的加入而發生變化。貝葉斯的三個核心概念：

• 先驗概率（prior）：在沒有進行任何嘗試時一個人對某件事的信念，例如人群中X病毒攜帶率為3%；
• 似然（likelihood）：觀測到的數據，例如A檢測結果為陽性；
• 后驗（posterior）：通過觀測到的數據對某件事的信念進行修正，例如A攜帶X病毒的概率為75%。

貝葉斯法則：根據數據更新你關于模型的概率的公式
$$P(MODEL\mid DATA)\propto P(MODEL)\times P(DATA\mid MODEL)$$

---

Why R

眾多貝葉斯計算引擎都可以在R中直接調用。

貝葉斯計算引擎：
貝葉斯引擎是用于貝葉斯統計建模和推斷的工具，這些引擎能幫助用戶構建復雜的概率模型，并利用馬爾可夫鏈蒙特卡羅（MCMC）等方法進行推斷。

• BUGS
• WinBUGS/OpenBUGS
• JAGS（后續會更新相關文章專門介紹JAGS）
• Stan

用于貝葉斯分析/計算的R包：

1. MCMCpack: 包含針對常用模型的 MCMC 算法，通常在社會和行為科學中使用
2. coda: 包含一套用于收斂診斷和輸出分析的函數，可用于總結、繪圖和診斷 MCMC 樣本的收斂性
3. rjags: 實現在 R 中完全交互式地使用JAGS的 R 包
4. r2jags: 基于 rjags，但提供了更為簡便的接口，簡化了 JAGS 的調用流程
5. runjags: 一個包含許多 JAGS 增強功能的 R 包，包括自定義 JAGS 模塊
6. rstan: Stan 的 R 語言接口。它允許用戶方便地從 R 中擬合 Stan 模型并訪問輸出結果
7. Nimble: 一個提供通用 MCMC 系統的 R 包，允許對用 BUGS/JAGS 模型語言編寫的模型進行自定義 MCMC。用戶可以選擇采樣器并編寫新的采樣器。模型和采樣器通過生成的C++ 代碼自動編譯以提高速度
8. LaplacesDemon: 一個旨在提供完整貝葉斯環境的 R 包，包括眾多 MCMC算法、具有多個優化算法的拉普拉斯近似、大量示例、數十種額外的概率分布以及眾多 MCMC 診斷工具等

---

理論基礎

單參數貝葉斯模型：僅包含一個待估計參數的概率模型的R語言建模實現，一個直觀的例子：總體比例p的貝葉斯推斷，希望從樣本數據中得到p的點估計和區間估計（不確定性）。本文將從總體比例p估計問題出發，探討單參數貝葉斯建模過程。

貝葉斯分析的三步走：

• 指定先驗分布
• 寫出似然函數
• 計算參數的后驗分布

上面三個步驟中，我們需要研究的問題：

1. 先驗分布：如何合理表達先驗信念

• 選擇什么先驗分布方便推導？（無信息先驗、層級先驗、共軛先驗等）

2. 似然函數：如何用分布描述數據生成過程

• 使用什么分布建模？（正態、二項、泊松等）

3. 后驗分布：如何計算后驗分布和后驗預測

• 是否能夠解析計算？（有共軛形式時可以解析；否則需數值方法如 MCMC）

對于總體比例估計問題，我們的似然函數就是二項分布，即獨立重復試驗中成功次數。若我們觀察到y 次成功，在 n 次試驗中：$y～Binomial(n,p)$，其似然函數為：
$$L(p)=P(y\mid p)=(n,y)p^y(1?p{)}^{n?y}$$

貝葉斯模型：
$$p(𝜃,y)=𝜋(𝜃)p(y|𝜃)$$

• $𝜋(𝜃)$是先驗分布
• $p(y|𝜃)$是似然函數
• 貝葉斯模型即𝜃和y的聯合分布$p(𝜃,y)$

我們關心的是參數的后驗分布$𝜋(𝜃\mid y)$，它可以由貝葉斯公式計算出：
$$𝜋(𝜃\mid y)=p(𝜃,y)/p(y)=𝜋(𝜃)p(y\mid 𝜃)/p(y)$$
其中𝑝(𝑦) 不依賴于𝜃，當給定觀測數據時可視作常數，從而可以看作是一個歸一化常數，用于將𝜋(𝜃 ∣ 𝑦)化成有效的概率密度函數。

---

一、三種先驗分布和對應后驗的計算

構建先驗分布的常見方法是假設先驗分布屬于特定的參數密度函數族。然后選擇先驗分布的參數（稱為超參數 hyper parameters），使先驗分布盡可能地表示先驗信念。

在研究總體概率p時，有三種典型先驗：1. 離散先驗；2. Beta先驗（共軛先驗）；3. 直方圖先驗。本文下面將詳細講解這三種先驗以及對應后驗的計算。

1. 離散先驗

假設𝑝只能取有限個值（如 0.1, 0.2, …, 0.9），我們給這些值一個先驗概率分布：$P(p=p_i)=w_i$， $w_i$求和為1。

離散先驗分布：

p = seq(0.05, 0.95, by = 0.1)
prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior = prior/sum(prior)
plot(p, prior, type = "h", ylab="Prior Probability")

• p為概率值向量
• prior為p對應的先驗概率向量，歸一化處理
• type = "h"為直方圖選項，用于繪制概率質量函數（pmf）

輸出：

已知觀測數據后計算后驗分布：

library(LearnBayes)
data = c(11, 16)
post = pdisc(p, prior, data)
round(cbind(p, prior, post),2)

• data為觀測數據：11次成功，16次失敗
• pdisc用于計算離散先驗下的后驗概率，三個參數分別是概率點向量、先驗概率、觀測數據

輸出：

對比先驗和后驗：

library(lattice)
PRIOR=data.frame("prior",p,prior)
POST=data.frame("posterior",p,post)
names(PRIOR)=c("Type","P","Probability")
names(POST)=c("Type","P","Probability")
data=rbind(PRIOR,POST)
xyplot(Probability~P|Type,data=data, layout=c(1,2), type="h", lwd=3,col="black")

• data.frame用于創建一個三列的數據框，names用于給數據框命名
  rbind表示將兩個數據框按行合并
  xyplot是lattice包的核心函數，用于繪制散點圖或線圖，這里Probability~P|Type表示縱軸為Probability,橫軸為P ，按照Type分面板繪制

輸出：

2.Beta先驗（共軛先驗）

從上圖介紹我們可以得到為什么要選擇Beta先驗：它和似然函數（二項分布）有相同形式的密度函數，這樣似然*先驗得到的后驗還是這個形式。

問題：如何確定Beta分布的兩個參數？

• 策略1: 畫出來好多條密度函數，找一條滿足期望的
• 策略2: 根據想要的均值/樣本量確定
  - $Beta(\alpha ,\beta )$分布的均值為$\alpha /(\alpha +\beta )$
  - $Beta(\alpha ,\beta )$分布的等效先驗樣本量為$\alpha +\beta ?2$（其中$\alpha ?1$次成功，$\beta ?1$次失敗）

在R中可以用LearnBayes包中的beta.select()函數找到與先驗知識（對分位數的認知）匹配的兩個超參數。

quantile1=list(p=.5,x=.3) # 中位數
quantile2=list(p=.9,x=.5) # 90%分位數
beta.select(quantile1,quantile2)

輸出：
[1] 3.26 7.19

• 先驗信息：研究者認為比例 𝑝 的中位數為 0.3，90% 分位數為 0.5
• 選出的$\alpha =3.26$，$\beta =7.19$

假設我們選出了一個合適的先驗$Beta(\alpha ,\beta )$，觀測$n$次中有$y$次成功，就可以用貝葉斯公式求出后驗：
$$\pi \propto p^y(1?p{)}^{n?y}{p}^{\alpha ?1}(1?p{)}^{\beta ?1}=p^{\alpha +y?1}(1?p{)}^{\beta +n?y?1}$$

這是$Beta(\alpha +y,\beta +n?y)$的核，這就是共軛性的含義：先驗分布和后驗分布來自同一種分布類型。

3. 直方圖先驗

就是每個區間內指定一個先驗概率：

midpt = seq(0.05, 0.95, by = 0.1) # 按照0.1的間隔劃分0.05到0.95的區間
prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0) # 每個區間的先驗概率
prior = prior/sum(prior) # 歸一化
curve(histprior(x,midpt,prior),from=0, to=1,ylab="Prior density",ylim=c(0,.3))

輸出：

計算后驗分布一樣，乘上似然函數的等價beta分布表示：

curve(histprior(x,midpt,prior) * dbeta(x,s+1,f+1),from=0, to=1,ylab="Posterior density")

輸出：

---

二. 后驗抽樣

后驗抽樣是指在獲取參數后驗分布后，從后驗分布中采樣參數，進而用于后續貝葉斯推斷。

1. 網格點采樣法

網格點采樣法，也稱為網格近似法 (Grid Approximation)，是一種通過在參數空間中定義一個離散的網格，然后計算每個網格點上的未歸一化后驗概率（即似然函數值乘以先驗概率密度值）來近似后驗分布的方法。一旦計算出這些概率，就可以將它們歸一化，使其總和為 1，從而得到一個離散化的后驗概率分布，從這個離散分布中抽取樣本。

p = seq(0, 1, length=500) # 網格點，500個0到1之間均勻的數值點
post = histprior(p, midpt, prior)*dbeta(p, s+1, f+1) # 以直方圖先驗為例，計算每個網格點的后驗概率
post = post/sum(post) # 后驗概率歸一化
ps = sample(p, 5000, replace = TRUE,prob = post) # 從post分布沖抽取5000個樣本
hist(ps, breaks= 30, xlab="p",main="")

輸出：

2. 其他方法

后續會專門寫一篇文章講各種MCMC方法。

---

三、貝葉斯推斷

通過確定先驗分布后得到的后驗分布包含了關于未知參數的所有當前信息，所有貝葉斯推斷都基于后驗分布，包括：

1. 參數估計
   • 點估計：后驗均值 (posterior mean)、眾數 (posterior mode / MAP)、中位數 (posterior median)
   • 離散程度：后驗方差
   • 可信區間(credible interval)：后驗分位數
2. 假設檢驗：后驗概率
3. 預測：估計潛在可觀測但當前未觀測的量
4. 模型選擇/變量選擇：Bayes factor,DIC

1. 參數估計

(1) 后驗均值

后驗分布的均值常用作$p$點估計，$Beta(\alpha +y,\beta +n?y)$的均值為$\alpha +y/\alpha +\beta +n$

(2) 后驗方差

后驗方差是后驗分布離散程度的一種度量，后驗方差越大，我們對參數的不確定性就越大。

$Beta(\alpha ,\beta )$的方差為：$\alpha \beta /(\alpha +\beta {)}^2(\alpha +\beta +1)$

對于 $Beta(\alpha ,\beta )$先驗和二項似然（n 次試驗中有 y 次成功），后驗分布為 $Beta(\alpha +y,\beta +n?y)$，從而后驗方差為$(\alpha +y)(\beta +n?y)/(\alpha +\beta +n{)}^2(\alpha +\beta +n+1)$

注意：對于均勻先驗和二項似然，后驗方差（幾乎）總是小于先驗方差

(3) 后驗區間

兩種類型的可信區間：

1. 等尾可信區間(equal-tail credible interval):

100(1 ? α)% 的等尾可信區間是從后驗分布的 α/2 分位數到 (1 ? α/2) 分位數的區間，若要構建 95% 的等尾可信區間，我們需要計算后驗分布的 0.025 和 0.975 分位數。

例如，求后驗均值的95%置信區間有兩種方法，一是直接根據beta分布來計算：

a = 3.26; b = 7.19 # 選出的先驗參數
s = 11; f = 16 # 觀測數據（似然）
qbeta(c(0.05, 0.95), a + s, b + f)

輸出：[1] 0.2555267 0.5133608

二是用模擬的方法：

ps = rbeta(1000, a + s, b + f)
quantile(ps, c(0.05, 0.95))

輸出：5% 95%
0.2490099 0.5152768

• beta表示從beta分布中生成指定個隨機數，
• 這里生成了1000個樣本，然后計算樣本的0.05和0.95分位數

對于貝葉斯學派，95% 可信區間的解釋是：觀察到觀測數據后，真實參數 p 落在該區間內的概率為 0.95

2. 最高后驗密度區間(highest posterior density, HPD):

區間內任意點的密度大于區間外任意點的密度，在包含指定概率的所有區間中長度最短（有時稱為最小長度可信區間）。

當后驗分布高度偏斜或多峰時，優于等尾可信區間

計算可以使用 TeachingDemos 包中的 R 函數hpd()或emp.hpd()

2. 假設檢驗

假設我們要檢驗關于參數 p 的以下假設：
$$H0:p\le 0.10v.sH1:p>0.10$$
我們只需計算后驗分布在這兩個區間上的概率，使用R函數 pbeta即可。

3. 預測

(1) 先驗預測

以離散先驗為例：

library(LearnBayes)
p = seq(0.05, 0.95, by=.1)
prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.8, 0.3, 0.1, 0.05)
prior = prior/sum(prior)
ns = 20
ys = 0:20
pred = pdiscp(p, prior, ns, ys)
round(cbind(0:20, pred), 3)

• 前面都和第一部分的離散先驗一樣
• ns表示未來樣本量，ys表示未來可能的成功個數（從0～20）
• pdiscp用于計算離散先驗下的后驗預測分布，四個參數分別為概率點向量、先驗概率、未來樣本量、未來成功次數

輸出：

(2) 后驗預測

以Beta先驗為例，有解析計算和模擬兩種方法：

1. 解析計算：本質就是更新beta分布的兩個參數

library(LearnBayes)
ab = c(3.26, 7.19)
ns = 20; ys = 0:20
pred = pbetap(ab, ns, ys)
round(cbind(0:20, pred), 3)

輸出：

• pbetap函數計算 Beta 先驗下的后驗預測分布（Beta-Binomial 分布），三個參數分別是Beta先驗的參數向量、未來試驗次數、需要預測的成功次數

2.基于模擬：模擬的方法適合任何先驗分布，下例展示了beta先驗的后驗模擬

m = 1000
a = 3.26; b = 7.19 # 先驗Beta分布的參數
s = 11; f = 16
ns = 20
p = rbeta(m, a+s, b+f)
y = rbinom(m, ns, p)
freq = table(y) # 統計各成功次數y出現的頻次
ys = as.integer(names(freq)) # 提取成功次數（去重）
predprob = freq/sum(freq) # y的概率分布
# 成功概率分布的直方圖繪制
plot(ys, predprob, type="h", xlab="y",
ylab="Predictive Probability",
cex.lab=1, cex.axis=1.0)
# 解析解驗證
library(LearnBayes)
TruePred1 = pbetap(c(a+s,b+f), ns, ys) # beta-二項后驗分布的概率
points(ys, TruePred1, type="p", col='red')

• rbeta從后驗分布中生成1000個p的隨機樣本
• rbinom對每個采樣的 p，生成 ns=20 次試驗的成功次數 y

輸出：

如果后驗是混合分布，那么后驗預測分布是一個連續混合分布或復合分布，下例中展示了用模擬的方法從混合正態分布0.3N(0, 1) + 0.7N(3, 1/2)中隨機抽樣：

n <- 10000; alpha <- 0.3
k <- sample(1:2, size = n, replace = T, prob=c(alpha, 1-alpha)) # 根據概率分到1組或2組
mu <- numeric(length(k)) 
sd <- numeric(length(k)) 
mu <- c(0,3)[k] # 根據k選擇均值（0或3）
sd <- c(1, 1/2)[k]# 根據k選擇方差（1或1/2）
x <- rnorm(n, mu, sd) # 生成混合樣本
hist(x, breaks='FD', prob=T, xlim=c(-4,6), ylim=c(0,0.8), main='Histogram of a mixture normal')
# 混合密度函數
p <- function(x,alpha){
alpha*dnorm(x, 0, 1)+(1-alpha)*dnorm(x, 3, 1/2)
}
curve(dnorm(x,0,1), col=4, lty=2, lwd=2, add=T) # N（0，1）
curve(dnorm(x,3,1/2), col=4, lty=4, lwd=2, add=T) # N（3，1/2）
curve(p(x,0.3), col=2, lty=1, lwd= 2, add=T) # 混合密度
legend("topleft", legend = c("N(0,1)","N(1,1/2)", "0.3N(0,1)+0.7N(1,1/2)"), cex=0.6, lty = c(2,4,1), col=c(4, 4, 2), lwd =c(2,2,2))

輸出：