####################
#	Question 1
####################


#	Accident data of Simar (1976) -- Example 3.2.3, Carlin and Louis page 60

y.vec		<- 0:7
count.vec	<- c(7840,1317,239,42,14,4,4,1)

cbind(y.vec,count.vec, round(count.vec/sum(count.vec),5))

(mn.y	<- sum(y.vec * count.vec) / sum(count.vec))
(var.y	<-sum((y.vec-mn.y)^2 * count.vec) / sum(count.vec))

#	This yields method of moments estimates for r and theta

(theta.MM	<- var.y/ mn.y - 1)
(r.MM		<- mn.y/theta.MM)


#	For lambda ~ Gamma( shape = r, scale = theta) Y is NB(size = r, prob = 1 / (1+theta))


(loglik.MM	<- sum(count.vec*dnbinom(y.vec, size = r.MM, prob = 1 / (1+theta.MM), log = T)))


r.vec		<- seq(0.55,1,length = 400)
theta.vec	<- seq(0.2,0.4,length = 400)
loglik.mat	<- array(dim=c(400,400))
for(i in 1:400)
for(j in 1:400)
	loglik.mat[i,j]	<- sum(count.vec*dnbinom(y.vec, size = r.vec[i], prob = 1 / (1+theta.vec[j]), log = T))


#	Draw likelihood contours

contour(r.vec,theta.vec,max(loglik.mat) - loglik.mat,levels = qchisq(c(0.5,0.90,0.95),2)/2,xlab = "r",ylab = "theta")
points(r.MM,theta.MM,col=2)


#	Find MLE

minus.loglik.theta.r	<- function(theta.r.vec) -sum(count.vec*dnbinom(y.vec, size = theta.r.vec[2], 
											prob = 1 / (1+theta.r.vec[1]), log = T))



minus.loglik.theta.r(c(theta.MM,r.MM))
(a	<- nlminb(c(theta.MM,r.MM), minus.loglik.theta.r))
points(a$par[2],a$par[1],col=3)


#	Are the CI's correct? (1. Study test statistic distribution)

theta.tst	<- 0.30
r.tst		<- 0.70
n		<- sum(count.vec)
y.smp		<- rnbinom(n, size = r.tst,prob = 1 / (1+theta.tst))
table(y.smp)

minus.loglik.theta.r.y.smp	<- function(theta.r.vec) - sum(dnbinom(y.smp, size = theta.r.vec[2], prob = 1 / (1+theta.r.vec[1]), log = T))
minus.loglik.theta.r.y.smp(c(theta.tst,r.tst))

p.vec	<- rep(NA,100)
for(i in 1:100)
{
	y.smp			<- rnbinom(n, size = r.tst,prob = 1 / (1+theta.tst))
	minus.max.loglik	<- nlminb(c(theta.tst,r.tst), minus.loglik.theta.r.y.smp)$objective
	minus.tst.loglik	<- minus.loglik.theta.r.y.smp(c(theta.tst,r.tst))
	p.vec[i]			<- 1 - pchisq(2 * (minus.tst.loglik - minus.max.loglik),2)
}
plot(sort(p.vec))


#	Are the CI's correct? (2. Look at size of likelihood -- goodness of fit tests can only tell us that model isn't good!!)


nlminb(c(theta.MM,r.MM), minus.loglik.theta.r)	# max NPML yielded max loglik  =  - 5340.74 






####################
#	Question 2
####################


#	From sikum 3: Bayes risk expression is (n - (n-2)^2 * E[1/S]) for S ~ (1+tau^2_0)* chisquared RV with n df 
#	Which is equal:  n - (n-2)/(1 + tau^2_0)

#tau.0.sq	<- 3
tau.0.sq	<- 0.03
n		<- 10

mu.mat	<- array(rnorm(10^4 * n ,mean = 0,sd = sqrt(tau.0.sq)),dim=c(10^4,n))
y.mat		<- mu.mat	+ rnorm(10^4*n,mean = 0,sd = 1)
s.vec		<- apply(y.mat^2,1,sum)

y.MLE		<-  y.mat
y.JS		<- (1 - (n-2)/s.vec) * y.mat
y.Bayes	<-  (tau.0.sq/(1+tau.0.sq)) * y.mat


risk.MLE	<- apply((mu.mat - y.MLE)^2,1,sum)
risk.JS	<- apply((mu.mat - y.JS)^2,1,sum)
risk.Bayes	<- apply((mu.mat - y.Bayes)^2,1,sum)

t.test(risk.MLE)			;   n 					#  = n
t.test(risk.JS)		 	;   n - (n-2)/(1+tau.0.sq)		# Between n and 2
t.test(risk.Bayes)		;   n * (tau.0.sq/(1+tau.0.sq))	# Between n and 0