It is relatively easy to set up a Gibbs sampling algorithm for the normal sampling problem when independent priors (of the conjugate type) are assigned to the mean and precision. Here we outline how to do this on R.

We start with an expression for the joint posterior of the mean

(Here S is the sum of squares of the observations about the mean.)

1. To start, we recognize the two conditional distributions.

- The posterior of
given is given by the usual updating formula for a normal mean and a normal prior. (Essentially this formula says that the posterior precision is the sum of the prior and data precisions and the posterior mean is a weighted average of the prior mean and the sample mean where the weights are proportional to the corresponding precisions. - The posterior of
given has a gamma form where the shape is given by and the scale is easy to pick up.

2. Now we’re ready to use R. I’ve written a short function that implements a single Gibbs sampling cycle. To understand the code, here are the variables:

– ybar is the sample mean, S is the sum of squares about the mean, and n is the sample size

– the prior parameters are (mu0, tau) for the prior and (a, b) for the precision

– theta is the current value of (

The function performs the simulations from the distributions [

`one.cycle=function(theta){`

mu = theta[1]; P = theta[2]

P1 = 1/tau0^2 + n*P

mu1 = (mu0/tau0^2 + ybar*n*P) / P1

tau1 = sqrt(1/P1)

mu = rnorm(1, mu1, tau1)

a1 = a + n/2

b1 = b + S/2 + n/2*(mu – ybar)^2

P = rgamma(1, a1, b1)

c(mu, P)

}

All there is left in the programming is some set up code (bring in the data and define the prior parameters), give a starting value, and collect the vectors of simulated draws in a matrix.