张金龙的博客分享 http://blog.sciencenet.cn/u/zjlcas 物种适应性、分布与进化

博文

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

已有 8724 次阅读 2014-11-13 03:54 |个人分类:统计分析|系统分类:科研笔记

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


张金龙

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)依据当前参数组合和起始点概率密度比值是否大于01之间的一个均匀分布的随机数来判断是否保留当前点。

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/




https://blog.sciencenet.cn/blog-255662-843026.html

上一篇:线性模型参数的极大似然估计
下一篇:herblabel R程序包打印植物标本标签 (2017年1月9日更新)

5 余国志 斯幸峰 王谢 许格希 廖晓琳

该博文允许注册用户评论 请点击登录 评论 (4 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-7-4 17:21

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部