R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法
原文鏈接：http://tecdat.cn/?p=23236?
原文出處：拓端數(shù)據(jù)部落公眾號
最近我們被客戶要求撰寫關(guān)于貝葉斯MCMC的研究報告，包括一些圖形和統(tǒng)計輸出。
什么是頻率學(xué)派？
在頻率學(xué)派中，觀察樣本是隨機(jī)的，而參數(shù)是固定的、未知的數(shù)量。
概率被解釋為一個隨機(jī)過程的許多觀測的預(yù)期頻率。
有一種想法是 "真實(shí)的"，例如，在預(yù)測魚的生活環(huán)境時，鹽度和溫度之間的相互作用有一個回歸系數(shù)？
什么是貝葉斯學(xué)派？
在貝葉斯方法中，概率被解釋為對信念的主觀衡量。
所有的變量--因變量、參數(shù)和假設(shè)都是隨機(jī)變量。我們用數(shù)據(jù)來確定一個估計的確定性（可信度）。
這種鹽度X溫度的相互作用反映的不是絕對的，而是我們對魚的生活環(huán)境所了解的東西（本質(zhì)上是草率的）。
目標(biāo)
頻率學(xué)派
保證正確的誤差概率，同時考慮到抽樣、樣本大小和模型。
缺點(diǎn)：需要對置信區(qū)間、第一類和第二類錯誤進(jìn)行復(fù)雜的解釋。
優(yōu)點(diǎn)：更具有內(nèi)在的 "客觀性 "和邏輯上的一致性。
貝葉斯學(xué)派
分析更多的信息能在多大程度上提高我們對一個系統(tǒng)的認(rèn)識。
缺點(diǎn)：這都是關(guān)于信仰的問題! ...有重大影響。
優(yōu)點(diǎn): 更直觀的解釋和實(shí)施，例如，這是這個假設(shè)的概率，這是這個參數(shù)等于這個值的概率?？赡芨咏谌祟愖匀坏亟忉屖澜绲姆绞?。
實(shí)際應(yīng)用中：為什么用貝葉斯
具有有限數(shù)據(jù)的復(fù)雜模型，例如層次模型，其中
實(shí)際的先驗(yàn)知識非常少
貝葉斯法則：
一些典型的貝葉斯速記法。
注意:
貝葉斯的最大問題在于確定先驗(yàn)分布。先驗(yàn)應(yīng)該是什么？它有什么影響？
目標(biāo):
計算參數(shù)的后驗(yàn)分布：π（θ|X）。
點(diǎn)估計是后驗(yàn)的平均值。
一個可信的區(qū)間是
你可以把它解釋為一個參數(shù)在這個區(qū)間內(nèi)的概率 。
計算
皮埃爾-西蒙-拉普拉斯（1749-1827）（見：Sharon Bertsch McGrayne: The Theory That Would Not Die)
有些問題是可分析的，例如二項(xiàng)式似然-貝塔先驗(yàn)。
但如果你有很多參數(shù)，這是不可能完成的操作
如果你有幾個參數(shù)，而且是奇數(shù)分布，你可以用數(shù)值乘以/整合先驗(yàn)和似然（又稱網(wǎng)格近似）。
盡管該理論可以追溯到1700年，甚至它對推理的解釋也可以追溯到19世紀(jì)初，但它一直難以更廣泛地實(shí)施，直到馬爾科夫鏈蒙特卡洛技術(shù)的發(fā)展。
MCMC
MCMC的思想是對參數(shù)值θi進(jìn)行 "抽樣"。
回顧一下，馬爾科夫鏈?zhǔn)且粋€隨機(jī)過程，它只取決于它的前一個狀態(tài)，而且（如果是遍歷的），會生成一個平穩(wěn)的分布。
技巧 "是找到漸進(jìn)地接近正確分布的抽樣規(guī)則（MCMC算法）。
有幾種這樣的（相關(guān)）算法。
Metropolis-Hastings抽樣
Gibbs 抽樣
No U-Turn Sampling (NUTS)
Reversible Jump
一個不斷發(fā)展的文獻(xiàn)和工作體系!
Metropolis-Hastings 算法
開始:
跳到一個新的候選位置:
計算后驗(yàn):
如果
如果
轉(zhuǎn)到第2步
Metropolis-Hastings: 硬幣例子
你拋出了5個正面。你對θ的最初 "猜測 "是
MCMC:
r
p.old <- prior *likelihood while(length(thetas) <= n){ ?theta.new <- theta + rnorm(1,0,0.05) ?p.new <- prior *likelihood ?if(p.new > p.old | runif(1) < p.new/p.old){ ? ?theta <- theta.new ? ?p.old <- p.new ?}
畫圖:
r
hist(thetas[-(1:100)] )curve(6*x^5 )
采樣鏈：調(diào)整、細(xì)化、多鏈
那個 "朝向 "平穩(wěn)的初始過渡被稱為 "預(yù)燒期"，必須加以修整。
怎么做？用眼睛看
采樣過程（顯然）是自相關(guān)的。
如何做？通常是用眼看，用acf()作為指導(dǎo)。
為了保證你收斂到正確的分布，你通常會從不同的位置獲得多條鏈（例如4條）。
有效樣本量
MCMC 診斷法
R軟件包幫助分析MCMC鏈。一個例子是線性回歸的貝葉斯擬合（α,β,σ
r
plot(line)
預(yù)燒部分:
r
plot(line[[1]], start=10)
MCMC診斷法
查看后驗(yàn)分布（同時評估收斂性）。
r
density(line)
參數(shù)之間的關(guān)聯(lián)性，以及鏈內(nèi)的自相關(guān)關(guān)系
r
levelplot(line[[2]])acfplot(line)
統(tǒng)計摘要
運(yùn)行MCMC的工具（在R內(nèi)部）
邏輯Logistic回歸：嬰兒出生體重低
r
logitmcmc(low~age+as.factor(race)+smoke )
r
plot(mcmc)
MCMC與GLM邏輯回歸的比較
MCMC與GLM邏輯回歸的比較
對于這個應(yīng)用，沒有很好的理由使用貝葉斯建模，除非--你是 "貝葉斯主義者"。 你有關(guān)于回歸系數(shù)的真正先驗(yàn)信息（這基本上是不太可能的）。
一個主要的缺點(diǎn)是 先驗(yàn)分布棘手的調(diào)整參數(shù)。
但是，MCMC可以擬合的一些更復(fù)雜的模型（例如，層次的logit MCMChlogit）。
Metropolis-Hastings
Metropolis-Hastings很好，很簡單，很普遍。但是對循環(huán)次數(shù)很敏感。而且可能太慢，因?yàn)樗罱K會拒絕大量的循環(huán)。
Gibbs 采樣
在Gibbs吉布斯抽樣中，你不是用適當(dāng)?shù)母怕式邮?拒絕，而是用適當(dāng)?shù)臈l件概率在參數(shù)空間中行進(jìn)。 并從該分布中抽取一次。
然后你從新的條件分布中抽取下一個參數(shù)。
比Metropolis-Hastings快得多。有效樣本量要高得多!
BUGS（OpenBUGS，WinBUGS）是使用吉布斯采樣器的貝葉斯推理。
JAGS是 "吉布斯采樣器"
其他采樣器
漢密爾頓蒙特卡洛（HMC）--是一種梯度的Metropolis-Hastings，因此速度更快，對參數(shù)之間的關(guān)聯(lián)性更好。
No-U Turn Sampler（NUTS）--由于不需要固定的長度，它的速度更快。這是STAN使用的方法（見http://arxiv.org/pdf/1111.4246v1.pdf）。
(Hoffman and Gelman 2011)
其他工具
你可能想創(chuàng)建你自己的模型，使用貝葉斯MC進(jìn)行擬合，而不是依賴現(xiàn)有的模型。為此，有幾個工具可以選擇。
BUGS / WinBUGS / OpenBUGS (Bayesian inference Using Gibbs Sampling) - 貝葉斯抽樣工具的鼻祖（自1989年起）。WinBUGS是專有的。OpenBUGS的支持率很低。
JAGS（Just Another Gibbs Sampler）接受一個用類似于R語言的語法編寫的模型字符串，并使用吉布斯抽樣從這個模型中編譯和生成MCMC樣本?？梢栽赗中使用rjags包。
Stan（以Stanislaw Ulam命名）是一個類似于JAGS的相當(dāng)新的程序--速度更快，更強(qiáng)大，發(fā)展迅速。從偽R/C語法生成C++代碼。安裝：http://mc-stan.org/rstan.html**
Laplace’s Demon?所有的貝葉斯工具都在R中：?http://www.bayesian-inference.com/software
STAN
要用STAN擬合一個模型，步驟是:
為模型生成一個STAN語法偽代碼（在JAGS和BUGS中相同
運(yùn)行一個R命令，用C++語言編譯該模型
使用生成的函數(shù)來擬合你的數(shù)據(jù)
STAN示例--線性回歸
STAN代碼是R（例如，具有分布函數(shù)）和C（即你必須聲明你的變量）之間的一種混合。每個模型定義都有三個塊。
1.數(shù)據(jù)塊:
?int n; // ?vector[n] y; // Y 向量
這指定了你要輸入的原始數(shù)據(jù)。在本例中，只有Y和X，它們都是長度為n的（數(shù)字）向量，是一個不能小于0的整數(shù)。
2. 參數(shù)塊
?real beta1; ?// slope
這些列出了你要估計的參數(shù)：截距、斜率和方差。
3. 模型塊
? ?sigma ~ inv_gamma(0.001, 0.001); ? ?yhat[i] <- beta0 + beta1 * (x[i] - mean(x));} ? ?y ~ normal(yhat, sigma);
注意:
你可以矢量化，但循環(huán)也同樣快
有許多分布（和 "平均值 "等函數(shù)）可用
請經(jīng)常參閱手冊！?https://github.com/stan-dev/stan/releases/download/v2.9.0/stan-reference-2.9.0.pdf
2. 在R中編譯模型
你把你的模型保存在一個單獨(dú)的文件中， 然后用stan_model()命令編譯這個模型。
這個命令是把你描述的模型，用C++編碼和編譯一個NUTS采樣器。相信我，自己編寫C++代碼是一件非常非常痛苦的事情（如果沒有很多經(jīng)驗(yàn)的話），而且它保證比R中的同等代碼快得多。
注意：這一步可能會很慢。
3. 在R中運(yùn)行該模型
這里的關(guān)鍵函數(shù)是sampling()。還要注意的是，為了給你的模型提供數(shù)據(jù)，它必須是列表的形式
模擬一些數(shù)據(jù)。
r
X <- runif(100,0,20)Y <- rnorm(100, beta0+beta1*X, sigma)
進(jìn)行取樣!
r
sampling(stan, Data)
這里有大量的輸出，因?yàn)樗嬎懔?/p>
r
print(fit, digits = 2)
MCMC診斷法
為了應(yīng)用coda系列的診斷工具，你需要從STAN擬合對象中提取鏈，并將其重新創(chuàng)建為mcmc.list。
r
extract(stan.fitalply(chains, 2, mcmc)
最受歡迎的見解
1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)
2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真
4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
7.R語言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析
8.R語言隨機(jī)搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
