这篇文章并非原创,只是对下面PPT的总结和理解,作为入门性文字说明。
http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc.pdf
MH算法:Metropolis Hasting Algorithm
Gibbs抽样:Gibbs Sampling
上述两个算法是两个典型的马尔科夫链蒙特卡洛方法(Markov Chain Mento Carlo Method),作为随机抽样方法,普遍用于多元随机变量的抽样。通俗的说,通过迭代方法从目标随机概率分布中抽取样本。
关于蒙特卡洛方法和马尔科夫链,以及大数法则和中心极限定理,这边的证明略去。PPT中有详细解释。
1. Gibbs 抽样
Gibbs抽样的本质就是从已知的联合概率分布p(θ1, θ2,....|Y)(后验分布)推导出满条件分布,然后从满条件分布中抽取样本——p(θj|θ-j, Y)
这里有个问题:为什么可以从联合概率分布中求解满条件分布?
因为根据Hammersley-Clifford Theorem,任何联合概率可以由他的条件概率计算得来。
有了以上基本概念后,Gibbs抽样可以总结,先求取满条件分布,再对满条件分布抽烟。
满条件分布求取:
a)先对后验分布进行整理,忽略相关的常数
b)分别对参数θ1,θ2,.......进行条件分布求解
c)正规化条件分布
Gibbs Sampler: 假设有三个参数θ1,θ2,θ3,后验分布为p(θ|Y)
a)给参数
θ 初始值,记作
θ(0),
θ取自于初始转移分布得来。
b)对任意一个参数进行基于满条件分布的参数估计。为方便起见,我们从θ1开始,从p(θ1|θ2(0),θ3(0),Y)中得到θ1(1)
c)同理,从p(θ2|θ1(1),θ3(0),Y)中得到θ2(1)
d)从p(θ3|θ1(1),θ2(1),Y)中得到θ3(1)
e)反复b,c,d步骤,直到求取到M个值之后停止。
f)可以对样本进行burn-in和thinning操作
例子:
一个核电厂有十个水泵,已知数据各水泵出故障次数,以及观察到出故障的时间
y <- c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22)
t <- c(94, 16, 63, 126, 5, 31, 1, 1, 2, 10)
建立关于水泵故障次数的泊松分布,λi 为单位时间每个水泵的故障次数。似然函数有:
∏i=1~10 Poisson(λi *ti)
假设Gamma(α, β) 是 λ 的先验分布。α=1.8. 所有的 λi都来自于这个分布。
假设Gamma(γ , δ) 是 β 的先验分布。γ=0.01 , δ=1
所以有:
p(λ, β|y, t) ∝ (∏i=1~10 Poisson(λi *ti)*Gamma(α, β))*Gamma(γ , δ)
经过整理有:
p(λ, β|y, t) ∝ ( ∏i=1~10 λi^(yi +α−1)*e^(−(ti +β)λi))*(β^(10α+γ−1)e^(−δβ))
通过上述联合分布,可以求得满条件分布:
p(λi |λ−i , β, y, t) ∝ λi^(yi +α−1)*e^(−(ti +β)λi)
p(β|λ, y, t)∝ β^(10α+γ−1)e^−β(δ+Sumi=1~10 λi)
gibbs<-function(n.sims,beta.start,alpha,gamma,delta,y,t,burnin=0,thin=1){
beta.draws<-c()
lambda.draws<-matrix(NA,nrow=n.sims,ncol=length(y))
beta.cur<-beta.start
lambda.update<-function(alpha,beta,y,t){
rgamma(length(y),y+alpha,t+beta)
}
beta.update<-function(alpha,gamma,delta,lambda,y){
rgamma(1,length(y)*alpha+gamma,delta+sum(lambda))
}
for(i in 1:n.sims){
lambda.cur<-lambda.update(alpha=alpha,beta=beta.cur,y=y,t=t)
beta.cur<- beta.update(alpha=alpha,gamma=gamma,delta=delta,lambda=lambda.cur,y=y)
if(i>burnin&&(i-burnin)%%thin==0){
lambda.draws[(i-burnin)%/%thin,]<-lambda.cur
beta.draws[(i-burnin)/thin]<-beta.cur
}
}
return(list(lambda.draws=lambda.draws,beta.draws=beta.draws))
}
其中lambda.update 方法和beta.update 方法是Gibbs Sampler过程。循环语句是对样本进行burn-in和thinning操作。
最后求取平均值,即所求参数λ, β的期望值。
posterior <- gibbs(n.sims = 10000, beta.start = 1, alpha = 1.8, gamma = 0.01, delta = 1, y = y, t = t)
colMeans(posterior$lambda.draws)
mean(posterior$beta.draws)
分享到:
相关推荐
吉布斯采样是统计学中用于马尔科夫蒙特卡洛的一种算法,用于在难以直接采样时从某一多变量概率分布中近似抽取样本序列。文档内有例子和代码以及运行结果。
这是一份英文文档,是比较详细地介绍LDA的一篇论文,写的比较基础,比LDA首创作者的论文好懂一些,最后的参数推理采用的是Gibbs抽样
gibbs算法示例,在统计学和统计物理学中,gibbs抽样是马尔可夫链蒙特卡尔理论(MCMC)中用来获取一系列近似等于指定多维概率分布(比如2个或者多个随即变量的联合概率分布)观察样本的算法。
利用条件概率和Gibbs抽样技术为分布估计算法构造通用概率模型利用条件概率和Gibbs抽样技术为分布估计算法构造通用概率模型
精要的证明推理,有关:LDA,主题模型,数学推导,GIbbs抽样 英文版的,但是还是很容易看懂 很短
我的课程作业……包括Metropolis,Metropolis Hastings, Laplace Approximation, Gibbs,Bayesian liner regression,Bayesian logistic regression的原理简单介绍和算法,水平有限一定会有错,发这就是为了保存...
基于EM算法和Gibbs抽样的污染模型的参数估计.pdf基于EM算法和Gibbs抽样的污染模型的参数估计.pdf
基于Gibbs抽样的贝叶斯随机波动模型分析,朱慧明,赵锐,随机波动性是经济金融时间序列中的一个普遍现象,它在金融风险管理研究中具有重要地位。通过分析随机波动模型的统计结构,推断了
针对现有动态面板数据分析中存在偶发参数和没有考虑模型参数的不确定性风险问题,提出了基于Gibbs抽样算法的贝叶斯随机系数动态面板数据模型。假设初始值服从平稳分布,自回归系数服从Logit正态分布的条件下,设计了...
gibbs_shiny Gibbs 采样器演示(Shiny 应用程序) 在查看
吉布斯采样matlab代码
信号与系统的连续时间信号的频域分析实验,有代码,有截图,有总结
实验4 周期信号的合成、分解与Gibbs现象
吉布斯采样matlab代码
研究简单回归模型中响应变量受到另一随机变量序列污染时,模型参数和污染系数的估计,运用 EM 算法给出了两类污染数据回归模型的参数的极大似然估计( MLE),并用 Gibbs抽样的方法给出了未知参数的 Monte-C arlo估计.
吉布斯抽样 用均值场法计算推论和推论。 此外,通过Jupyter Notebook的变量消除方法可以计算出准确的结果。 ================================================== ============================ [[Bayesian_roximate...
为解决基因数据变量多、抽样收敛速度慢的问题, 本文将 向前向后Gibbs算法应用到模型参数的马尔可夫链蒙特卡罗抽样估计中, 以加快收敛.在对多形性胶质母细胞瘤基 因数据分析中, 本文方法能有效识别出DNA拷贝数异常的...
异常值会使统计分析误差增大,为了识别这些异常值,本文给出了基于Gibbs抽样识别可加异常值的方法,并用我国人民币对美元汇率的月度数据进行实证研究。
运用基于Gibbs抽样的Markov Chain M0me Carlo贝叶斯方法,估计稳健ARMA模型参数,同步确定观测值中异常点的位置,辨别异常点类型。并利用我国人口自然增长数据进行仿真分析,研究结果表明:贝叶斯方法能够有效地识别...
采用Gibbs采样方法用于图像分割,涉及到K—means聚类的初始化以及形态学处理