Multavariate Gibbs Samples

1. Algoritmo no R:

library(mnormt)
library(ggplot2)
library(gridExtra)

# (C.) Anne Sabourin, 2009
#a ###########################################
T=500 ;p=5 ;r=0.25
X=cur=rnorm(p)
for (t in 1 :T){
  for (j in 1 :p){
    m=sum(cur[-j])/(p-1)
    cur[j]=rnorm(1,(p-1)*r*m/(1+(p-2)*r),
                 sqrt((1+(p-2)*r-(p-1)*r^2)/(1+(p-2)*r)))
  }
  X=cbind(X,cur)
}

# no ggplot

X = X[,-1]
colnames(X) = paste(1:T);
rownames(X) = paste("X", paste(1:p), sep="")

qq = seq(-4,4, len=500)
teorico = dnorm(qq)

X = rbind(X, teorico, qq)

val=as.data.frame(t(X))
dados <- ggplot(data=val)

cores = c("blue", "red", "green", "black")
xx <- list()
for(i in 1:p){
  xx[[i]] <- dados + geom_histogram(aes(x=val[,i], y=..density..), fill=cores[i], alpha=0.5, colour="black") +
    geom_line(aes(x=qq,y=teorico))
}

do.call("grid.arrange", c(xx))

2. Funçao em Rcpp

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix funcaoA(int T, int p, double r) {
  
  RNGScope scope;
  NumericMatrix X(T+1,p);
  X(0,_) = rnorm(p);
  NumericVector z;
  double m;
  NumericVector cur;

  for(int i = 0; i < T; i++){
    for(int j = 0; j < p; j++){
      z = X(i,_);
      m = std:: accumulate(z.begin(),z.end(), 0.0) - X(i,j);
      m = m/(p-1);
      cur = rnorm(1,(p-1)*r*m/(1+(p-2)*r),sqrt((1+(p-2)*r-(p-1)*pow(r,2))/(1+(p-2)*r)));
      X(i,j) = cur[0];
    }
  }

//  NumericMatrix X2 = X(Range(0,T-1),Range(0,p-1));
  
  return X;
}

T=500 ;p=5 ;r=0.25
X = funcaoA(T,p,r)
X = t(X)

X = X[,-1]
colnames(X) = paste(1:T);
rownames(X) = paste("X", paste(1:p), sep="")

qq = seq(-4,4, len=500)
teorico = dnorm(qq)

X = rbind(X, teorico, qq)

val=as.data.frame(t(X))
dados <- ggplot(data=val)

cores = c("blue", "red", "green", "black")
xx <- list()
for(i in 1:p){
  xx[[i]] <- dados + geom_histogram(aes(x=val[,i], y=..density..), fill=cores[i], alpha=0.5, colour="black") +
    geom_line(aes(x=qq,y=teorico))
}

do.call("grid.arrange", c(xx))