# 1c dat <- read.csv("growth.csv")[,2] # alpha*lambda^alpha/(x+lambda)^(alpha+1) # L() = alpha^n*lambda^(alpha*n)/prod(x+lamda)^(alpha+1) ll <- function(dat,alpha,lambda){ n <- length(dat) n*log(alpha)+n*alpha*log(lambda)-(alpha+1)*sum(log(dat+lambda)) } ll(dat,3,3) lprior <- function(x)dexp(x,1/3,log=TRUE) nrun <- 11000 nburn <- 1000 alphavec <- numeric(nrun) lambdavec <- numeric(nrun) accept <- 0 prop_cov <- matrix(c(.7,1.03,1.03,1.85),2,2)/50 library(mnormt) alphavec[1] <- 1 lambdavec[1] <- 1 for(i in 2:nrun){ parsnew <- exp(rmnorm(1,t(log(c(alphavec[i-1],lambdavec[i-1]))),prop_cov)) alph <- ll(dat,parsnew[1],parsnew[2])+sum(lprior(parsnew))+sum(log(parsnew))- ll(dat,alphavec[i-1],lambdavec[i-1])-lprior(alphavec[i-1])-lprior(lambdavec[i-1])-log(alphavec[i-1])-log(lambdavec[i-1]) alphavec[i] <- alphavec[i-1] lambdavec[i] <- lambdavec[i-1] if(log(runif(1)) < alph){ alphavec[i] <- parsnew[1] lambdavec[i] <- parsnew[2] accept <- accept + 1 } } accept/nrun plot.ts(alphavec) plot.ts(lambdavec) acf(alphavec) acf(lambdavec) cov(log(cbind(alphavec,lambdavec))) mean(alphavec[(nburn+1):nrun]) mean(lambdavec[(nburn+1):nrun]) # Credible Intervals quantile(alphavec[(nburn+1):nrun],c(.025,.975)) quantile(lambdavec[(nburn+1):nrun],c(.025,.975))