# Gibbs Sampler for Normal and Inverse Gamma Distribution in R

I'm trying to implement a Gibbs sampler for the following conditional distributions using R:

This is the code I have in R so far:


gibbs = function(N, v0, alpha, beta){

# Data
X = c(x1,...,xn);

# Parameters
n = length(X);
xbar = mean(X);
sigma2 = var(X);

# Starting values
mu0 = v0[1];
sigma20 = v0[2];
out = NULL;

for (i in 1:N){
mu0 = rnorm(1, mean = ((n*xbar/sigma2 + mu0/sigma20) / (n/sigma2) + (1/sigma20)), sd = n/sigma2 + (1/sigma20)^-1);
sigma20 = rinvgamma(1, n + alpha, 0.5*(sum(X - mu0))^2 + beta);
out = rbind(out, c(mu0, sigma20))
}

return(out)

}



When I run the function using gibbs(1000, v0 = c(1,1), 0.001, 0.01), the code runs for 5 or so iterations before starting to output NaN for everything.

I'm new to R, and to MCMC methods, so any help would be greatly appreciated.

Thanks!

• The code does not implement the given formulas, see the second argument of 'rnorm'. – Yves May 24 at 12:40
• Thanks for catching that! – T.J May 24 at 13:18

Several conceptual and R-coding errors:

1. $$\mu_0$$ and $$\sigma_0$$ are hyper-parameters for the prior, hence should not be modified along iterations

2. the conditionals in the Gibbs sampler are about the sampling model parameters $$\mu$$ and $$\sigma$$ which means updating mu=rnorm(... and sigma2=rinvgamma(... in the R code

3. there was a factor 2 missing in n/2 in the shape of the rinvgamma call

4. The term sum(X - mu0))^2 is doubly wrong as it should be sum((X - mu)^2)

5. the power -1 is not correctly set in the standard deviation in the rnorm call, it should be

(n/sigma2 + 1/sigma20)^(-1/2)

6. the rinvgamma function in MCMCpack is parameterised in terms of scale and shape, not rate and shape, hence the second parameter is the inverse of what it should be:

sigma2 = rinvgamma(1, n/2 + alpha, 1/(0.5*sum((X - mu)^2) + beta))

which is the principal reason for the diverging chain in the original version.

The entire R code could thus be

gibbs = function(N, v0, alpha, beta){
# parameters
n = length(X)
nxbar = sum(X)
sigma2 = var(X)
# hyperparameters
mu0 = v0[1]
sigma20 = v0[2]
out = NULL
for (i in 1:N){
mu = rnorm(1,
mean = (nxbar/sigma2 + mu0/sigma20) / (n/sigma2) + (1/sigma20),
sd = 1/sqrt(n/sigma2 + 1/sigma20))
sigma2 = rinvgamma(1, n/2 + alpha, 1/(sum((X - mu)^2)/2 + beta))
out = rbind(out, c(mu, sigma2)) # yuk!!
}
return(out)}

• Hi, thank you so much for picking those out! Quick question: regarding the hyper-parameters, should I specify these based on the data? e.g. μ = mean(X) – T.J May 24 at 12:32
• This is not an issue with the MCMC but with the Bayesian analysis. On principle the prior hyperparameters should not depend on the data. – Xi'an May 24 at 12:37
• Note that rnorm comes from the stats rpackage which is a base package and can be thought of as part of R. The inverse gamma distribution does not exist yet in stats and to has to be taken from another R package. There are some differences across packages in the parameterizations for this distribution so it could help to mention the package used. – Yves May 24 at 13:22
• It's worth noting that out = rbind(out, c(mu, sigma2)) is an example of "growing objects", a very inefficient way to store MCMC iterations, see the R inferno for details. – alan ocallaghan May 24 at 13:33
• @alanocallaghan: thank you for mentioning the R inferno book, it is a great collection of horror stories. – Xi'an Jun 10 at 7:03