PoisGam <- function(sx, n, alpha, beta, thres, quant, xmax){ # ***************************************************************************** # Function to (1) Estimate posterior quantities for the poisson-gamma model, # (2) Plot likelihood, prior, and posterior. # # Inputs: # sx: \sum_{i=1}^{n}x_{i} # n: number of observations # alpha, beta: shape/rate parameters for gamma prior on lambda # thres: threshod for cumulative probability # quant: quantiles to compute # xmax: max x-axis value for lambda (for plot) # # Outputs: # post.mean: Posterior mean of lambda # post.quant: Posterior quantiles # cumprob: Cumulative posterior prob that lambda <= thres # post.samples: Posterior samples of lambda # XfGivenTheta: Posterior predictive samples # ***************************************************************************** # Plot likelihood (scaled) lambda <- seq(0,xmax,length=1000) plot(lambda, (10^17.4)* exp(-n*lambda) * (lambda^sx), type = "l", col = "black", xlab = expression(lambda), ylab = "probability (scaled)", xlim = c(0, xmax), main = paste("a=", alpha, "; b=", beta)) # Plot prior lines(lambda, dgamma(lambda, alpha, beta), col = "black", lty = 3) # Plot posterior lines(lambda, dgamma(lambda, alpha+sx, n+beta), col = "red") # Compute posterior mean post.mean <- (alpha + sx)/ (n + beta) # Compute quantiles post.quant <- qgamma(quant, alpha+sx, n+beta) abline(v=post.quant, col = "blue") # Compute cumulative probability cumprob <- pgamma(thres, alpha+sx, n+beta) abline(v=thres, col = "green") # Draw posterior samples post.samples <- rgamma(10000, alpha+sx, n+beta) post.mean.aprox <- mean(post.samples) post.quant.approx <- quantile(post.samples, quant) cumprob.approx <- length(which(post.samples<=thres)) / length(post.samples) # Draw posterior predictive samples using MC approximation ThetaGivenX <- rgamma(10000, alpha+sx, n+beta) XfGivenTheta <- rpois(10000, ThetaGivenX) legend("topright", c("Likelihood", "Prior", "Posterior", "Quantiles", "Cumulative Prob. Threshold"), col = c("black", "black", "red", "blue", "green"), lty = c(1,3,1,1,1)) return(list(post.mean = post.mean, post.quant = post.quant, cumprob = cumprob, post.samples = post.samples, XfGivenTheta = XfGivenTheta)) }