#---------------------------------------------------------------------
#  Script for solution to question 1
#---------------------------------------------------------------------



p.seq	<- seq(0,1,length = 100)

binom.gep.pval 	<- 1 - pbinom(8,12,p.seq)
binom.lep.pval 	<- pbinom(9,12,p.seq)

nbinom.gep.pval 	<- 1 - pnbinom(8,3,1-p.seq)
nbinom.lep.pval 	<- pnbinom(9,3,1-p.seq)

plot(p.seq,binom.gep.pval,type = "l",col=2,lwd = 2)
lines(p.seq,binom.lep.pval,col=3,lwd = 2)

lines(p.seq,nbinom.gep.pval,type = "l",col=2,lwd = 2)
lines(p.seq,nbinom.lep.pval,col=3,lwd = 2)
abline(h = 0.025,lty = 2)


approx(




#---------------------------------------------------------------------
#  Script for solution to question 5
#---------------------------------------------------------------------


library(LearnBayes)
data(cancermortality)

rbetabinom.ab <- function(n,size,a,b)   rbinom(n,size,rbeta(n,a,b))

y.vec    <- cancermortality$y
n.vec    <- cancermortality$n



#	2.2  Use learnbayes package functions to plot posterior probability of (eta, K)


# Beta-binomial model
#  y ~ binom(n,p)
#  p ~ beta(alpha,beta)

#  Parameterization
#  mean      --  eta =  alpha / (alpha + beta)     
#  precision --  K   =  alpha + beta 

# Prior
#	g(eta,K) \propto (eta^-1 *(1-eta)^-1 *(1+K)^-2



alpha.hat	<-  .5 + sum(y.vec)
beta.hat	<- .5 + sum(n.vec) - sum(y.vec)
eta.hat	<-  alpha.hat / (alpha.hat + beta.hat)
K.hat		<-  alpha.hat + beta.hat
theta.hat.1	<- log(eta.hat/(1-eta.hat))
theta.hat.2	<- log(K.hat)



#  2.3 Use learnbayes functions for rejection sampling from  posterior probability of (eta, K)


(fit <- laplace(betabinexch,c(-7,6),cancermortality))

(eta.mode	<- exp(fit$mode[1]) / (1+exp(fit$mode[1])))
(K.mode	<- exp(fit$mode[2]))

eta.mode * K.mode		#	a - parameter estimate
(1-eta.mode) * K.mode	#	b - parameter estimate


# Rejection Sampling

betabinT=function(theta,datapar)
{
	data <- datapar$data
	tpar <- datapar$par
	d    <- betabinexch(theta,data)-dmt(theta,mean=c(tpar$m),
							  S=tpar$var,df=tpar$df,log=TRUE)
return(d)
}


tpar    <- list(m=fit$mode,var=2*fit$var,df=4)
datapar <- list(data=cancermortality,par=tpar)
datapar <- list(par=tpar,data=cancermortality)
start <- c(-7,12)
fit1  <- laplace(betabinT,start,datapar)

fit1$mode
betabinT(fit1$mode,datapar)

theta <- rejectsampling(betabinexch,tpar,-569.28,10000,cancermortality)
dim(theta)

mycontour(betabinexch,c(-8,-4.5,3,16.5),cancermortality)
title(xlab="logit eta",ylab="log K")
points(theta[,1],theta[,2])


#	Find credible intervals for the city-cancer-rate (the p parameter) for the beta-binomial model


eta.smp		<- exp(theta[,1]) / (1+exp(theta[,1]))
K.smp			<- exp(theta[,2])
alpha.smp		<- eta.smp*K.smp
beta.smp		<- (1-eta.smp)*K.smp
n.smp			<- length(beta.smp)

cred.int	<- array(dim=c(2,21))

for(i in 1:20)
	cred.int[,i]	<- quantile(rbeta(n.smp, alpha.smp + y.vec[i],beta.smp + n.vec[i]-y.vec[i]),c(0.025,0.975)) 
	
cred.int[,21]	<- quantile(rbeta(n.smp, alpha.smp,beta.smp),c(0.025,0.975)) 


par(mfcol=c(1,1))
plot(n.vec, y.vec / n.vec,log="x",xlab = "Number of patients at risk", ylab = "Cancer rate",pch = "+",col=4,cex = 1.5)
abline(h = (.5 + sum(y.vec)) / (1+ sum(n.vec)),col=2,lty=1,lwd = 1)
abline(h = cred.int[1,21],col=2,lty=1,lwd = 1)
abline(h = cred.int[2,21],col=2,lty=1,lwd = 1)
segments(n.vec,cred.int[1,1:20],n.vec,cred.int[2,1:20],col = 5)