R中的MCMC Gibbs采样器

Question

首先，我不擅长编码(而且我不是编码员(，尤其是对图表的编码，这就是我需要帮助的原因。出于我个人的目的，我想玩MCMC吉布斯采样，我发现了以下MATLAB代码：

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

然而，比起MATLAB，我更喜欢R。我认为我自己转换得很好代码的最大部分：

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)