是Gibbs抽样吧
这里有一个例子，你看看有没有用
一个核电厂有十个水泵，已知数据各水泵出故障次数，以及观察到出故障的时间
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.startlambda.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.curbeta.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)
扫码加好友，拉您进群
全部版块
我的主页
收藏
