参数估计_MCMC-模型参数估计

MCMC方法的目的是获得服从高维分布的样本，理论涉及平稳分布马尔科夫链转移概率等，还是比较麻烦且不好懂的，但好在网上已有不少讲解得比较详细的。

对于统计计算而言，获得高维分布样本后可以用于计算高维空间的积分。对于统计模型而言，获得高维分布样本后可以用于估计参数。网上大部分讲解理论后给出的是一个估计

参数的例子，对于如何具体用到模型中如简单的回归却还是模糊。

这里分别使用MCMC中的Gibbs抽样、Metropolis-Hasting算法对简单回归模型

中的参数

进行参数估计。给出MCMC用于模型（贝叶斯估计）的一个例子，其它复杂模型使用MCMC估计参数时可类似该过程使用。最后给出R语言中MCMCpack包使用mcmc进行参数估计。

1.用Gibbs抽样

资料来源 Gibbs sampling for Bayesian linear regression in Python，下面的代码会改为R的。

假设模型为，

似然函数，

假设三个未知参数的先验分布，

Gibbs抽样需要获得这三个参数的后验分布，

上面三个分布也就是Gibbs抽样中满条件分布，依次循环抽样至平稳，即为3个参数的分布。

下面推导这些满条件分布，

表示所有样本的联合分布。

一个问题是上式左右是正比连接的，而非等号，貌似没法求。实际上求一个分布时，我们并不需要得到完整密度函数，只需要一些项即可，如某种正态分布

，只需知道

的值即可获得该分布的期望方差（和正态分布形式对比即可知），即得到该分布。左右再加个对数，就只需要知道右边

的系数即可。

更详细可看贝叶斯估计共轭先验分布和分布的核的概念。

对右边的式子取对数（此时

就是似然函数），取出我们关心的含

的项，得到，

由上式得到

的系数为

，

的系数为

，由这个两个系数得到

后验分布的均值方差，得到后验分布，

先验分布是正态分布，关于其后验分布是否仍是正态，看分布的核的概念，实际就是看看概率函数乘后是否仍是正态分布形式。

同理有，

对右边取对数，拿出仅和

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

和真实还是很接近的。且参数接近正态分布，这是因为正态样本关于均值的共轭先验分布为正态分布，故后验分布为正态分布，即上图。