library(MASS)
library(coda)
#Bayesian analysis of simulated data using a gamma likelihood function
#and Jeffreys prior on gamma dist parameters
#derivation of Jeffreys not shown
#simulate data from gamma(3,6) distibution
set.seed(848)
n <- 30
a <- 3
b <- 6
data<-rgamma(n, shape=a, rate=b)
prod = prod(data) #for liklihood ratios
sum = sum(data) #for liklihood ratios
data
cov = matrix(0,2,2)
cov[,1] = c(1,0)
cov[,2] = c(0,1)
#drive our parameters around a paramter space by sampling
#from a bivariate_normal(mu, cov) distribution
#parameters: mu->vector containing position of (a,b) parameters at step in MCMC
# cov-> var/cov matrix (assumed)
#return: new vector of parameters (a,b) to be proposed for next step of MCMC
driver <- function(mu, cov){
return(mvrnorm(n = 1, mu, cov))
}
#compute ratio of likelihood functions for two different sets of (a,b)
#assuming data comes from gamma dist
#Parameters: a,b -> current shape and rate parameters in MCMC
# propa, probb -> proposed shape and rate parameters
#returns likelihood ratio, used to calculate acceptance rate in MCMC sim
likehoodratio<-function(a, b, propa, propb){
return((((propb^propa)/(b^a))^n)*prod^(propa-a)*
((gamma(a)/gamma(propa))^n)*exp((b-propb)*sum))
}
#compute the prior ratio based on current and proposed values of a and b
#hard coded Jeffery's prior, could have used trigamma function instead
#Parameters: a,b -> current shape and rate parameters in MCMC
# propa, probb -> proposed shape and rate parameters
#returns pior ratio, used to calculate acceptance rate in MCMC sim
priorratio<-function(a,b,propa,propb){
numerator = 0
den = 0
for(i in 0:100){
numerator = numerator+(a)/((b^2)*(i+a)^2)
den = den+(propa)/((propb^2)*(i+propa)^2)
}
numerator = sqrt(numerator)
den = sqrt(den)
return(den/numerator)
}
#parent function to compute acceptance rate for MCMC
#that is to be compared to a uniform random variable on the interval 0,1
#Parameters: cura,curb -> current shape and rate parameters in MCMC
# propa, probb -> proposed shape and rate parameters
#returns alpha, used to dictate next step in MCMC sim
dist <- function(cura, curb, propa, propb){
l <- likehoodratio(cura, curb, propa, propb)
p <- priorratio(cura, curb, propa, propb)
return(l*p)
}
#metropolis-hastings algorithm
#Parameters: init -> initial starting point for MCMC sim (irrelevant)
# chain -> empty matrix to retain accepted vectors of (a, b)
#returns: matrix with posterior samples of f(a, b)
metro_hast <- function(init, chain){
acc = 0 #number of acceptances
for(i in 1:(length(chain)/2)){ #interate through empty nx2 matrix
prop <- driver(init, cov) #propose new values of a and b from driver
if(0<prop[1]&&0<prop[2]){ #values must be positive for gamma dist
p_1 <- dist(init[1], init[2], prop[1], prop[2]) #acceptance ratio
if(runif(1)<p_1){
init<-prop #accept proposed step with probability p_1
acc = acc+1 #note we accpepted something
chain[i,]<-init} #store values
}
}
print(acc/(length(chain)/2)) #print % of accepted steps
chain=na.omit(chain) #clear chain of NAs
return(chain) #return MCMC walk
}
init <- c(5,5) #initialize
chain <- matrix(NA, nrow = 100000, ncol = length(init))
mcmc <- metro_hast(init, chain) #run algorithm, accepted 32.401% of the time
##
## Iterations = 1:31402
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 31402
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## [1,] 4.095 0.978 0.005519 0.02401
## [2,] 7.682 1.927 0.010874 0.04915
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## var1 2.392 3.396 4.021 4.722 6.181
## var2 4.367 6.293 7.545 8.921 11.768
