# PLOTS IN LECTURE NOTES AND EXERCISES

# ==========================================================================================

# One proportion examples from CLAYTON & HILLS:
# =============================================

# Likelihood
#====================
x=c(1,1,0,1,0,0,0,1,0,0)
m=mean(x)
p=seq(m-0.3,m+0.3,0.01)
L<-rep(0,length(p))
for (i in 1:length(p)) L[i]=prod(dbinom(x,1,p[i]))
plot(p,L,type="l")
abline(v=0.4)

# Likelihood ratio (compared to most likely)
#========================================
x=c(1,1,0,1,0,0,0,1,0,0)
m=mean(x)
p=seq(m-0.3,m+0.3,0.001)
LR<-rep(0,length(p))
ml<-prod(dbinom(x,1,0.4))
for (i in 1:length(p)) LR[i]=prod(dbinom(x,1,p[i]))/ml
plot(p,LR,type="l")
abline(v=0.4,h=0.258)
arrows(0.178,0.258,0.178,0)
arrows(0.655,0.258,0.655,0)

# Likelihood ratio for large sample size
#========================================
x=rep(c(1,0),c(20,30))
m=mean(x)
p=seq(m-0.3,m+0.3,0.001)
LR<-rep(0,length(p))
ml<-prod(dbinom(x,1,0.4))
for (i in 1:length(p)) LR[i]=prod(dbinom(x,1,p[i]))/ml
plot(p,LR,type="l")
abline(v=0.4,h=0.258)
arrows(0.291,0.258,0.291,0)
arrows(0.516,0.258,0.516,0)

# Log likelihood
#================
x=c(1,1,0,1,0,0,0,1,0,0)
m=mean(x)
p=seq(m-0.3,m+0.3,0.01)
l<-rep(0,length(p))
ml<-prod(dbinom(x,1,0.4))
for (i in 1:length(p)) l[i]=sum(log(dbinom(x,1,p[i])))-log(ml)
plot(p,l,type="l",ylab="l")
abline(v=0.4)
abline(v=0.4,h=-1.353)
arrows(0.178,-1.353,0.178,-3)
arrows(0.655,-1.353,0.655,-3)

# ==========================================================================================

library(lattice)

# Exercise 1 on binned data
#====================

# 2-d version
mu<-seq(-10,100,1)
sigma=seq(0.1,10,0.5)
g<-expand.grid(mu=mu,sigma=sigma)
count<-0
for (i in 1:length(mu)) { 
for (j in 1:length(sigma)) { 
count<-count+1
g$L[count]<-(pnorm(4,mu[i],sigma[j])-pnorm(3,mu[i],sigma[j]))*(pnorm(6,mu[i],sigma[j])-pnorm(4,mu[i],sigma[j]))*(1-pnorm(5,mu[i],sigma[j]))*(pnorm(3.6,mu[i],sigma[j]))
} }
g[g$L==max(g$L),]
wireframe(L~mu*sigma,data=g,drape = TRUE, colorkey = TRUE,scale=list(arrows=FALSE),abline=g[g$L==max(g$L),])


# tumbler mortality:
#====================

# 2-d solution
a<-read.table("x://teaching//bios602 (2009)//tumbler.txt",header=T)
rate<-seq(0.01,0.3,0.01)
shape<-seq(0.1,3,0.1)
g<-expand.grid(shape=shape,rate=rate)
count<-0
for (i in 1:length(shape)) { 
for (j in 1:length(rate)) { 
count<-count+1
p<-(pgamma(a$x,shape[i],rate[j])-pgamma(a$x-1,shape[i],rate[j]))/(1-pgamma(a$x-1,shape[i],rate[j]))
g$L[count]<-sum(a$n*log(p)+(a$N-a$n)*log(1-p))
}
}
g[g$L==max(g$L),]
wireframe(L~shape+rate,data=g,drape = TRUE, colorkey = TRUE,scale=list(arrows=FALSE),abline=g[g$L==max(g$L),])