1. 蒙特卡洛估計
若$\theta$是要估計的參數，$y_{1},...,y_{n}$是從分布$p(y_{1},...,y_{n}|\theta) $中采樣的樣本值，假定我們從后驗分布$p(\theta|y_{1},...,y_{n})$中獨立隨機采樣$S$個$\theta$值，則$$ \theta^{(1)},...,\theta^{(S)}\sim^{i.i.d.}p(\theta|y_{1},...,y_{n}) $$
那么我們就能夠通過樣本$\{\theta^{(1)},...,\theta^{(S)}\} $對$p(\theta|y_{1},...,y_{n}) $進行估計，且隨着$S$增加，估計精度也會提高。樣本$\{\theta^{(1)},...,\theta^{(S)}\} $的經驗分布就被稱為$p(\theta|y_{1},...,y_{n}) $的$Monte\;Carlo $估計。
另外，令$g(\theta) $為任意函數，根據大數定理，如果$\theta^{(1)},...,\theta^{(S)} $是從$p(\theta|y_{1},...,y_{n}) $中獨立同分布采樣的，那么有：$$ \frac{1}{S}\sum_{s=1}^Sg(\theta^{(s)})\rightarrow E[g(\theta)|y_{1},...,y_{n}]=\int g(\theta)p(\theta|y_{1},...,y_{n})d\theta \; as\;S\rightarrow \infty $$
當$S\rightarrow \infty $，
| $\bar{\theta}=\sum_{s=1}^S \theta^{(s)}/S\rightarrow E[\theta|y_{1},...,y_{n}] $ $\sum_{s=1}{S}(\theta^{(s)}-\bar{\theta})^2/(S-1)\rightarrow Var[\theta|y_{1},...,y_{n}] $ $\#(\theta^{(s)}\leq c)/S\rightarrow Pr(\theta\leq c|y_{1},...,y_{n}) $ $the\;empirical\;of\;\{\theta^{(1)},...,\theta^{(S)}\}\rightarrow p(\theta|y_{1},...,y_{n})$ $the\; median\; of\; \{\theta^{(1)},...,\theta^{(S)}\}\rightarrow \theta_{1/2} $ $the\;\alpha-percentile\;of\;\{\theta^{(1)},...,\theta^{(S)}\}\rightarrow \theta_{\alpha} $ |
2. Gibbs采樣
$Monto Carlo $方法是通過產生偽隨機數，使之服從一個概率分布$\pi(X) $。當$X$是一維的情況下，這很容易做到，但當$X\in R^k $，往往要么樣本不獨立，要么不符合$\pi$。解決該問題目前最常用的是$MCMC $方法，其樣本的產生與馬氏鏈有關。基於條件分布的迭代采樣是另外一種重要的$MCMC $方法，最著名的就是$Gibbs $采樣。
已知$p(\theta|\sigma^2,y_{1},...,y_{n}) $和$p(\sigma^2|\theta,y_{1},...,y_{n}) $分別是$\theta $和$\sigma^2$的全條件分布，給定當前參數$\phi^{(s)}=\{\theta^{(s)},\tilde{\sigma}^{2(s)}\} $，我們按下列方式生成新參數：
| 1. $sample\;\theta^{(s+1)}\sim p(\theta|\tilde{\sigma}^{2(s)},y_{1},...,y_{n}) $ 2. $sample\;\tilde{\sigma}^{2(s+1)}\sim p(\tilde{\sigma}^2|\theta^{(s+1)},y_{1},...,y_{n}) $ 3. $let\;\phi^{(s+1)}=\{\theta^{(s+1)},\tilde{\sigma}^{2(s+1)}\} $ |
最終生成一個相關馬氏序列$\{\phi^{(1)},\phi^{(2)},...,\phi^{(S)}\} $。
3 Methopolis-Hastings算法
當先驗分布是共軛或是半共軛的時候，后驗分布可以用$Monte\; Carlo $方法或是$Gibbs$采樣的方法進行估計。但當先驗分布不能采用共軛分布的時候，參數的全條件分布就沒有標准形式，因此很難采用$Gibbs$采樣進行估計。$Methopolis-Hastings $算法作為一種一般的后驗概率估計算法，可以用於任何先驗分布和采樣模型。
讓我們考慮一般的情況，我們的采樣模型為$Y\sim p(y|\theta) $，先驗分布為$p(\theta) $，雖然對大多數問題來說，我們可以用$p(\theta|y)=p(\theta)p(y|\theta)/\int p(\theta')p(y|\theta')d\theta' $來計算$p(\theta|y) $，但是通常情況下，積分是很難算的。如果我們能夠從$p(\theta|y) $中采樣$\theta^{(1)},...,\theta^{(S)}\sim^{i.i.d.} p(\theta|y) $，那么就可以采用$Monte\; Carlo $方法來估計后驗，使得$$ E[g(\theta)|y]\approx \frac{1}{S}\sum_{s=1}^S g(\theta^{(s)}) $$
但是如果我們不能夠直接從$p(\theta|y) $中采樣呢？其實在估計后驗分布的過程中，最重要的事情不是從$p(\theta|y) $中獲取獨立同分布的樣本，而是要構造一個很大的數據集使得對於任何兩個不同的參數來說有$$ \frac{\#\{\theta^{(s)}\;in\;the\;collection\;=\;\theta_{a}\}}{\#\{\theta^{(s)}\;in\;the\;collection\;=\;\theta_{b}\}}\approx \frac{p(\theta_{a}|y)}{p(\theta_{b}|y)} $$
假定現在的工作集合是$\{\theta^{(1)},...,\theta^{(s)}\} $，我們考慮是否將新參數$\theta^{(s+1)} $加到該集合當中。基本思路是:$$ r=\frac{p(\theta^*|y)}{p(\theta^{(s)|y})}=\frac{p(y|\theta^*)p(\theta^*)}{p(y)}\frac{p(y)}{p(y|\theta^{(s)})p(\theta^{(s)})}=\frac{p(y|\theta^*)p(\theta^*)}{p(y|\theta^{(s)})p(\theta^{(s)})} $$
如果$r>1 $，因為$\theta^{(s)} $已經在我們的數據集當中，由於$\theta^* $擁有更高的概率，我們應該接受$\theta^* $到我們的數據集當中，即$\theta^{(s+1)}=\theta^* $。
如果$r<1 $，我們以$r$或$1-r$的概率來接受$\theta^* $或$\theta^{(s)} $。
我們在抽取$\theta^* $時要服從對稱分布$J(\theta^*|\theta^{(s)}) $，即$J(\theta_{b}|\theta_{a})=J(\theta_{a}|\theta_{b}) $，通常很好選取，如$J(\theta^*|\theta^{(s)}) = normal(\theta^{(s)},\sigma^2)$
綜上，$Methopolis $算法如下
|
$Methopolis-Hastings $算法如下:
利用該算法估計$p_{0}(u,v) $
| 1. 更新$U$ a) 采樣$u^*\sim J_{u}(u|u^{(s)},v^{(s)}) $ b) 計算接受率$$r=\frac{p_{0}(u^*,v^{(s)})}{p_{0}(u^{(s)},v^{(s)})}\times \frac{J_{u}(u^{(s)}|u^*,v^{(s)})}{J_{u}(u^*|u^{(s)},v^{(s)})} $$ c) 以$\min(1,r) $和$\max(0,1-r) $的概率將$u^{(s+1)} $設為$u^* $和$u^{(s)} $ 2. 更新$V$ a) 采樣$v^*\sim J_{v}(v|u^{(s+1)},v^{(s)}) $ b) 計算接受率$$r=\frac{p_{0}(u^{(s+1)},v^*)}{p_{0}(u^{(s+1)},v^{(s)})}\times \frac{J_{v}(v^{(s)}|u^{(s+1)},v^* )}{J_{v}(v^*|v^{(s+1)},v^{(s)})} $$ c) 以$\min(1,r) $和$\max(0,1-r) $的概率將$v^{(s+1)} $設為$v^* $和$v^{(s)} $ |
$Methopolis $算法是$Methopolis-Hastings $算法的一種特殊形式，$Gibbs $抽樣算法也是$Methopolis-Hastings $算法。
參考文獻：Hoff, Peter D. A first course in Bayesian statistical methods. Springer Science & Business Media, 2009.