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

enter image description here

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)) 



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.


  • 1
    $\begingroup$ The code does not implement the given formulas, see the second argument of 'rnorm'. $\endgroup$ – Yves May 24 at 12:40
  • $\begingroup$ Thanks for catching that! $\endgroup$ – 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!!
  • $\begingroup$ 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) $\endgroup$ – T.J May 24 at 12:32
  • $\begingroup$ This is not an issue with the MCMC but with the Bayesian analysis. On principle the prior hyperparameters should not depend on the data. $\endgroup$ – Xi'an May 24 at 12:37
  • 4
    $\begingroup$ 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. $\endgroup$ – Yves May 24 at 13:22
  • 6
    $\begingroup$ 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. $\endgroup$ – alan ocallaghan May 24 at 13:33
  • 1
    $\begingroup$ @alanocallaghan: thank you for mentioning the R inferno book, it is a great collection of horror stories. $\endgroup$ – Xi'an Jun 10 at 7:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.