代码是对原文内容的一个检验，用的R语言，写的有点潦草见谅见谅，原文链接: 概率与统计——正态分布的共轭分布.

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

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

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

d = 1
niter = 10^5
mu = rep(0,niter)
mu[1]=X[1]

for(i in 2 : niter){
  mu.p = mu[i-1]+rnorm(1,0,d)
  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))#不要对s_1和s_0的公式进行化简，否则会造成e的幂非常大而出现计算错误
  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,mu,type="l",xlab="x",ylab="mu",main="MCMC",xlim=c(0,niter),ylim=c(min(mu),max(mu)),col="blue")

在这里插入图片描述

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

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

#利用MCMC对正态分布(期望已知，方差未知)的后验分布进行模拟phi=sigma^2
phi_0=45
t_0 = 15
S = 724
S_0 = 585
n=20
d = 3
niter = 10^5
phi = rep(0,niter)
phi[1] = phi_0

for(i in 2 : niter){
  phi.p = phi[i-1]+rnorm(1,0,d)
  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")