#---------------------------------------------------------------------
#  1. Plots for Example 2.3
#---------------------------------------------------------------------




theta.vec <- seq(-10,10,length=20001)

plot(theta.vec,theta.vec^2,xlim=c(-3,3),ylim=c(0,3),type="l",xlab="Theta",ylab="Frequentist risk")
lines(theta.vec,.5^2+(1-.5)^2*theta.vec^2,lty=2,col=2)
lines(theta.vec,.8^2+(1-.8)^2*theta.vec^2,lty=3,col=3)
lines(theta.vec,1^2+(1-1)^2*theta.vec^2,lty=4,col=4)
lines(theta.vec,1.2^2+(1-1.2)^2*theta.vec^2,lty=5,col=5)
lines(theta.vec,1.3^2+(1-1.3)^2*theta.vec^2,lty=6,col=6)



par(mfcol=c(1,3))
theta.sd	<- c(0.2,1,5)

for(i in 1:3)
{
	theta.smp	<- rnorm(10^3,0,theta.sd[i])
	y.smp		<- theta.smp + rnorm(10^3,0,1)
 	plot(y.smp,theta.smp,xlab = "observation -- Y", ylab = "parameter -- theta",xlim = c(-2,2),ylim = c(-4,4))
	abline(0,1,col=3,lwd = 2)
	abline(0,theta.sd[i]^2/(1+theta.sd[i]^2),col=2,lwd = 2)
}




#---------------------------------------------------------------------
#  2. Hierarchical Bayesian model
#---------------------------------------------------------------------


library(LearnBayes)
data(cancermortality)

#?cancermortality


#	2.0 Function that generates negtive-binomial samples

summary(rbeta(10^5,2,3))
summary(rbinom(10^5,20,rbeta(10^5,2,3)))

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



#	2.1 Scatter plot of data

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

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 = qbeta(.025,.5 + sum(y.vec),.5 + sum(n.vec) - sum(y.vec)),col=2,lty=2,lwd = 2)
abline(h = qbeta(.975,.5 + sum(y.vec),.5 + sum(n.vec) - sum(y.vec)),col=2,lty=2,lwd = 2)
abline(h = (.5 + sum(y.vec)) / (1+ sum(n.vec)),col=2,lty=1,lwd = 1)


#	2.2 Overlay with predictive 95% CI



y.qt.025	<- rep(NA,20)
y.qt.975	<- rep(NA,20)

for(i in 1:20)
{
	y.pred.smp	<- rbetabinom.ab(10^4,n.vec[i], .5 + sum(y.vec[-i]),.5 + sum(n.vec[-i]) - sum(y.vec[-i]))
	y.qt.025[i]	<- quantile(y.pred.smp,0.025)
	y.qt.975[i]	<- quantile(y.pred.smp,0.975)
}

segments(n.vec,y.qt.975/n.vec,n.vec,y.qt.025/n.vec)


#	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

mycontour(betabinexch,c(-8,-4.5,3,16.5),cancermortality)
title(xlab="logit eta",ylab = "log K")


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)

points(theta.hat.1,theta.hat.2,pch = "+",col=4,cex = 2)


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


(fit <- laplace(betabinexch,c(-7,6),cancermortality))
points(fit$mode[1],fit$mode[2],pch=3)

npar <- list(m=fit$mode,v=fit$var)
mycontour(lbinorm,c(-8,-4.5,3,16.5),npar,add = T,col=2)
title(xlab="logit eta", ylab="log K")

(a <- fit$mode[1] + 1.96 * sqrt(fit$var[1,1]))
(b <- fit$mode[1] - 1.96 * sqrt(fit$var[1,1]))

exp(a) / (1 + exp(a))
exp(b) / (1 + exp(b))


(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 posterior means of the rates for the beta-binomial model

n.smp	<- dim(theta)[1]

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 = qbeta(.025,.5 + sum(y.vec),.5 + sum(n.vec) - sum(y.vec)),col=2,lty=2,lwd = 2)
abline(h = qbeta(.975,.5 + sum(y.vec),.5 + sum(n.vec) - sum(y.vec)),col=2,lty=2,lwd = 2)
abline(h = (.5 + sum(y.vec)) / (1+ sum(n.vec)),col=2,lty=1,lwd = 1)

p.post.mean	<- rep(0,20)

for(i in 1:n.smp)
{	
	eta.smp		<- exp(theta[i,1]) / (1+exp(theta[i,1]))
	K.smp			<- exp(theta[i,2])
	alpha.smp		<- eta.smp*K.smp
	beta.smp		<- (1-eta.smp)*K.smp
	p.post.mean		<-  p.post.mean + ((alpha.smp + y.vec) / (alpha.smp + beta.smp + n.vec)) / n.smp
}


points(n.vec,p.post.mean,,pch = "+",col=3,cex = 1.5)
segments(n.vec,y.vec/n.vec,n.vec,p.post.mean,col = 5)