基于贝叶斯法的线性模型参数估计

张金龙

jinlongzhang01@gmail.com

在探讨一些统计问题时， 人们有时候希望基于现有数据， 得到某些参数的概率密度分布。以便对参数的估计有更加深入的了解。此外,有时 模型过于相当复杂，精确的公式难以给出解析解， 人们希望用模拟的方法，基于现有的数据， 探讨某些参数的分布。 这里就必须要用到贝叶斯推断。

介绍贝叶斯推断， 则绕不过贝叶斯公式， 这是一个表示参数集theta以及数据y关系的概率公式。

#### 根据贝叶斯公式

p(theta|y) = {p(y|theta)*p(theta)}/p(y)

其中theta为参数集， y为数据：

p(y|theta) 为给定参数下， 数据y出现的概率， 即Likelihood。

p(theta) 为参数集theta的先验分布，可假设为任何概率密度分布。

p（y）则称为normalizing constant， 一般并不关心其取值。

但是， 要获得给定数据下的参数的分布却并不容易，p（y）是极难获取的值。 幸运的是， 在统计学家发明了一种基于随机抽样的方法， 可以获得参数的近似分布， 而不必考虑p（y）。 这几种有两种最常用技术，都是基于蒙特卡罗马尔科夫链。

如果一个过程的参数， 在下一个时刻的状态仅与当前时刻有关， 而与之前的状态无关， 则该过程具有马尔科夫性。换句话说， 给定一个模型， 如果要估计的一系列参数theta1, theta2, theta3, ..., theta0, 当前的取值为 theta1_present, theta2_present, theta3_present, ..., thetaN_present, 下一个时刻的取值为theta1', theta2', theta3', ..., thetaN', 每个参数在未来取值只与对应的参数在当前时刻的取值有关， 这样一个过程， 就是一个马尔科夫过程。 记录下所有时刻参数的变化， 就构成了马尔科夫链。现代统计学中，在求复杂贝叶斯积分时， 用到马尔科夫链以及一些随机抽样的方法， 获得参数的近似分布。 常用的有两种方法， 分别称为Gibbs Sampling和Metrapolis Sampling。

马尔科夫链每变化一次， 称为一代。在控制马尔科夫链的状态变化以及抽样函数设置合理的情况下， 蒙特卡罗马尔科夫链进行的代数越多， 则所得结果越接近最大似然估计的结果。

本文就尝试用后者对简单的线性回归进行参数估计。当然， 笔者作为初学者， 理解上难免出现错误， 还请各位读者指正。

以下是Metrapolis Hastings算法的简要介绍

（1）在参数空间指定一套参数的初始值，作为起始点。初始值可以随机给出。

（2）按照参数的先验概率分布， 在初始值的基础上， 每个参数进行一定量的随机变化，计算当前参数组合的概率密度。

（3）依据当前参数组合和起始点概率密度比值是否大于0到1之间的一个均匀分布的随机数来判断是否保留当前点。

（4）若当前点的概率密度大于该随机数，就称这种状态为接受状态，此时，在满足参数概率分布的前提下，继续随机抽取参数的组合，作为下一点，计算下一点的概率密度，并计算下一点概率密度和概率密度的比值，并继续循环。

（5）若当前点不能被接受，则继续在满足参数概率分布的前提下，继续生成随机数，作为新的参数组合，直到参数组合能够被接受为止。

更详细的理论介绍， 请参阅 http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

以下用实例， 即用MCMC （基于Metropolis Hastings 抽样）估计线性回归模型的参数

### 第一步 首先生成要拟合的数据

set.seed(12345)     ### 设定随机数种子

x = 1:60            ### 生成x数据， 共60个点

beta0_actual = 20   ### beta0真实值为20

beta1_actual = 5    ### beta1真实值为5

sigma_actual = 15   ### sigma真实值为15

y = x*5 + 20 + rnorm(60, 0, 15) ## 生成y数据， 这里假设已知 beta1 = 5， beta0 = 20， mu = 0, sd = 15

plot(y ~ x)         ### 绘图

# 本套模拟数据所应的 beta0, beta1

linmod <- lm(y ~ x)

summary(linmod)

### 第二步， 写出 似然函数， 因线性回归一般模式均可写成

### y = beta1*x + beta0 + error,

### 而error符合 N(0, sigma)

### 所有error的概率乘积， 即联合概率密度， 称为Likelihood

### 因为Likelihood数值往往很小，在计算机中为了便于计算， 一般都取自然对数，

### 因此此时dnorm()之外是相加sum()

LogLikelihood <- function(beta0, beta1, sigma){

  R = y - x*beta1 - beta0  

  sum(dnorm(R, mean = 0, sigma, log = TRUE))

}

#### LogLikelihood(50, 6, 1)

#### 第三步给出先验分布

#### beta0， beta1， sigma等所要估计参数的先验分布， 先验分布与取值以及概率密度与统计分布无关。

#### 这里为了演示方便beta0prior，beta1prior 取自均匀分布。prior三者同时发生， 因此取联合概率密度。

#### min，max的范围，也是根据参数的可能取值范围而估计。

prior <- function(beta0, beta1, sigma){

   beta0prior  = dunif(beta0, min=-100, max=100, log = TRUE)

   beta1prior  = dunif(beta1, min=-100, max=100, log = TRUE)

   sigmaprior  = dnorm(sigma, mean = 0, sd = 10, log = TRUE)

   return(beta0prior + beta1prior + sigmaprior)

}

### prior(5,5,1)

### 写出Proposalfunction， 控制马尔科夫链状态变化的强弱， 即每一次变化， 数值的大小。  Proposalfunction一般取对称分布的概率密度，一般取正态分布，这是决定坐标点（beta0, beta1, sigma）在空间移动快慢的函数， 其中sd的取值， 决定了链的冷热程度。

proposalfunction <- function(beta0, beta1, sigma){

   return(rnorm(3, mean = c(beta0, beta1, sigma), sd= c(0.1, 0.5, 0.3)))

}

### 后验分布函数， 后验分布。 因为概率密度均已经取log， 所以这里概率相乘用加法即可。  

posterior <- function(beta0, beta1, sigma){

  return (LogLikelihood(beta0, beta1, sigma) + prior(beta0, beta1, sigma))

}

MCMC <- function(intialvalue, generations){

   chain = matrix(rep(NA, ((generations+1) * 3)), ncol = 3)  ## 为马尔科夫链准备一个矩阵， 三列， 分别放置beta0, beta1, sigma的值

   colnames(chain) <- c("beta0", "beta1", "sigma")            ## 为列命名

   chain[1,] = intialvalue                                     ## 读取初始值

   for (i in 1:generations){                                  ## 开始循环， 循环次数按照代数而定

       proposal = proposalfunction(chain[i,1], chain[i,2], chain[i,3])  ### 围绕当前状态， 生成一套新的参数， 即让状态变化一次

       probab = exp(posterior(proposal[1], proposal[2], proposal[3]) - posterior(chain[i,1], chain[i,2], chain[i,3]))   ### 计算新状态与旧状态概率的比值， 因为概率之前都已经取对数， 因此这里取减号。

       if (runif(1) < probab){  ## 生成一个0到1之间的随机数，如果该随机数小于新旧状态的比值， 则从记录下新的状态， 作为下一个时刻的起始状态。

           chain[i+1,] = proposal

       }else{  ### 否则， 状态不变， 代数加1

           chain[i+1,] = chain[i,]

       }

   }

   return(chain)

}

res.chain1 <- MCMC(intialvalue = c(4,20,1), generations = 100000) ###初始值设定为 beta0 = 4, beta1 = 20, sigma = 1, 运行第一个马尔科夫链

res.chain2 <- MCMC(intialvalue = c(4,20,1), generations = 100000) ###初始值设定为 beta0 = 4, beta1 = 20, sigma = 1, 运行第二个马尔科夫链

####

par(mfrow = c(1, 3))

plot(res.chain1[,1], res.chain1[,2], type = "l")  ### 可以清楚得看出， 随机取值点逐渐偏离初始值， 最后稳定在一个值附近

plot(res.chain1[,1], res.chain1[,3], type = "l")

plot(res.chain1[,2], res.chain1[,3], type = "l")

图1 估计的参数的Burning In及马尔可夫链收敛后的变化

######################################################

############ 用coda程序包分析马尔可夫链 ##############

library(coda)

mcmc.res <- mcmc(res.chain1)

summary(mcmc.res)

plot(mcmc.res)

#### 马尔可夫链是否收敛， 进入平稳状态？ #############

#### 以下几个判别标准可以帮我们诊断， 读者可以参阅相应函数的说明

mcmcchains = mcmc.list(mcmc(res.chain1), mcmc(res.chain2))

plot(mcmcchains)

图2 trace plot以及概率密度分布

gelman.diag(mcmcchains)

gelman.plot(mcmcchains)

geweke.plot(mcmcchains)

geweke.diag(mcmcchains)

图3 用gelman.plot检验马尔可夫链是否收敛

#### 查看Densityplot图， 显示所求参数的概率分布情况

autocorr.plot(mcmc(res.chain1))

densityplot(mcmcchains)

par(mfrow = c(1, 3))

densplot(mcmcchains)

par(mfrow = c(1, 3))

traceplot(mcmc(res.chain1))

xyplot(mcmc(res.chain1))

###################################################

#### 当然， 从初始值开始， 一直到马尔可夫链到达平稳状态之前，这部分的结果是不能用的， 要去除的代数，常常约定俗成，如运行100万代， 500万代，舍去之前的10%。 为了减少抽象中自相关的影响， 还会每隔多少代取一个状态等。 这些内容， 绝大部分贝叶斯软件都有现成的函数。

beta0_actual = 20   ### beta0真实值为 20

beta1_actual = 5    ### beta1真实值为5

sigma_actual = 15   ### sigma真实值为15

burnIn = 20000

acceptance = 1 - mean(duplicated(res.chain1[-(1:burnIn),]))

par(mfrow = c(2,3))

hist(res.chain1[-(1:burnIn),1],nclass=50, , main="Posterior of beta0", xlab="True value = red line" )

abline(v = mean(res.chain1[-(1:burnIn),1]), col = "blue")

abline(v = linmod$coefficients[1], col="red" )

hist(res.chain1[-(1:burnIn),2],nclass=50, main="Posterior of beta1", xlab="True value = red line")

abline(v = mean(res.chain1[-(1:burnIn),2]), col = "blue")

abline(v = linmod$coefficients[2], col="red" )

hist(res.chain1[-(1:burnIn),3],nclass=50, main="Posterior of sigma", xlab="True value = red line")

abline(v = mean(res.chain1[-(1:burnIn),3]) , col = "blue")

abline(v = sigma_actual, col="red" ) ## 注意这里用到的sigma并非来自数据本身。

plot(res.chain1[-(1:burnIn),1], type = "l", xlab="Generations" , ylab = "Beta0", main = "Chain values of beta0", )

abline(h = linmod$coefficients[1], col="red" )

plot(res.chain1[-(1:burnIn),2], type = "l", xlab="Generations" , ylab = "Beta1", main = "Chain values of beta1", )

abline(h = linmod$coefficients[2], col="red" )

plot(res.chain1[-(1:burnIn),3], type = "l", xlab="Generations" , ylab = "Sigma", main = "Chain values of sd", )

abline(h = sigma_actual, col="red" ) ## 注意这里用到的sigma并非来自数据本身。

图 4 参数的Bayes估计值与真实值的比较

参考：

Florian Hartig / December 9, 2011

MCMC chain analysis and convergence diagnostics with coda in R

https://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/