R中的MCMC Gibbs采样器

Question

首先，我不太擅长编码（而且我不是编码员）——尤其是编码图表——这就是我需要帮助的原因。出于我个人的目的，我想玩 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,1], x[1:stepsToDisplay,2], pch=16, col="black", type="b", lty=2,cex = 1)
lines(x[1,1], x[1,2], pch=16, col="green", type="b", lty=2,cex = 1)
text(x[1:stepsToDisplay,1], x[1:stepsToDisplay,2], 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,1),x(t+1,1)],[x(t,2),x(t+1,2)],'k-');
    h2 = plot(x(t+1,1),x(t+1,2),'ko');
end

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

Answer 1

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

#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,1],x[t+1,1]),c(x[t,2],x[t,2]), lty=2,cex = 1)
    lines(c(x[t+1,1],x[t+1,1]),c(x[t,2],x[t+1,2]), lty=2,cex = 1)
    points(x[t+1,1],x[t+1,2], pch=16, col="black", lty=2,cex = 1)
    text(x[t+1,1], x[t+1,2], labels=t, cex= 0.7, pos=4)
}
points(x[1,1],x[1,2], pch=16, col="green", lty=2,cex = 1)
text(x[1,1], x[1,2], labels="start", cex= 0.7, 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)