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

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

介绍贝叶斯推断， 则绕不过贝叶斯公式， 这是一个表示参数集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 抽样）估计线性回归模型的参数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

### 第一步 首先生成要拟合的数据
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估计值与真实值的比较

参考：