MCMC方法的目的是獲得服從高維分布的樣本，理論涉及平穩(wěn)分布馬爾科夫鏈轉(zhuǎn)移概率等，還是比較麻煩且不好懂的，但好在網(wǎng)上已有不少講解得比較詳細(xì)的。

對(duì)于統(tǒng)計(jì)計(jì)算而言，獲得高維分布樣本后可以用于計(jì)算高維空間的積分。對(duì)于統(tǒng)計(jì)模型而言，獲得高維分布樣本后可以用于估計(jì)參數(shù)。網(wǎng)上大部分講解理論后給出的是一個(gè)估計(jì)

參數(shù)的例子，對(duì)于如何具體用到模型中如簡單的回歸卻還是模糊。

這里分別使用MCMC中的Gibbs抽樣、Metropolis-Hasting算法對(duì)簡單回歸模型

中的參數(shù) 進(jìn)行參數(shù)估計(jì)。給出MCMC用于模型（貝葉斯估計(jì)）的一個(gè)例子，其它復(fù)雜模型使用MCMC估計(jì)參數(shù)時(shí)可類似該過程使用。最后給出R語言中MCMCpack包使用mcmc進(jìn)行參數(shù)估計(jì)。

1.用Gibbs抽樣

資料來源 Gibbs sampling for Bayesian linear regression in Python，下面的代碼會(huì)改為R的。

假設(shè)模型為，

似然函數(shù)，

假設(shè)三個(gè)未知參數(shù)的先驗(yàn)分布，

Gibbs抽樣需要獲得這三個(gè)參數(shù)的后驗(yàn)分布，

上面三個(gè)分布也就是Gibbs抽樣中滿條件分布，依次循環(huán)抽樣至平穩(wěn)，即為3個(gè)參數(shù)的分布。

下面推導(dǎo)這些滿條件分布，

表示所有樣本的聯(lián)合分布。

一個(gè)問題是上式左右是正比連接的，而非等號(hào)，貌似沒法求。實(shí)際上求一個(gè)分布時(shí)，我們并不需要得到完整密度函數(shù)，只需要一些項(xiàng)即可，如某種正態(tài)分布

，只需知道 的值即可獲得該分布的期望方差（和正態(tài)分布形式對(duì)比即可知），即得到該分布。左右再加個(gè)對(duì)數(shù)，就只需要知道右邊 的系數(shù)即可。

更詳細(xì)可看貝葉斯估計(jì)共軛先驗(yàn)分布和分布的核的概念。

對(duì)右邊的式子取對(duì)數(shù)（此時(shí)

就是似然函數(shù)），取出我們關(guān)心的含 的項(xiàng)，得到，

由上式得到

的系數(shù)為 ， 的系數(shù)為 ，由這個(gè)兩個(gè)系數(shù)得到 后驗(yàn)分布的均值方差，得到后驗(yàn)分布，

先驗(yàn)分布是正態(tài)分布，關(guān)于其后驗(yàn)分布是否仍是正態(tài)，看分布的核的概念，實(shí)際就是看看概率函數(shù)乘后是否仍是正態(tài)分布形式。

同理有，

對(duì)右邊取對(duì)數(shù)，拿出僅和

相關(guān)的項(xiàng)，

獲得

的系數(shù)為 ， 的系數(shù)為 ，得到后驗(yàn)分布的期望方差，從而有后驗(yàn)分布為，

取對(duì)數(shù)，取出僅含

的項(xiàng)，

注意

服從gamma分布，gamma分布取對(duì)數(shù)為 ，所以這里取 的系數(shù)為 ， 的系數(shù)為 ，根據(jù)這兩個(gè)系數(shù)得到gamma分布后驗(yàn)分布為，

#模擬數(shù)據(jù) x = -15:15 y <- 3 + 5 * x + rnorm(n=length(x),mean=0,sd=7)plot(x,y)

#參數(shù)設(shè)定 N = 31mu0 = 0 tau0 = 1mu1 = 0 tau1 = 1alpha = 2 beta = 1#后驗(yàn)分布抽樣函數(shù) #beta_0 sample_beta_0 = function(beta1,tau){mean = (tau0*mu0 + tau*sum(y - beta1*x))/(tau0+tau*N)sd = sqrt(1/(tau0+tau*N))rnorm(1,mean = mean, sd = sd)}#beta_1 sample_beta_1 = function(beta0,tau){mean = (tau1*mu1+tau*sum((y-beta0)*x))/(tau1+tau*sum(x^2))sd = sqrt(1/(tau1+tau*sum(x^2)))rnorm(1,mean = mean, sd = sd)}#tau sample_tau = function(beta0,beta1){alpha = alpha+N/2beta = beta + sum(((y-beta0-beta1*x)^2)/2)rgamma(1,shape=alpha,rate = beta)}#迭代過程 beta0_r = c(0) #三個(gè)參數(shù)初始值設(shè)為0,0,2 beta1_r = c(0) tau_r = c(2)n = 1e4for(i in 1:n){beta0_r = c(beta0_r,sample_beta_0(beta1_r[i], tau_r[i]))beta1_r = c(beta1_r,sample_beta_1(beta0_r[i+1],tau_r[i]))tau_r = c(tau_r,sample_tau(beta0_r[i+1],beta1_r[i+1]))}tau_r = sqrt(1/tau_r) #sigma^2=1/tau，獲得sigma

畫圖看看，

h = 5000:npar(mfrow=c(2,3))plot(beta0_r[h]) abline(h=3,col="red")plot(beta1_r[h]) abline(h=5,col="red")plot(tau_r[h]) abline(h=7,col="red")hist(beta0_r[h]) hist(beta1_r[h]) hist(tau_r[h])

和最小二乘的結(jié)果比較看看，

> #看看最小二乘的結(jié)果 > model = lm(y~x) > mean(model$residuals^2) [1] 51.31251 > modelCall: lm(formula = y ~ x)Coefficients: (Intercept) x 6.990 4.945 > > #看看Gibbs的結(jié)果 > beta0 = mean(beta0_r[h]) > beta1 = mean(beta1_r[h]) > beta0 [1] 2.126934 > beta1 [1] 4.802838 > > mean((y-beta0-beta1*x)^2) #單從均方誤效果比最小二乘差 [1] 76.57444> mean(tau_r[h]) #tau的估計(jì)也還行 [1] 8.579605

真實(shí)

為3，最小二乘得到7，Gibbs得到2.1，都比較遠(yuǎn)，當(dāng)然這又樣本隨機(jī)的原因。 都比較接近真實(shí)。

2.Metropolis-Hasting 算法

摘自A simple Metropolis-Hastings MCMC in R，這里做下解釋。

真實(shí)模型為

#模擬數(shù)據(jù) trueA <- 5 trueB <- 0 trueSd <- 10 sampleSize <- 31# create independent x-values x <- (-(sampleSize-1)/2):((sampleSize-1)/2) # create dependent values according to ax + b + N(0,sd) y <- trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd)plot(x,y, main="Test Data")

訓(xùn)練數(shù)據(jù)

，假設(shè)模型為 ，且已知 。需要訓(xùn)練得到的參數(shù) 。

貝葉斯估計(jì)中假設(shè)真實(shí)參數(shù)不是固定的常數(shù)，而是服從某種分布。這里假設(shè)各個(gè)參數(shù)的先驗(yàn)分布為，

。加上樣本信息的后驗(yàn)分布為，

正比符號(hào)去掉了無關(guān)的

#樣本的似然函數(shù) likelihood <- function(param){a = param[1]b = param[2]sd = param[3]pred = a*x + bsinglelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)sumll = sum(singlelikelihoods)return(sumll) }# 參數(shù)的先驗(yàn)分布似然 prior <- function(param){a = param[1]b = param[2]sd = param[3]aprior = dunif(a, min=0, max=10, log = T)bprior = dnorm(b, sd = 5, log = T)sdprior = dunif(sd, min=0, max=30, log = T)return(aprior+bprior+sdprior) }#聯(lián)合后驗(yàn)分布似然 posterior <- function(param){return (likelihood(param) + prior(param)) }######## Metropolis 算法 ################ proposalfunction <- function(param){ #建議密度函數(shù)return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3))) }run_metropolis_MCMC <- function(startvalue, iterations){chain = array(dim = c(iterations+1,3)) #按列分開了參數(shù)chain[1,] = startvaluefor (i in 1:iterations){proposal = proposalfunction(chain[i,])probab = exp(posterior(proposal) - posterior(chain[i,])) #前面取了對(duì)數(shù)，這里取指數(shù)if (runif(1) < probab){ #使用 mcmc 接受-拒絕樣本，獲得(beta_0,beta_1,sigma)的多維聯(lián)合樣本chain[i+1,] = proposal}else{chain[i+1,] = chain[i,]}}return(chain) }startvalue = c(4,0,10) chain = run_metropolis_MCMC(startvalue, 10000)burnIn = 5000 acceptance = 1-mean(duplicated(chain[-(1:burnIn),]))par(mfrow = c(2,3))hist(chain[-(1:burnIn),1],nclass=30, , main="Posterior of a", xlab="True value = red line" ) abline(v = mean(chain[-(1:burnIn),1])) abline(v = trueA, col="red" )hist(chain[-(1:burnIn),2],nclass=30, main="Posterior of b", xlab="True value = red line") abline(v = mean(chain[-(1:burnIn),2])) abline(v = trueB, col="red" )hist(chain[-(1:burnIn),3],nclass=30, main="Posterior of sd", xlab="True value = red line") abline(v = mean(chain[-(1:burnIn),3]) ) abline(v = trueSd, col="red" )plot(chain[-(1:burnIn),1], type = "l", xlab="True value = red line" , main = "Chain values of a", ) abline(h = trueA, col="red" )plot(chain[-(1:burnIn),2], type = "l", xlab="True value = red line" , main = "Chain values of b", ) abline(h = trueB, col="red" )plot(chain[-(1:burnIn),3], type = "l", xlab="True value = red line" , main = "Chain values of sd", ) abline(h = trueSd, col="red" )

3.MCMCmetrop1R()函數(shù)

使用貝葉斯估計(jì)的一個(gè)主要問題是，

中分母難求，辦法是當(dāng)作以 為概率求期望，實(shí)際就是積分，那么只要用mcmc方法對(duì) 抽樣，針對(duì) 求均值即可。

在R的MCMCpack包里面，只要寫出似然函數(shù)+先驗(yàn)即可。

仍然以一般回歸為例，對(duì)數(shù)似然為

與 的無關(guān)項(xiàng)去掉了假設(shè)了 ，當(dāng)然 也可以做個(gè)先驗(yàn)分布，這里為方便假設(shè)已知。

先驗(yàn)分布，假設(shè)各個(gè)參數(shù)相互獨(dú)立，先驗(yàn)分布均用正態(tài)分布，

，對(duì)數(shù)為 ，每個(gè) 都是正態(tài)分布。當(dāng)然也可以用個(gè)3維正態(tài)分布做先驗(yàn)。library(MCMCpack)x1 = runif(100) x2 = runif(100)x_data = cbind(1,x1,x2) y_data = 3 + x1 + 5*x2 + rnorm(100)log_fun = function(beta,x,y){dim(beta) = c(3,1)loglike = sum(-(y- x %*% beta)^2) #對(duì)數(shù)似然prior = sum(log(sapply(beta,dnorm))) #給個(gè)先驗(yàn)分布的對(duì)數(shù)，程序會(huì)自動(dòng)執(zhí)行mcmcloglike + prior }m = MCMCmetrop1R(log_fun, theta.init=c(0,0,0),x=x_data, y=y_data,mcmc=4000, burnin=500) #參數(shù)x,y和log_fun中的x,y對(duì)應(yīng)#模型診斷 raftery.diag(m)plot(m)summary(m)Iterations = 501:4500 Thinning interval = 1 Number of chains = 1 Sample size per chain = 4000 1. Empirical mean and standard deviation for each variable,plus standard error of the mean:Mean SD Naive SE Time-series SE [1,] 2.881 0.1965 0.003106 0.01008 [2,] 1.182 0.2418 0.003824 0.01318 [3,] 5.024 0.2425 0.003834 0.013702. Quantiles for each variable:2.5% 25% 50% 75% 97.5% var1 2.4936 2.750 2.871 3.021 3.266 var2 0.7219 1.019 1.176 1.344 1.657 var3 4.5547 4.855 5.023 5.190 5.500

和真實(shí)還是很接近的。且參數(shù)接近正態(tài)分布，這是因?yàn)檎龖B(tài)樣本關(guān)于均值的共軛先驗(yàn)分布為正態(tài)分布，故后驗(yàn)分布為正態(tài)分布，即上圖。

總結(jié)

以上是生活随笔為你收集整理的参数估计_MCMC-模型参数估计的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。