微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

R 中的 MCMC Gibbs 采样器

如何解决R 中的 MCMC Gibbs 采样器

首先,我非常不擅长编码(而且我不是编码员)——尤其是编码图表——这就是我需要一些帮助的原因。出于我个人的目的,我想使用 MCMC Gibbs 采样,我发现了以下 MATLAB 代码

https://theclevermachine.wordpress.com/2012/11/05/mcmc-the-gibbs-sampler/

但是我喜欢 R 多过 MATLAB。我认为我自己转换得很好,代码的最大部分:

library("phonTools")
Nsamples<-5000
mu<-c(0,0) #moyenne cible
rho<-c(0.8,0.8) #rho_21 rho_12

#initialisation de l'échantillon de Gibbs
propSigma<-1
minn<-c(-3,-3)
maxx<-c(+3,+3)

#on initialise les échantillons
x<-phonTools::zeros(Nsamples,2)
x[1,1]<-runif(1,min = minn[1],max = maxx[1])
x[1,2]<-runif(1,min = minn[2],max = maxx[2])

dims<- 1:2 #index dans chaque dimesion

#on exécute l'échantillonnage de Gibbs
t<-1
while (t < Nsamples) {
   t<-t + 1
   T<-c(t-1,t)
   for (iD in 1:2) { #on boucle sur les dimensions
     #on met à jour les échantillons
     nIx<-(dims!=iD)
     #moyenne conditionnelle
     muCond <- mu[iD] + rho[iD]*(x[T[iD],nIx]-mu[nIx]);
     #variance conditionnelle
     varCond <- sqrt(1-rho[iD]^2)
     x[t,iD] <-rnorm(1,mean=muCond,sd=varCond)
   }
}

#on affiche le graph
stepsTodisplay<-10
plot(x[,1],x[,2],main = "Gibbs Sampling",xlab = "x_1",ylab = "x_2",col="red",pch=19,cex = 0.5)


lines(x[1:stepsTodisplay,x[1:stepsTodisplay,pch=16,col="black",type="b",lty=2,cex = 1)
lines(x[1,x[1,col="green",cex = 1)
text(x[1:stepsTodisplay,labels=1:5,cex= 0.7,pos=3)

legend("bottomright",legend=c("Samples","1st 50 samples","x(t=0)"),col=c("red","black","green"),pch = c(16,16,16),cex=0.8)

我被困在从 MATLAB 转换以下视觉部分(对于习惯用 R 绘制图形的人来说,很可能很简单):

% CONDITIONAL STEPS/SAMPLES
hold on;
for t = 1:50
    plot([x(t,1),x(t+1,1)],[x(t,2),x(t,2)],'k-');
    plot([x(t+1,'k-');
    h2 = plot(x(t+1,'ko');
end

非常感谢任何帮助或改进建议

解决方法

我想我找到了解决方案。如果有人感兴趣,我会得到以下信息:

#pour les ellipses de confiance plus tard
library("car")

Nsamples<-500
mu<-c(0,0) #moyenne cible
rho<-c(0.8,0.8) #rho_21 rho_12

#initialisation de l'échantillon de Gibbs
propSigma<-1
minn<-c(-3,-3)
maxx<-c(+3,+3)

#on initialise les échantillons
#x<-phonTools::zeros(Nsamples,2)
x<-matrix( ncol=2,rep( 0,len=2*Nsamples))
x[1,1]<-runif(1,min = minn[1],max = maxx[1])
x[1,2]<-runif(1,min = minn[2],max = maxx[2])

dims<- 1:2 #index dans chaque dimesion

#on exécute l'échantillonnage de Gibbs
t<-1
while (t < Nsamples) {
   t<-t + 1
   T<-c(t-1,t)
   for (iD in 1:2) { #on boucle sur les dimensions
     #on met à jour les échantillons
     nIx<-(dims!=iD)
     #moyenne conditionnelle
     muCond <- mu[iD] + rho[iD]*(x[T[iD],nIx]-mu[nIx]);
     #variance conditionnelle
     varCond <- sqrt(1-rho[iD]^2)
     x[t,iD] <-rnorm(1,mean=muCond,sd=varCond)
   }
}

#on affiche le graph
stepsToDisplay<-5
car::dataEllipse(x[,1],x[,2],xlab = "x_1",ylab = "x_2",col="red",pch=19,cex = 0.5,levels=c(0.70,0.85,0.95,0.99),fill=TRUE,fill.alpha=0.15,lty=1,lwd=1,main="Bivariate Gibbs Sampler")

for(t in 1:stepsToDisplay){
    lines(c(x[t,x[t+1,1]),c(x[t,x[t,2]),lty=2,cex = 1)
    lines(c(x[t+1,cex = 1)
    points(x[t+1,pch=16,col="black",cex = 1)
    text(x[t+1,labels=t,cex= 0.7,pos=4)
}
points(x[1,x[1,col="green",cex = 1)
text(x[1,labels="start",pos=4)

legend("bottomright",legend=c("Samples",paste("1st",stepsToDisplay,"samples"),"x(t=0)"),col=c("red","black","green"),pch = c(16,16,16),cex=0.8)

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。