lambdarw <- function(lambdaold, failures, time, mu, sigma2, rwsd) { # ************************************************************************* # # FILE: lambdarw.r # # Date Created: April 20, 2004 # Author: Mark E. Irwin # # Contents: Modification of the Pump failure example of Gaver and O'Muircheartaigh # (1987) as proposed by Gelfand and Smith (1990). Analysis by Metropolis # within Gibbs. # # p_i | lambda_i ~ Poisson(lambda_i t_i) # lambda_i | mu, sigma2 ~ logN(mu, sigma2) # mu | nu, tau ~ N(nu, tau2) # sigma2 | delta, gamma ~ IGamma(gamma, delta) # # Samples new lambda by M-H with a random walk proposal on log lambda_i # # Revision History # Date Name Changes/Reasons # # ************************************************************************* # # Input: # # lambdaold: current lambda # failures: number of pump failures # time: observation times for the pumps # mu: mean parameter for prior on log(lambda) # sigma2: variance for prior on log(lambda) # rwsd: standard deviation of random walk on log lambda for proposal distribution # # Output: # # lambda: simulated pump failure rates [npumps] # # ************************************************************************* # Initialize sampler variables npump <- length(failures) ntime <- length(time) if (npump != ntime) { stop('Number of pumps does not equal number of times.') } sigma <- sqrt(sigma2) lambda <- lambdaold for (i in 1:npump) { # Generate proposal prop <- rlnorm(1, log(lambdaold[i]), rwsd) # calculate acceptance probability HR <- dpois(failures[i], prop * time[i]) / dpois(failures[i], lambdaold[i] * time[i]) HR <- HR * dlnorm(prop, mu, sigma) / dlnorm(lambdaold[i], mu, sigma) HR <- HR * dlnorm(lambdaold[i], log(prop), rwsd) / dlnorm(prop, log(lambdaold[i]), rwsd) acc <- min(HR, 1) # Accept proposal? u <- runif(1) lambda[i] <- ifelse(u <= acc, prop, lambdaold[i]) } lambda }