Gibbs Sampler for Poisson distribution ((Robert and Casella, 10.17)

Gibbs Sampler is used a lot in Bayesian inference. The following is an R function that output the mean and standard deviation of lambda1 to lambda10 and beta.

Gibbs.pump <- function(alpha = 1.8, r = 0.01, delta = 1, simulation = 1000, burn = 100) {
    ##  data
    y <- c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22)
    t <- c(94, 16, 63, 126, 5, 31, 1, 1, 2, 10)
    lambda <- numeric(10)

    ##  initial values
    beta <- 0.01

    ##  burn-in session
    for (i in 1:(burn + 1)) {       
        for (j in 1:10) {
            lambda[j] <- rgamma(1, y[j] + alpha, t[j] + beta)
        }

        beta <- rgamma(1, 10 * alpha + r, delta + sum(lambda))      
    }

    ##  markov chain session
    lambda2 <- matrix(NA, nrow = simulation, ncol = 10, byrow = TRUE)
    lambda2[1, ] <- lambda

    beta2 <- numeric(simulation)
    beta2[1] <- beta

    for (i in 2:simulation) {
        beta2[i] <- rgamma(1, 10 * alpha + r, delta + sum(lambda2[i-1, ]))

        for (j in 1:10) {
            lambda2[i, j] <- rgamma(1, y[j] + alpha, t[j] + beta2[i])
        }
    }

    ##  returned values
    lambda.mean <- apply(lambda2, 2, mean)
    lambda.sd <- apply(lambda2, 2, sd)
    lambda.dataframe <- data.frame(mean = lambda.mean, sd = lambda.sd)

    beta.mean <- mean(beta2)
    beta.sd <- sd(beta2)
    beta.dataframe <- data.frame(mean = beta.mean, sd = beta.sd)

    res <- list(lambda = lambda.dataframe, beta = beta.dataframe)
    return(res)
}

We can run the function and specify the burn-in number and simulation times.

> Gibbs.pump(burn = 1000, simulation = 1e4)
$lambda
         mean         sd
1  0.07051077 0.02698162
2  0.15207612 0.09187032
3  0.10421768 0.03991478
4  0.12289529 0.03097761
5  0.64916456 0.30120514
6  0.62272301 0.13671542
7  0.85414118 0.54445976
8  0.85574205 0.54394850
9  1.35889842 0.60607078
10 1.92709361 0.40330237

$beta
      mean        sd
1 2.406072 0.6915125

We can also calculate the other properties of the posteriors, like median and confidence interval, etc. The refined code is as follows.

Gibbs.pump <- function(alpha = .7, r = 0.01, delta = .01, simulation = 1e4, burn = 1e3) {
    ##  data
    y <- c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22)
    t <- c(94.32, 15.72, 62.88, 125.76, 5.24, 31.44, 1.048, 1.048, 2.096, 10.48)
    lambda <- numeric(10)

    ##  initial values
    beta <- 0.5

    ##  burn-in session
    for (i in 1:(burn + 1)) {       
        for (j in 1:10) {
            lambda[j] <- rgamma(1, y[j] + alpha, t[j] + beta)
        }

        beta <- rgamma(1, 10 * alpha + r, delta + sum(lambda))      
    }

    ##  markov chain session
    lambda2 <- matrix(NA, nrow = simulation, ncol = 10, byrow = TRUE)
    lambda2[1, ] <- lambda

    beta2 <- numeric(simulation)
    beta2[1] <- beta

    for (i in 2:simulation) {
        beta2[i] <- rgamma(1, 10 * alpha + r, delta + sum(lambda2[i-1, ]))

        for (j in 1:10) {
            lambda2[i, j] <- rgamma(1, y[j] + alpha, t[j] + beta2[i])
        }
    }

    ##  returned values
    lambda.mean <- apply(lambda2, 2, mean)
    lambda.sd <- apply(lambda2, 2, sd)
    lambda.median <- apply(lambda2, 2, median)
    lambda.lwr <- apply(lambda2, 2, function(x) {quantile(x, .025)})
    lambda.upr <- apply(lambda2, 2, function(x) {quantile(x, .975)})
    lambda.dataframe <- data.frame(mean = lambda.mean, 
                                                     sd = lambda.sd,
                                                     median = lambda.median,
                                                     X2.5 = lambda.lwr,
                                                     X97.5 = lambda.upr)

    beta.mean <- mean(beta2)
    beta.sd <- sd(beta2)
    beta.median <- median(beta2)
    beta.lwr <- quantile(beta2, .025)
    beta.upr <- quantile(beta2, .975)
    beta.dataframe <- data.frame(mean = beta.mean,
                                                  sd = beta.sd,
                                                  median = beta.median,
                                                  X2.5 = beta.lwr,
                                                  X97.5 = beta.upr, row.names = "1")

    res <- list(lambda = lambda.dataframe, beta = beta.dataframe)
    return(res)
}

Running the function we get

> Gibbs.pump()
$lambda
         mean         sd     median        X2.5     X97.5
1  0.05991165 0.02478223 0.05625226 0.021838162 0.1167618
2  0.10153485 0.07851149 0.08210976 0.009136116 0.3037535
3  0.08920885 0.03784539 0.08350120 0.031003551 0.1760821
4  0.11614286 0.03024431 0.11357654 0.064803399 0.1817103
5  0.58649314 0.31205208 0.53555263 0.139882174 1.3345694
6  0.60671504 0.13816254 0.59560073 0.362908849 0.9037621
7  0.82511310 0.66919682 0.65190496 0.075633463 2.5556600
8  0.82876692 0.67192202 0.65214217 0.072787714 2.5980293
9  1.50726343 0.72729109 1.39112058 0.446992670 3.2691806
10 1.96348289 0.41518719 1.93861873 1.237573701 2.8563891

$beta
      mean        sd   median      X2.5   X97.5
1 1.090288 0.4761821 1.009036 0.3956896 2.21233

GaryBAYLOR / R-code

Gibbs Sampler for Poisson distribution ((Robert and Casella, 10.17) #11