正态后验分布的MCMC方法

$$\begin{gathered} 代码是对原文例子的一个验证，用的R语言，使用的Metropolis-Hastings方法，写的有点潦草见谅见谅， \\ 原文链接 \end{gathered}$$: 概率与统计——正态分布的共轭分布.

一. 正态分布后验分布的MCMC方法

1. 期望$\mu$未知，方差$\phi$已知的情形

#利用MCMC对正态分布(期望未知，方差已知)的后验分布进行模拟
n   =100 		   #生成100个样本
var = 10		   #总体正态分布的已知方差
X = rnorm(n,23,var)#以23为期望的正态分布生成n个样本，即假定真实期望为23

theta_0=23		 #总体期望所服从的先验分布的均值
phi_0 = 2.25     #总体期望所服从的先验分布的方差
niter = 10^5     #链中状态转移步数
mu = rep(0,niter)#存储每一步后的总体期望值
mu[1]=X[1]		 #总体期望值初始状态赋值

for(i in 2 : niter){
  mu.p = mu[i-1]+rnorm(1,0,1)#利用一个标准正态分布来驱动总体期望值状态转移
  m = rep(mu.p,n)
  s_1 = exp(-(mu.p-theta_0)^2/(2*phi_0))*exp(-(t(X - m)%*%(X - m))[1]/(2*var))
  m = rep(mu[i-1],n)
  s_0 = exp(-(mu[i-1]-theta_0)^2/(2*phi_0))*exp(-(t(X - m)%*%(X - m))[1]/(2*var))
  r = s_1/s_0
  flip = rbinom(1,1,min(r,1))
  mu[i]=if(flip==1)mu.p else mu[i-1]
}
par(mfrow=c(1,2))
h_mu=mu[-(1:niter/2)]
hist(h_mu,breaks=100)
print(mean(h_mu))
x=seq(1,niter,1)
plot(x,h_mu,type="l",xlab="x",ylab="mu",main="MCMC",xlim=c(0,niter),ylim=c(min(h_mu),max(h_mu)),col="blue")

说明
图示里的MCMC结果虽然比较符合均值为23的正态分布，但并不是每次都能得到比较理想的直方图，有时会出现期望值偏得比较多的情况，这些都与生成的样本向量有关

2. 期望$\mu$已知，方差$\phi$未知的情形

#利用MCMC对正态分布(期望已知，方差未知)的后验分布进行模拟phi=sigma^2
t_0 = 15		  #预估的先验分布上的自由度t_0
S_0 = 585         #预估的先验分布上的参数S_0

S = 724			  #样本数据与已知总体期望的差的平方和
n = 20			  #样本数据的个数

niter = 10^5	  #链中状态转移步数
phi = rep(0,niter)#存储每一步后的总体方差值
phi[1] = 45 	  #总体方差初始状态赋值


for(i in 2 : niter){
  phi.p = phi[i-1]+rnorm(1,0,3)#利用一个均值0，方差3正态分布来驱动总体期望值状态转移
  r = phi.p^(-n/2)*exp(-S/(2*phi.p)) * dchisq(S_0/phi.p,t_0)*(S_0/phi.p^2)/(phi[i-1]^(-n/2)*exp(-S/(2*phi[i-1])) * dchisq(S_0/phi[i-1],t_0)*(S_0/phi[i-1]^2))
  flip = rbinom(1,1,min(r,1))
  phi[i]=if(flip==1)phi.p else phi[i-1]
}

par(mfrow=c(1,2))
h_phi=phi[-(1:niter/2)]
hist(h_phi,breaks=100)
print(mean(h_phi))
x=seq(1,niter,1)
plot(x,phi,type="l",xlab="x",ylab="phi",main="MCMC",xlim=c(0,niter),ylim=c(min(phi),max(phi)),col="blue")

说明
方差的MCMC结果是比较稳定的，直方图近似于一个逆卡方分布，方差状态的均值整体稳定在39~40之间，与例题答案39.67比较接近