sampling a multimensional posterior distribution using MCMC Metropolis-Hastings algo in R

Question

I am quite new in sampling posterior distributions(so therefore Bayesian approach) using a MCMC technique based on Metropolis-Hastings algorithm. I am using the mcmc library in R for this. My distribution is multidimensionnal. In order to check if this metro algorithm works for multivaiate distribution I did it successfully on a multidimensional student-t distribution (package mvtnorm, function dmvt). Now I want to apply the same thing to my multivariate distribution (2 vars x and y) but it doesn't work; I get an error : Error in X[, 1] : incorrect number of dimensions

Here is my code:

library(mcmc)
library(mvtnorm)
my.seed <- 123

logprior<-function(X,...)
{
      ifelse( (-50.0 <= X[,1] & X[,1]<=50.0) & (-50.0 <= X[,2] & X[,2]<=50.0), return(0), return(-Inf))
}

logpost<-function(X,...)
{
      log.like <- log( exp(-((X[,1]^2 + X[,2]^2 - 4)/10 )^2) * sin(4*atan(X[,2]/X[,1])) )
      log.prior<-logprior(X)
      log.post<-log.like + log.prior # if flat prior, the posterior distribution is the likelihood one
      return (log.post)
}

x <- seq(-5,5,0.15)
y <- seq(-5,5,0.15)
X<-cbind(x,y)

#out <- metrop(function(X) dmvt(X, df=3, log=TRUE), 0, blen=100, nbatch=100) ; this works
out <- metrop(function(X) logpost(X), c(0,0), blen=100, nbatch=100)
out <- metrop(out)
out$accept

So I tried to respect the same kind of format than for the MWE, but it doesn't work still as I got the error mentioned before. Another thing, is that applying logpost to X works perfectly.

Thanks in advance for your help, best

score 0 · Answer 1 · answered Nov 16 '18 at 17:57

The metrop function passes individual samples, and therefore a simple vector to logpost, not a matrix (which is what X is). Hence, the solution is to change X[,1] and X[,2] to X[1] and X[2], respectively.

I ran it like this, and it leads to other issues (X[2]/X[1] is NaN for the initialization), but that has more to do with your specific likelihood model and is out of the scope of your question.

sampling a multimensional posterior distribution using MCMC Metropolis-Hastings algo in R

1 Answers1