R語言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型|附代碼數(shù)據(jù)
全文下載鏈接：http://tecdat.cn/?p=11974
最近我們被客戶要求撰寫關(guān)于WinBUGS 的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。 R2WinBUGS軟件包提供了從R調(diào)用WinBUGS的便捷功能。它自動(dòng)以WinBUGS可讀的格式寫入數(shù)據(jù)和腳本，以進(jìn)行批處理（自1.4版開始）。WinBUGS流程完成后，可以通過程序包本身將結(jié)果數(shù)據(jù)讀取到R中（這提供了推斷和收斂診斷的緊湊圖形摘要），也可以使用coda程序包的功能對(duì)輸出進(jìn)行進(jìn)一步分析。?
WinBUGS軟件可從http://www.mrc-bsu.cam.ac.uk/bugs/免費(fèi)獲得。??
R是一種“用于數(shù)據(jù)分析和圖形處理的語言”，是一種實(shí)現(xiàn)該語言的開放源代碼和免費(fèi)提供的統(tǒng)計(jì)軟件包，請(qǐng)參見http://www.R-project.org/。? R和R2WinBUGS可從CRAN 獲得，即http://CRAN.R-Project.org或其鏡像之一。?如果可以使用Internet連接，則可以在R命令提示符下鍵入install.packages（“ R2WinBUGS”）來安裝R2WinBUGS。別忘了用library(R2WinBUGS)
例子
學(xué)校數(shù)據(jù)
學(xué)術(shù)能力測(cè)驗(yàn)（SAT）衡量高中生的能力，來幫助大學(xué)做出入學(xué)決定。我們的數(shù)據(jù)來自1970年代后期進(jìn)行的一項(xiàng)實(shí)驗(yàn)，來自八所高中的SAT-V（學(xué)業(yè)能力測(cè)試語言）。SAT-V是由教育測(cè)試服務(wù)局管理的標(biāo)準(zhǔn)多項(xiàng)選擇測(cè)試。該服務(wù)對(duì)所選學(xué)校中每所學(xué)校的教練計(jì)劃的效果很感興趣。
相關(guān)視頻
拓端
拓端
實(shí)現(xiàn)
R2WinBUGS軟件包的實(shí)現(xiàn)非常簡(jiǎn)單。main“函數(shù)bugs（） 由用戶調(diào)用。原則上，它是對(duì)?其中逐步調(diào)用的其他幾個(gè)函數(shù)的包，如下：
bugs.data.inits（）寫入數(shù)據(jù)文件' data.txt”和“ inits1.txt”，“ inits2.txt” ...進(jìn)入 工作目錄。
bugs.script（）寫入WinBUGS用于批處理的文件“ script.txt”。
bugs.run（）更新WinBUGS注冊(cè)表 ，調(diào)用WinBUGS，并使用 'script.txt' 以批處理模式運(yùn)行它。
bugs.sims（）如果參數(shù)codaPkg已設(shè)置為false（默認(rèn)值）才調(diào)用。
否則，bugs（）返回存儲(chǔ)數(shù)據(jù)的文件名。例如，這些可以通過打包的coda 導(dǎo)入，該軟件包提供了收斂診斷，蒙特卡洛估計(jì)的計(jì)算，跡線圖等功能。
bugs.sims（）函數(shù)將WinBUGS中的模擬讀取到R中，將其格式化，監(jiān)視收斂，執(zhí)行收斂檢查并計(jì)算中位數(shù)和分位數(shù)。它還為bugs（）本身準(zhǔn)備輸出。
這些功能不由用戶直接調(diào)用。參數(shù)將從bugs（）傳遞給其他函數(shù)。
例子?
我們將 R2WinBUGS提供的功能應(yīng)用于示例數(shù)據(jù)并分析輸出。
學(xué)校數(shù)據(jù)
示例數(shù)據(jù) ：
>?schools
為了對(duì)這些數(shù)據(jù)進(jìn)行建模，我們使用了Gelman等人提出的分層模型。我們假設(shè)每所學(xué)校的觀測(cè)估計(jì)值具有正態(tài)分布，且均值theta 和方差tau.y，逆方差為1 =σ.y2，其先驗(yàn)分布在（0,1000）上是均勻的。對(duì)于均值theta，我們采用另一個(gè)正態(tài)分布 平均為mu.theta和逆方差為tau.theta。有關(guān)其先驗(yàn)分布，請(qǐng)參見以下WinBUGS代碼：
model {for (j in 1:J){y[j] ~ dnorm (theta[j], tau.y[j])theta[j] ~ dnorm (mu.theta, tau.theta)tau.y[j] <- pow(sigma.y[j], -2)}mu.theta ~ dnorm (0.0, 1.0E-6)tau.theta <- pow(sigma.theta, -2)sigma.theta ~ dunif (0, 1000)}
點(diǎn)擊標(biāo)題查閱往期內(nèi)容
R語言用貝葉斯線性回歸、貝葉斯模型平均 (BMA)來預(yù)測(cè)工人工資
左右滑動(dòng)查看更多
01
02
03
04
此模型必須存儲(chǔ)在單獨(dú)的文件中，例如'schools.bug'2，在適當(dāng)?shù)哪夸浿?，例如c：/ schools /。在R中，用戶必須準(zhǔn)備bugs（）函數(shù)所需的數(shù)據(jù)輸入。這可以是包含每個(gè)數(shù)據(jù)向量名稱的列表，例如
>?J?<- nrow(schools)
使用這些數(shù)據(jù)和模型文件，我們可以運(yùn)行MCMC模擬以獲取theta， mu.theta和sigma.theta的估計(jì)值。在運(yùn)行之前，用戶必須確定要運(yùn)行多少個(gè)鏈 （n.chain = 3）和迭代次數(shù)（n.iter = 1000）。另外，用戶必須指定鏈的初始值，例如通過編寫函數(shù)：
> inits <- function(){+ list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),+ sigma.theta = runif(1, 0, 100))+ }
可以開始MCMC模擬，R中的參數(shù)bugs.directory必須指向WinBUGS的安裝目錄。可以通過print（schools.sim）方便地輸出school.sim對(duì)象中的結(jié)果。?
對(duì)于此示例，將獲得類似的結(jié)果
Inference for Bugs model at "c:/schools/schools.bug"3 chains, each with 1000 iterations (first 500 discarded)n.sims = 1500 iterations savedmean sd 2.5% 25% 50% 75% 97.5% Rhat n.efftheta[1] 11.1 9.1 -3.0 5.0 10.0 16.0 31.8 1.1 39theta[2] 7.6 6.6 -4.7 3.3 7.8 11.6 21.1 1.1 42theta[3] 5.7 8.4 -12.5 0.6 6.1 10.8 21.8 1.0 150theta[4] 7.1 7.0 -6.6 2.7 7.2 11.5 21.0 1.1 42theta[5] 5.1 6.8 -9.5 0.7 5.2 9.7 18.1 1.0 83theta[6] 5.7 7.3 -9.7 1.0 6.2 10.2 20.0 1.0 56theta[7] 10.4 7.3 -2.1 5.3 9.8 15.3 25.5 1.1 27theta[8] 8.3 8.4 -6.6 2.8 8.1 12.7 26.2 1.0 64mu.theta 7.6 5.9 -3.0 3.7 8.0 11.0 19.5 1.1 35sigma.theta 6.7 5.6 0.3 2.8 5.1 9.2 21.2 1.1 46deviance 60.8 2.5 57.0 59.1 60.2 62.1 66.6 1.0 170pD = 3 and DIC = 63.8 (using the rule, pD = var(deviance)/2)For each parameter, n.eff is a crude measure of effective sample size,and Rhat is the potential scale reduction factor (at convergence, Rhat=1).DIC is an estimate of expected predictive error (lower deviance is better).
此外，用戶可以通過輸入plot（schools.sim）生成結(jié)果圖。結(jié)果圖如圖所示。在該圖中，左列顯示了以下內(nèi)容的快速摘要：
推論和收斂（所有參數(shù)的Rb都接近1.0，表明三個(gè)鏈的良好混合，因此近似收斂）；右列顯示每組參數(shù)的推論。從右欄中可以看到，R2WinBUGS使用 WinBUGS中的參數(shù)名稱將輸出構(gòu)造為標(biāo)量，向量和參數(shù)數(shù)組。
點(diǎn)擊文末?“閱讀原文”
本文選自《R語言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型》。
點(diǎn)擊標(biāo)題查閱往期內(nèi)容
R語言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析
python貝葉斯隨機(jī)過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實(shí)現(xiàn)
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機(jī)波動(dòng)率模型、序列蒙特卡羅SMC、M H采樣分析時(shí)間序列
R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機(jī)波動(dòng)率SV模型、粒子濾波、Metropolis Hasting采樣時(shí)間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實(shí)例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機(jī)森林算法預(yù)測(cè)心臟病
R語言中貝葉斯網(wǎng)絡(luò)（BN）、動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯(cuò)頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負(fù)擔(dān)能力數(shù)據(jù)集
R語言實(shí)現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測(cè)選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測(cè)模型
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進(jìn)行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
WinBUGS對(duì)多元隨機(jī)波動(dòng)率模型：貝葉斯估計(jì)與模型比較
R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)
