文章目录
练习:用程序实现正态分布均值、方差的后验分布抽样。
题目背景
记
X
i
,
i
=
1...
,
n
{X_i,i=1...,n}
Xi,i=1...,n来自一维正态
N
(
μ
,
σ
2
)
N(\mu,\sigma^2)
N(μ,σ2)的独立同分布样本,
μ
\mu
μ具有先验分布
N
(
a
,
b
)
N(a,b)
N(a,b),
σ
2
\sigma^2
σ2具有先验分布
i
n
v
G
a
m
m
a
(
c
,
d
)
invGamma(c,d)
invGamma(c,d).
似然函数:
π
(
X
∣
μ
,
σ
2
)
∝
(
σ
2
)
−
n
/
2
exp
{
−
∑
i
(
X
i
−
μ
)
2
2
σ
2
}
\pi(X|\mu,\sigma^2) \propto (\sigma^2)^{-n/2} \exp\{ - \sum_i \frac{(X_i-\mu)^2}{2\sigma^2}\}
π(X∣μ,σ2)∝(σ2)−n/2exp{−i∑2σ2(Xi−μ)2}
补充
i
n
v
G
a
m
m
a
(
α
,
λ
)
invGamma(\alpha,\lambda)
invGamma(α,λ)的密度函数形式:
p
(
x
;
α
,
λ
)
=
λ
α
Γ
(
α
)
x
−
(
α
+
1
)
exp
{
−
λ
x
}
,
x
≥
0
p(x;\alpha,\lambda)=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}x^{-(\alpha+1)}\exp \{ -\frac{\lambda}x \}, \quad x\ge 0
p(x;α,λ)=Γ(α)λαx−(α+1)exp{−xλ},x≥0
Gibbs抽样(详细公式推导)
μ
\mu
μ的条件(后验)分布:*
p
(
μ
∣
X
,
σ
2
)
∝
π
(
X
∣
μ
,
σ
2
)
π
(
μ
)
∝
exp
{
−
∑
i
(
X
i
−
μ
)
2
2
σ
2
}
exp
{
−
(
μ
−
a
)
2
2
b
}
∝
exp
{
−
(
μ
−
μ
n
)
2
2
σ
n
2
}
∼
N
(
μ
n
,
σ
n
2
)
p(\mu|X,\sigma^2) \propto \pi(X|\mu,\sigma^2)\pi(\mu) \propto \exp\{ - \sum_i \frac{(X_i-\mu)^2}{2\sigma^2} \} \exp\{ -\frac{(\mu-a)^2}{2b}\} \\ \propto \exp\{ -\frac{(\mu-\mu_n)^2}{2\sigma_n^2} \} \sim N(\mu_n,\sigma_n^2)
p(μ∣X,σ2)∝π(X∣μ,σ2)π(μ)∝exp{−i∑2σ2(Xi−μ)2}exp{−2b(μ−a)2}∝exp{−2σn2(μ−μn)2}∼N(μn,σn2)
μ
n
=
(
∑
i
X
i
σ
2
+
a
b
)
/
(
n
σ
2
+
1
b
)
\mu_n=(\frac{\sum_i X_i}{\sigma^2}+\frac{a}{b})/(\frac{n}{\sigma^2} + \frac{1}{b})
μn=(σ2∑iXi+ba)/(σ2n+b1),
σ
n
2
=
(
n
σ
2
+
1
b
)
−
1
\sigma_n^2=(\frac{n}{\sigma^2}+\frac{1}{b})^{-1}
σn2=(σ2n+b1)−1.
σ
2
\sigma^2
σ2的条件(后验)分布:
p
(
σ
2
∣
X
,
μ
)
∝
π
(
X
∣
μ
,
σ
2
)
π
(
σ
2
)
∝
σ
−
n
exp
{
−
∑
i
(
X
i
−
μ
)
2
2
σ
2
}
σ
−
2
(
c
+
1
)
exp
{
−
d
σ
2
}
∝
σ
−
2
(
n
/
2
+
c
+
1
)
exp
{
−
1
σ
2
[
∑
i
(
X
i
−
μ
)
2
2
+
d
]
}
∼
i
n
v
G
a
m
m
a
(
c
n
,
d
n
)
p(\sigma^2|X,\mu) \propto \pi(X|\mu,\sigma^2)\pi(\sigma^2) \\ \propto \sigma^{-n}\exp\{ - \sum_i \frac{(X_i-\mu)^2}{2\sigma^2} \} \sigma^{-2(c+1)}\exp \{ -\frac{d}{\sigma^2}\} \\ \propto \sigma^{-2(n/2+c+1)} \exp\{-\frac{1}{\sigma^2}[\sum_i \frac{(X_i-\mu)^2}{2}+d ] \} \sim invGamma(c_n,d_n)
p(σ2∣X,μ)∝π(X∣μ,σ2)π(σ2)∝σ−nexp{−i∑2σ2(Xi−μ)2}σ−2(c+1)exp{−σ2d}∝σ−2(n/2+c+1)exp{−σ21[i∑2(Xi−μ)2+d]}∼invGamma(cn,dn)
c
n
=
n
/
2
+
c
c_n=n/2+c
cn=n/2+c,
d
n
=
∑
i
(
X
i
−
μ
)
2
2
+
d
d_n=\sum_i \frac{(X_i-\mu)^2}{2}+d
dn=∑i2(Xi−μ)2+d
Gibbs采样R代码实现
μ \mu μ的更新函数
#给定观测值和其他参数,更新参数mu
updt_mu_r = function(X,sigma2,prior){ #输入X,上一步迭代中的sigma,超参数a,b的先验值
n = length(X)
sumX = sum(X)
tmp = 1 / (n/sigma2 + 1/prior[2]) #mu的后验分布方差
mu_n = ( sumX/sigma2 + prior[1]/prior[2] ) * tmp #mu的后验分布均值
sigma_n = sqrt(tmp) #mu的后验分布标准差
return( rnorm(1,mu_n,sigma_n) ) #返回一个从mu的后验分布中的抽样值即为mu的更新值
}
σ 2 \sigma^2 σ2的更新函数
updt_sigma2_r = function(X,mu,prior){ #输入X,上一步迭代中的mu,超参数c,d的先验值
n = length(X)
shape = 0.5*n + prior[1] #sigma^2的后验分布形状参数
scale = sum((X-mu)^2)
scale = scale/2 + prior[2] #sigma^2的后验分布尺度参数
return( 1/rgamma(1,shape,scale) ) #返回一个从sigma^2的后验分布中的抽样值即为sigma^2的更新值
}
主函数
mcmc_Gaussian_r = function(X, n_sample, burn_in, thin, init, prior_mu, prior_sigma){
#输入X,采样的量,燃烧期,保存的间隔数,参数初始值,mu的先验分布中的超参数值,sigma2的先验分布中的超参数值。
mu_sample = rep(0,n_sample)
sigma2_sample = rep(0,n_sample)
#赋予初值
mu = init[1]
sigma2 = init[2]
#燃烧期
for(i in 1:burn_in){ #反复更新采样后验均值和方差
mu = updt_mu_r(X, sigma2, prior_mu)
sigma2 = updt_sigma2_r(X,mu,prior_sigma)
}
if(thin<=0) thin=1 #thin不能小于1,表示每隔thin-1次迭代记录一次数据。
for(i in 1:n_sample){
for(j in 1:thin){
mu = updt_mu_r(X, sigma2, prior_mu)
sigma2 = updt_sigma2_r(X,mu,prior_sigma)
}
mu_sample[i] = mu
sigma2_sample[i] = sigma2
}
out = list(mu = mu_sample, sigma2 = sigma2_sample)
return(out) #返回列表,out[[1]]是mu的模拟值,out[[2]]是sigma2的模拟值
}
案例求解
#模拟从一个均值为1,方差为4的正态分布中抽样
a = rnorm(5000,1,2)
res2 = mcmc_Gaussian_r(a,n_sample = 1000, burn_in = 100, thin = 1,
init = c(10,10), prior_mu = c(0,100),prior_sigma = c(1,1))
#返回结果res为一个列表
部分结果展示:
> res2[[1]][1:10]
[1] 0.9787088 0.9541495 1.0188802 0.9628411 1.0027464 0.9689442 0.9722291 0.9486590
[9] 0.9947923 0.9592428
> res2[[2]][1:10]
[1] 3.840518 3.901282 3.968341 3.875591 3.974695 3.921036 3.888012 4.023784 3.890451
[10] 3.919828
很接近真实的均值和方差。
(补充)RCPP代码
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
void updt_mu(double &mu, NumericVector X, double sigma2, NumericVector prior){
int n = X.size();
double sumX = 0;
double mu_n = 0;
double sigma_n = 0;
double tmp = 0;
for(int i = 0; i < n; i++){
sumX += X[i];
}
tmp = 1 / (n/sigma2 + 1/prior[1]);
mu_n = ( sumX/sigma2 + prior[0]/prior[1] ) * tmp;
sigma_n = sqrt(tmp);
mu = rnorm(1,mu_n,sigma_n)[0];
}
// [[Rcpp::export]]
void updt_sigma2(double &sigma2, NumericVector X, double mu, NumericVector prior){
int n = X.size();
double shape = 0;
double scale = 0;
shape = 0.5*n + prior[0];
for(int i = 0; i < n; i++){
scale += (X[i] - mu) * (X[i] - mu);
}
scale = scale/2 + prior[1];
sigma2 = 1/rgamma(1,shape,1/scale)[0] ;
}
// [[Rcpp::export]]
double loglik(NumericVector X, double mu, double sigma2){
int n = X.size();
double out = 0;
for(int i=0; i < n; i++){
out -= (X[i] - mu)*(X[i] - mu);
}
out = out/2/sigma2;
out -=-0.5*n * log(2*3.1415926*sigma2);
return out;
}
// [[Rcpp::export]]
List mcmc_Gaussian(NumericVector X, int n_sample, int burn_in, int thin,
NumericVector init, NumericVector prior_mu, NumericVector prior_sigma){
NumericVector mu_sample(n_sample);
NumericVector sigma2_sample(n_sample);
NumericVector loglk(n_sample);
double mu = init[0];
double sigma2 = init[1];
for(int i = 0; i < burn_in; i++){
updt_mu(mu, X, sigma2, prior_mu);
updt_sigma2(sigma2, X, mu, prior_sigma);
}
if(thin<=0) thin=1;
for(int i = 0; i < n_sample; i++){
/*
if(i%100==0){
std::cout << "iteration: " << i << " mu: " << mu << std::endl;
}
*/
for(int j = 0; j < thin; j++){
updt_mu(mu, X, sigma2, prior_mu);
updt_sigma2(sigma2, X, mu, prior_sigma);
}
mu_sample[i] = mu;
sigma2_sample[i] = sigma2;
loglk[i] = loglik(X, mu, sigma2);
}
return List::create(
_["mu"] = mu_sample,
_["sigma2"] = sigma2_sample,
_["loglik"] = loglk
);
}
MH采样代码
MH采样原理不细说。
相同的案例。
直接上代码,来对比一下两种采样方法。
不好意思,这里只上RCPP代码,RCPP可以找资料学,不想学也不影响从代码中观看算法的逻辑!因为我也是RCPP入门级别选手~
小tips:老师说凡是遇到求密度,尤其是联合密度值,在代码中一定一定要取对数运算(不是对数变换,就是取对数后按照对数的计算规则进行后续一切运算),原因是密度在0~1之间,连续相乘后会非常接近0,很可能导致数值溢出啥的。
// [[Rcpp::export]]
//求X样本的联合概率密度的对数的函数,输入X,mu,sigma2,输出
double loglik(NumericVector X, double mu, double sigma2){
int n = X.size(); //初始化变量n为数据X的维度大小
double out = 0; //初始化变量out=0
for(int i=0; i < n; i++){ //i从0到n-1的循环
out -= (X[i] - mu)*(X[i] - mu); //这里在求X的联合密度的exp()内的一部分
}
out = out/2/sigma2; //这里输出X的联合密度的exp()内的完整部分
out -= 0.5*n * log(2*3.1415926*sigma2); //这里右边的式子计算的是X密度的exp()外面的数取对数,整个式子输出完整的X对数密度值
return out;
}
// [[Rcpp::export]]
//MH更新参数的函数,注意到MH一次性更新所有参数,而Gibbs是逐个参数进行更新。
//该函数输入上一步迭代中的mu,sigma2,数据X,学习率lr,实际更新次数count。
void updt_para_MH(double &mu, double &sigma2, NumericVector X, double lr, double &count){
//同时更新mu和sigma2,lr为学习率。更新步长采样自正态随机数(乘上学习率)
NumericVector z0 = rnorm(2) * lr;
double mu_new = mu + z0[0];
double sigma2_new = sigma2 + z0[1];
double acc = 0;
if(sigma2_new > 0){ //必须保证采样空间是正确的:sigma2为正数
NumericVector test = runif(1); //用来进行接受拒绝采样的随机数
double logL_new = loglik(X,mu_new,sigma2_new); //计算新的对数似然
double logL_old = loglik(X,mu,sigma2); //上一步参数对应的对数似然
if( log(test[0]) < (logL_new - logL_old) ){ //如果均匀分布随机数小于新旧密度之比,则接受采样值
acc = 1; //确定更新
count++; //count记录真实的更新次数
}
}
if(acc==1){ //确定更新
mu = mu_new;
sigma2 = sigma2_new;
}
}
// [[Rcpp::export]]
//主函数:输入X,采样的次数,燃烧期,采样间隔记录的步数,参数初始值,学习率;输出参数MH采样值
List mcmc_Gaussian_MH(NumericVector X, int n_sample, int burn_in, int thin,
NumericVector init, double lr){
//定义以及初始化变量
NumericVector mu_sample(n_sample);
NumericVector sigma2_sample(n_sample);
NumericVector loglk(n_sample);
double mu = init[0];
double sigma2 = init[1];
double count = 0; //这里初始化count=0
//利用上一个参数更新函数进行参数的更新
for(int i = 0; i < burn_in; i++){
//updt_para_MH(mu,sigma2,X,lr,); //偷懒省去燃烧期,更新过程与下面类似~
}
if(thin<=0) thin=1; //thin必须为正整数
for(int i = 0; i < n_sample; i++){ //循环n_sample次
for(int j = 0; j < thin; j++){
updt_para_MH(mu,sigma2,X,lr,count); //更新参数!
if(i%100==0){ #i每100进行下述操作
if(count > 30){ #如果实际更新的次数大于30次进行下述操作
lr = lr * 1.1; #学习率增大0.1倍
}else if(count < 20){ #否则学习率缩小0.1倍
lr = lr * 0.9;
}
//上面的操作实际上是随着采样过程的进行,每100次采样中,当参数实际更新次数<20次时,学习率就衰减为原来的0.9倍;实际更新次数在20~30时,不改变学习率;实际更新次数>30时,学习率增大为原来的1.1倍.原因可能是怕刚开始走太快偏离方向于是慢慢走,后来趋于稳定变化很慢于是加大步伐?不懂这些实操技术,有些玄妙~
std::cout << count/100 << " lr: " << lr << std::endl; //不用管,用来在控制台输出学习率的变化情况~
count = 0; //count重置为0(该变量的时候加了"&"有关,像是内置参数)
}
}
mu_sample[i] = mu; //记录
sigma2_sample[i] = sigma2;
loglk[i] = loglik(X, mu, sigma2);
}
//定义输出列表
return List::create(
_["mu"] = mu_sample,
_["sigma2"] = sigma2_sample,
_["loglik"] = loglk
);
}
总结:
当参数的后验分布有显式的表达时,MCMC抽样不是必要的,在进行推断的时候有可能根据分布直接进行计算(如Gibbs);当后验分布没有显式表达时(归一化参数在很多时候是算不出来的),可以通过MCMC算法进行抽样。