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

對于統計計算而言，獲得高維分布樣本后可以用于計算高維空間的積分。對于統計模型而言，獲得高維分布樣本后可以用于估計參數。網上大部分講解理論后給出的是一個估計

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

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

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

1.用Gibbs抽樣

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

假設模型為，

似然函數，

假設三個未知參數的先驗分布，

Gibbs抽樣需要獲得這三個參數的后驗分布，

上面三個分布也就是Gibbs抽樣中滿條件分布，依次循環抽樣至平穩，即為3個參數的分布。

下面推導這些滿條件分布，

表示所有樣本的聯合分布。

一個問題是上式左右是正比連接的，而非等號，貌似沒法求。實際上求一個分布時，我們并不需要得到完整密度函數，只需要一些項即可，如某種正態分布

，只需知道 的值即可獲得該分布的期望方差（和正態分布形式對比即可知），即得到該分布。左右再加個對數，就只需要知道右邊 的系數即可。

更詳細可看貝葉斯估計共軛先驗分布和分布的核的概念。

對右邊的式子取對數（此時

就是似然函數），取出我們關心的含 的項，得到，

由上式得到

的系數為 ， 的系數為 ，由這個兩個系數得到 后驗分布的均值方差，得到后驗分布，

先驗分布是正態分布，關于其后驗分布是否仍是正態，看分布的核的概念，實際就是看看概率函數乘后是否仍是正態分布形式。

同理有，

對右邊取對數，拿出僅和

相關的項，

獲得

的系數為 ， 的系數為 ，得到后驗分布的期望方差，從而有后驗分布為，

取對數，取出僅含

的項，

注意

服從gamma分布，gamma分布取對數為 ，所以這里取 的系數為 ， 的系數為 ，根據這兩個系數得到gamma分布后驗分布為，

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

#參數設定 N = 31mu0 = 0 tau0 = 1mu1 = 0 tau1 = 1alpha = 2 beta = 1#后驗分布抽樣函數 #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) #三個參數初始值設為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])

和最小二乘的結果比較看看，

> #看看最小二乘的結果 > model = lm(y~x) > mean(model$residuals^2) [1] 51.31251 > modelCall: lm(formula = y ~ x)Coefficients: (Intercept) x 6.990 4.945 > > #看看Gibbs的結果 > 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的估計也還行 [1] 8.579605

真實

為3，最小二乘得到7，Gibbs得到2.1，都比較遠，當然這又樣本隨機的原因。 都比較接近真實。

2.Metropolis-Hasting 算法

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

真實模型為

#模擬數據 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")

訓練數據

，假設模型為 ，且已知 。需要訓練得到的參數 。

貝葉斯估計中假設真實參數不是固定的常數，而是服從某種分布。這里假設各個參數的先驗分布為，

。加上樣本信息的后驗分布為，

正比符號去掉了無關的

#樣本的似然函數 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) }# 參數的先驗分布似然 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) }#聯合后驗分布似然 posterior <- function(param){return (likelihood(param) + prior(param)) }######## Metropolis 算法 ################ proposalfunction <- function(param){ #建議密度函數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)) #按列分開了參數chain[1,] = startvaluefor (i in 1:iterations){proposal = proposalfunction(chain[i,])probab = exp(posterior(proposal) - posterior(chain[i,])) #前面取了對數，這里取指數if (runif(1) < probab){ #使用 mcmc 接受-拒絕樣本，獲得(beta_0,beta_1,sigma)的多維聯合樣本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()函數

使用貝葉斯估計的一個主要問題是，

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

在R的MCMCpack包里面，只要寫出似然函數+先驗即可。

仍然以一般回歸為例，對數似然為

與 的無關項去掉了假設了 ，當然 也可以做個先驗分布，這里為方便假設已知。

先驗分布，假設各個參數相互獨立，先驗分布均用正態分布，

，對數為 ，每個 都是正態分布。當然也可以用個3維正態分布做先驗。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) #對數似然prior = sum(log(sapply(beta,dnorm))) #給個先驗分布的對數，程序會自動執行mcmcloglike + prior }m = MCMCmetrop1R(log_fun, theta.init=c(0,0,0),x=x_data, y=y_data,mcmc=4000, burnin=500) #參數x,y和log_fun中的x,y對應#模型診斷 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

和真實還是很接近的。且參數接近正態分布，這是因為正態樣本關于均值的共軛先驗分布為正態分布，故后驗分布為正態分布，即上圖。

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。