pumpgibbs <- function(failures, time, alpha, delta, gamma, niter, betastart) { # ************************************************************************* # # FILE: pumpgibbs.r # # Date Created: April 13, 2004 # Author: Mark E. Irwin # # Contents: Pump failure example of Gaver and O'Muircheartaigh (1987) as # proposed by Gelfand and Smith (1990). Analysis by Gibbs # Sampler. # # p_i | lambda_i ~ Poisson(lambda_i t_i) # lambda_i | beta ~ Gamma(alpha, beta) # beta | delta, gamma ~ IGamma(gamma, delta) # # Revision History # Date Name Changes/Reasons # # ************************************************************************* # # Input: # # failures: number of pump failures # time: observation times for the pumps # alpha: shape parameter for prior on lambda # delta: scale parameter for prior on beta # gamma: shape parameter for prior on beta # niter: number of Gibbs samples # betastart: starting value for beta, the scale parameter for prior on # lambda # # Output: # # lambda: simulated pump failure rates [niter x npumps] # beta: simulated scale parameter for prior on lambda [niter] # # ************************************************************************* # Initialize sampler variables npump <- length(failures) ntime <- length(time) if (npump != ntime) { stop('Number of pumps does not equal number of times.') } lambda <- matrix(0, nrow=niter, ncol=npump) beta <- rep(0, niter) betaold <- betastart # Run sampler for (iter in 1:niter) { lambdaold <- rgamma(npump, alpha + failures, scale = 1 / (time + 1 / betaold)) lambda[iter, ] <- lambdaold betainv <- rgamma(1, gamma + npump * alpha, scale = 1 / (sum(lambdaold) + delta)) betaold <- 1 / betainv beta[iter] <- betaold } list(lambda=lambda, beta=beta) }