I'm no expert on MCMC methods, so I'm not going to even attempt to explain how this works, but I did see how the different algorithms behave, and it is rather insightful. For the gamma distribution, while I always parameterize it using shape and scale, for the purposes of Gibbs sampling, I used shape/rate. The nice thging about Gibbs sampling is that as long as your can sample from a distribution which is proportional to the actual conditional distribution, you're fine, and do not have to calculate the constants of proportionality that bring the distribution back to an integral of 1.
In the case of the gamma distribution, and using a flat prior, it is the likelihood which dominates the conditional distribution. Using the shape/rate parameterization,
, the likelihood over the observations 
is:
When looking for a conditional distribution for alpha, this does not resolve into anything remotely nice (although see Fink (1997) "A Compendium of Conjugate Priors" (pdf) for more) so Gibbs sampling will have to give way to a Metropolis algorithm. I'm not sophisticated enough to use a non-symmetrical proposal distribution, so it is not a Metropolis-Hastings algorithm, maybe in a future post. However, when looking at beta, and discarding the x's and alphas, the core of the likelihood is rather gamma-like. Actually, it is proportional to a gamma distribution with the rate equal to the sum of the observed data points and the shape equal to one more than the product of the number of observations and the underlying alpha.
With that observation, and properly writing the algorithm (which took more time than I care to admit), a comparison between using Gibbs sampling or straight Metropolis on the scale parameter can be made.
Here is a "pure" Metropolis algorithm to calculate the posterior distribution, assuming the input data is gamma (this and all "pretty" R code created by Pretty R at inside-R.org :
MHGam <- function (x, n.iter=10000, n.burnin=0){ alpha <- matrix(.5, n.iter) theta <- matrix(10, n.iter) acca<- 0 accq <- 0 lnpost <- function(a, q, z){ NLL <- (a-1)*log(z)-(z/q)-a*log(q)-lgamma(a) return(sum(NLL)) } for (i in 2:n.iter) { #sample theta acc <- 1 theta[i] <- theta[i-1] + rnorm(1, mean=0, sd=100) if (theta[i]<=0) {acc=0; theta[i] <- theta[i-1]} LL <- lnpost(alpha[i-1], theta[i], x) if ((LL - lnpost(alpha[i-1], theta[i-1], x)) < log(runif(1))) { theta[i] <- theta[i-1] acc <- 0 } accq <- accq + acc #sample alpha acc <- 1 alpha[i] <- alpha[i-1] + runif(1, min=-.5, max=.5) if (alpha[i]<=0) {acc=0; alpha[i] <- alpha[i-1]} if ((lnpost(alpha[i], theta[i], x) - LL) < log(runif(1))) { alpha[i] <- alpha[i-1] acc <- 0} acca <- acca+acc if (i%%100==0) {print(c(i, alpha[i], theta[i], acca/i, accq/i))} } amean <- mean(alpha[n.burnin+1:n.iter]) qmean<-mean(theta[n.burnin+1:n.iter]) Fit <- structure(list(alpha=alpha, theta=theta, amean=amean, qmean=qmean)) return(Fit) }
The above code will just hand off the parameter sampling between alpha and theta (which is the scale and not rate) and will, hopefully, eventually settle down. I added two extra parameters for iterations and burnin, and eventually, I could add another, outer, for-next loop for multiple chains.
The next code takes advantage of Gibbs sampling for the scale parameter:
MHGamGib <- function (x, n.iter=10000, n.burnin=0){ alpha <- matrix(.5, n.iter) theta <- matrix(10, n.iter) acca<- 0 lnpost <- function(a, q, z){ NLL <- (a-1)*log(z)-(z/q)-a*log(q)-lgamma(a) return(sum(NLL)) } for (i in 2:n.iter) { #sample theta theta[i] <- 1/rgamma(1, rate=sum(Test), shape=alpha[i-1]*length(x)+1) #sample alpha acc <- 1 alpha[i] <- alpha[i-1] + runif(1, min=-.5, max=.5) if (alpha[i]<=0) {acc=0; alpha[i] <- alpha[i-1]} if ((lnpost(alpha[i], theta[i], x) - lnpost(alpha[i-1], theta[i], x)) < log(runif(1))) { alpha[i] <- alpha[i-1] acc <- 0} acca <- acca+acc if (i%%100==0) {print(c(i, alpha[i], theta[i], acca/i))} } amean <- mean(alpha[n.burnin+1:n.iter]) qmean<-mean(theta[n.burnin+1:n.iter]) Fit <- structure(list(alpha=alpha, theta=theta, amean=amean, qmean=qmean)) return(Fit) }
Note that I'm still parameterizing using the scale by sampling the rate and using its reciprocal.
To set up, run, and compare the models, I used the following code:
set.seed(22) Test <- rgamma(100, shape=2.5, scale=1e5) MHT<-system.time(MH<-MHGam(Test, n.iter=50000)) MHGT<-system.time(MHG<-MHGamGib(Test, n.iter=50000))
The Gibbs version runs about 40% faster, comparing MHT with MHGT (at least on my laptop), which stands to reason as Gibbs sampling is Metropolis sampling where every variate is accepted. What is more fascinating is how fast the parameter values converge:
Using the following code to overlay the trace plots for theta:
We see:
What is also surprising is how the better estimates for theta affect the random-walk Metropolis algorithm for alpha. Using the following code:
We get the following trace:
