###################################################
### chunk number 1: 
###################################################
source("formats.R")
library(lattice)
options(width=65,show.signif.stars=FALSE)
options(SweaveHooks=list(fig=function() {
ltheme <- canonical.theme(color = FALSE)     ## in-built B&W theme
ltheme$strip.background$col <- "transparent" ## change strip bg
ltheme$plot.symbol$col <- "black"
ltheme$axis.text$cex <- 1.2
ltheme$par.xlab.text$cex <- 1.2
ltheme$par.ylab.text$cex <- 1.2
lattice.options(default.theme = ltheme)      ## set as default
}))


###################################################
### chunk number 2: 
###################################################
statlab <- read.csv("Data/statlab.csv")


###################################################
### chunk number 3: girls.hist eval=FALSE
###################################################
## # data(statlab,package="SMIR")
## girls <- subset(statlab,sex=="girl")
## library(lattice)
## print(histogram(~c.b.wgt,data=girls,ylab="Freq/0.6 units of birthweight"))


###################################################
### chunk number 4: 
###################################################
# data(statlab,package="SMIR")
girls <- subset(statlab,sex=="girl")
library(lattice)
print(histogram(~c.b.wgt,data=girls,ylab="Freq/0.6 units of birthweight"))


###################################################
### chunk number 5: 
###################################################
summary(girls$c.b.wgt)


###################################################
### chunk number 6: girls.hist eval=FALSE
###################################################
## print(histogram(~c.b.wgt,data=girls,breaks=seq(2.0,12,by=0.1),
## xlab="Birth weight in pounds", 
## main='Histogram of birthweight of girls'))


###################################################
### chunk number 7: 
###################################################
print(histogram(~c.b.wgt,data=girls,breaks=seq(2.0,12,by=0.1),
xlab="Birth weight in pounds", 
main='Histogram of birthweight of girls'))


###################################################
### chunk number 8: girls.ogive eval=FALSE
###################################################
## n <- with(girls,table(c.b.wgt))
## y <- as.numeric(names(n))
## cn <- cumsum(n)
## cn <- cn/(sum(n)+1)
## print(xyplot(cn~y,type='p',xlab='birthweight',ylab='cumulative proportion'))


###################################################
### chunk number 9: 
###################################################
n <- with(girls,table(c.b.wgt))
y <- as.numeric(names(n))
cn <- cumsum(n)
cn <- cn/(sum(n)+1)
print(xyplot(cn~y,type='p',xlab='birthweight',ylab='cumulative proportion'))


###################################################
### chunk number 10: girls.ccurve eval=FALSE
###################################################
## curve(pnorm(x,7.22,1.14),from=2.3,to=10.9,
## xlab='birthweight',ylab='cumulative proportion',las=1)
## points(y,cn)


###################################################
### chunk number 11: 
###################################################
curve(pnorm(x,7.22,1.14),from=2.3,to=10.9,
xlab='birthweight',ylab='cumulative proportion',las=1)
points(y,cn)


###################################################
### chunk number 12: girls.normal.deviate eval=FALSE
###################################################
## #nd <- qnorm(cn)
## curve(qnorm(pnorm(x,7.22,1.14)),from=2.3,to=10.9,
## xlab='birthweight',ylab='normal deviate',lwd=2)
## points(y,qnorm(cn),lwd=2)


###################################################
### chunk number 13: 
###################################################
#nd <- qnorm(cn)
curve(qnorm(pnorm(x,7.22,1.14)),from=2.3,to=10.9,
xlab='birthweight',ylab='normal deviate',lwd=2)
points(y,qnorm(cn),lwd=2)


###################################################
### chunk number 14: 
###################################################
source('NPL.bands.R')


###################################################
### chunk number 15:  eval=FALSE
###################################################
## library(SMIR)


###################################################
### chunk number 16: girls_cdf eval=FALSE
###################################################
## print(xyplot(lower+upper~x,
##              data=NPL.bands(girls$c.b.wgt),
##              pch=20,
## xlab='birthweight',ylab='cumulative proportion',
## panel=function(x,y,...){
## panel.xyplot(x,y,...)
## panel.curve(pnorm(x,7.22,1.14),lwd=1.5)}))#,from=min(x),to=max(x),lwd=2)
## #}))


###################################################
### chunk number 17: 
###################################################
print(xyplot(lower+upper~x,
             data=NPL.bands(girls$c.b.wgt),
             pch=20,
xlab='birthweight',ylab='cumulative proportion',
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.curve(pnorm(x,7.22,1.14),lwd=1.5)}))#,from=min(x),to=max(x),lwd=2)
#}))


###################################################
### chunk number 18: girls.bounds eval=FALSE
###################################################
## print(xyplot(qnorm(lower)+qnorm(upper)~x,
##              data=NPL.bands(girls$c.b.wgt),
##              pch=20,xlab="Birth weight (pounds)",
##              ylab="normal deviate", panel=function(x,y,...){
##                panel.xyplot(x,y,...)
##                panel.curve(qnorm(pnorm(x,7.22,1.14)),lwd=1.5)}))


###################################################
### chunk number 19: 
###################################################
print(xyplot(qnorm(lower)+qnorm(upper)~x,
             data=NPL.bands(girls$c.b.wgt),
             pch=20,xlab="Birth weight (pounds)",
             ylab="normal deviate", panel=function(x,y,...){
               panel.xyplot(x,y,...)
               panel.curve(qnorm(pnorm(x,7.22,1.14)),lwd=1.5)}))


###################################################
### chunk number 20: 
###################################################
env <-  NPL.bands(girls$c.b.wgt,conf.level=0.95)
print(xyplot(qnorm(lower)+qnorm(upper)~x,data=env,col=1,
pch=20,xlab="Birth weight (pounds)",
ylab="normal deviate", panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.curve(qnorm(0.09745*pnorm(x,9.132,0.8645)+
            0.885*pnorm(x,7.078,0.8645)+
            0.017*pnorm(x,3.892,0.8645)),
            from=min(x),to=max(x),lwd=2) }))


###################################################
### chunk number 21: boys_birthweight
###################################################
boys <- subset(statlab,sex=="boy")
minx <- min(boys$c.b.wgt)-0.1
maxx <- max(boys$c.b.wgt)+0.1
plota <- histogram(~c.b.wgt,data=boys,
                   ylab='Frequency/0.1 units of weight',
panel=function(x){
panel.histogram(x,breaks=seq(minx,maxx,by=0.1),type="count")
panel.text(12,35,"(a)",font=2)})
#
n <-  table(boys$c.b.wgt)
c <- cumsum(n) 
c <-  c/sum(n) 
y <- sort(unique(boys$c.b.wgt))
plotb <- xyplot(c~y,xlab="birthweight",pch=20,
                ylab="cumulative proportion",col="black",
                panel=function(x,y,...){
                  panel.xyplot(x,y,...)
                  panel.curve(pnorm(x,7.65,1.12),from=min(x),to=max(x))
                panel.text(6,0.8,"(b)",font=2)})
#
env <-  NPL.bands(boys$c.b.wgt,conf.level=0.95) 
plotc <- xyplot(qnorm(lower)+qnorm(upper)~x,data=env,
                xlab="birthweight",ylab="normal deviates",
                col="black",pch=20,
                panel=function(x,y,...){
                  panel.xyplot(x,y,...)
                  panel.curve(qnorm(pnorm(x,7.65,1.12)),
                              from=min(x),to=max(x),lwd=2)
                panel.text(6,3,"(c)",font=2)})
#
plotd <- xyplot(qnorm(lower)+qnorm(upper)~x,data=env,
                xlab="birthweight",ylab="normal deviates",
                col="black",pch=20,
                panel=function(x,y,...){
                  panel.xyplot(x,y,...)
                  panel.curve(qnorm(0.059*pnorm(x,6.96,1.082)+ 
                                    0.939*pnorm(x,7.685,1.082)+
                                    0.0016*pnorm(x,13.99,1.082)),
                              from=min(x),to=max(x),lwd=2)
                panel.text(6,3,"(d)",font=2)})
#
print(plota,position=c(0,0.49,0.51,1),more=TRUE)
print(plotb,position=c(0.49,0.49,1,1),more=TRUE)
print(plotc,position=c(0,0,0.51,0.51),more=TRUE)
print(plotd,position=c(0.49,0,1,0.51))


###################################################
### chunk number 22: hist.boy.mother.wgt eval=FALSE
###################################################
## boys <- subset(statlab,sex=="boy")
## library(lattice)
## print(histogram(~m.b.wgt,ylab='Freq/units of weight',
## main="Histogram of mother's weight at diagnosis of pregnancy, boys",
## breaks=seq(min(boys$m.b.wgt)-1,max(boys$m.b.wgt)+1,by=1),type="count",
## xlab="Mother's weight",data=boys))


###################################################
### chunk number 23: 
###################################################
boys <- subset(statlab,sex=="boy")
library(lattice)
print(histogram(~m.b.wgt,ylab='Freq/units of weight',
main="Histogram of mother's weight at diagnosis of pregnancy, boys",
breaks=seq(min(boys$m.b.wgt)-1,max(boys$m.b.wgt)+1,by=1),type="count",
xlab="Mother's weight",data=boys))


###################################################
### chunk number 24: boys.mum.wgt.mixture
###################################################
boys <- subset(statlab,sex=="boy")
library(npmlreg)
bmwgt.mix2 <- alldist(log(m.b.wgt)~1,data=boys,k=2,random.distribution="np",verbose=FALSE,plot.opt=0)


###################################################
### chunk number 25: mums_weight
###################################################
#boys <- subset(statlab,sex=="boy")
minx <- min(boys$m.b.wgt)-1
maxx <- max(boys$m.b.wgt)+1
meanx <- mean(boys$m.b.wgt)
sdx <- sd(boys$m.b.wgt)
meanlogx <- mean(log(boys$m.b.wgt))
sdlogx <- sd(log(boys$m.b.wgt))
#
n <-  table(boys$m.b.wgt)
c <- cumsum(n) 
c <-  c/sum(n) 
y <- as.numeric(names(n))
meanm1 <- coef(bmwgt.mix2)[1]
meanm2 <- coef(bmwgt.mix2)[2]
sdmix <- bmwgt.mix2$sdev$sdev
massm1 <- bmwgt.mix2$masses[1]
massm2 <- bmwgt.mix2$masses[2]
#
par(mfrow=c(2,2),mar=c(4,4,2,1.5))
#
plot(c~y,xlab="Mothers weight",pch=20,
                ylab="cumulative proportion",col="black")
curve(pnorm(x,meanx,sdx),lwd=1.5,add=TRUE)
text(200,0.6,"(a)",font=2)
#
#plota <- xyplot(c~y,xlab="Mothers weight",pch=20,
#                ylab="cumulative proportion",col="black",
#                panel=function(x,y,...){
#                  panel.xyplot(x,y,...)
#                  panel.curve(pnorm(x,meanx,sdx),lwd=1.5)
#                panel.text(200,0.8,"(a)",font=2)})
#
env <-  NPL.bands(boys$m.b.wgt,conf.level=0.95) 
plot(qnorm(lower)~x,data=env,ylim=c(-4,4),ylab="normal deviates",
xlab="mother's weight",col="black",pch=20)
points(qnorm(upper)~x,data=env,pch=20,col="black")
curve(qnorm(pnorm(x,meanx,sdx)) ,lwd=1.5,add=TRUE)
text(125,3,"(b)",font=2)
#
#plotb <- xyplot(qnorm(lower)+qnorm(upper)~x,data=env,
#                xlab="mother's weight",ylab="normal deviates",
#                col="black",pch=20,
#                panel=function(x,y,...){
#                  panel.xyplot(x,y,...)
#                  panel.curve(qnorm(pnorm(x,meanx,sdx)),
#                              from=min(x),to=max(x),lwd=1.5)
#                  panel.text(150,3,"(b)",font=2)})
#
env <-  NPL.bands(log(boys$m.b.wgt),conf.level=0.95) 
plot(qnorm(lower)~x,data=env,ylim=c(-4,4),
                xlab="mother's weight (log scale)",
                ylab="normal deviates",
                col="black",pch=20,axes=FALSE)
axis(1,at=log(c(100,150,200,250)),labels=c(100,150,200,250))
axis(2)
axis(3,at=log(c(100,150,200,250)),labels=FALSE)
axis(4,at=c(-4,-2,0,2,4),labels=FALSE)
box()
points(qnorm(upper)~x,data=env,pch=20)
curve(qnorm(pnorm(x,meanlogx,sdlogx)) ,lwd=1.5,add=TRUE)
text(log(125),3,"(c)",font=2)
#
#plotc <- xyplot(qnorm(lower)+qnorm(upper)~x,data=env,
#                xlab="mother's weight (log scale)",
#                ylab="normal deviates",
#                col="black",pch=20,
#                scales=list(x=list(log="e",
#                             at=seq(100,250,by=25),
#                             labels=seq(100,250,by=25),
#                              cex=0.8)),
#)),
#                panel=function(x,y,...){
#                  panel.xyplot(x,y,...)
#                  panel.curve(qnorm(pnorm(exp(x),meanx,sdx)),lwd=1.5)
#                  panel.text(log(125),3,"(c)",font=2)})
#
plot(qnorm(lower)~x,data=env,ylim=c(-4,4),
                xlab="mother's weight (log scale)",
                ylab="normal deviates",
                col="black",pch=20,axes=FALSE)
axis(1,at=log(c(100,150,200,250)),labels=c(100,150,200,250))
axis(2)
axis(3,at=log(c(100,150,200,250)),labels=FALSE)
axis(4,at=c(-4,-2,0,2,4),labels=FALSE)
box()
points(qnorm(upper)~x,data=env,pch=20,col="black")
curve(qnorm(massm1*pnorm(x,meanm1,sdmix)+ 
                                      massm2*pnorm(x,meanm2,sdmix)),lwd=1.5,add=TRUE)
text(log(125),3,"(d)",font=2)
#
#plotd <- xyplot(qnorm(lower)+qnorm(upper)~x,data=env,
#                xlab="mother's weight (log scale)",ylab="normal deviates",
#                col="black",pch=20,
#                scales=list(x=list(log="e",
#                                         at=seq(100,250,by=25),
#                                         labels=seq(100,250,by=25),
#                              cex=0.8)),
#                panel=function(x,y,...){
#                  panel.xyplot(x,y,...)
#                  panel.curve(qnorm(massm1*pnorm(log(x),meanm1,sdmix)+ 
#                                      massm2*pnorm(log(x),meanm2,sdmix)),lwd=1.5)
#                  panel.text(log(125),3,"(d)",font=2)})
#
#print(plota,position=c(0,0.49,0.51,1),more=TRUE)
#print(plotb,position=c(0.49,0.49,1,1),more=TRUE)
#print(plotc,position=c(0,0,0.51,0.51),more=TRUE)
#print(plotd,position=c(0.49,0,1,0.51))


###################################################
### chunk number 26: rlike eval=FALSE
###################################################
## c <- 0.146
## curve(exp(-12.5*(0.4-x)^2),from=-0.62,to=(-0.6 + 0.02*100),
## xlab='mean',
## ylab='relative likelihood')
## abline(h=c,lty='dotted')


###################################################
### chunk number 27: 
###################################################
c <- 0.146
curve(exp(-12.5*(0.4-x)^2),from=-0.62,to=(-0.6 + 0.02*100),
xlab='mean',
ylab='relative likelihood')
abline(h=c,lty='dotted')


###################################################
### chunk number 28: normal.t.like eval=FALSE
###################################################
## library(lattice)
## print(xyplot(0:1~-0.5:1.5,
## ylab='relative likelihoods',xlab='mean',type="n",
## panel=function(...){
## panel.xyplot(...)
## panel.curve((1+(5*(0.4-x))^2/24)^(-12.5))
## panel.curve(exp(-12.5*(0.4-x)^2),lty=2)
## panel.abline(h=c,lty=3)
## }))


###################################################
### chunk number 29: 
###################################################
library(lattice)
print(xyplot(0:1~-0.5:1.5,
ylab='relative likelihoods',xlab='mean',type="n",
panel=function(...){
panel.xyplot(...)
panel.curve((1+(5*(0.4-x))^2/24)^(-12.5))
panel.curve(exp(-12.5*(0.4-x)^2),lty=2)
panel.abline(h=c,lty=3)
}))


###################################################
### chunk number 30: profile.marginal.L eval=FALSE
###################################################
## #sig <- seq(0.515,length=100,by=0.015)
## #prl <- (0.98/sig)^25*exp(12.5*(1-(0.96/sig^2)))
## #mrl <- (1/sig)^24*exp(12*(1-(1/sig^2)))
## print(xyplot(0:1~0.6:2,type="n",
## xlab=expression(sigma),ylab='relative likelihood',
## panel=function(...){
## panel.curve((0.98/x)^25*exp(12.5*(1-(0.96/x^2))))
## panel.curve((1/x)^24*exp(12*(1-(1/x^2))),lty="dotted")
## panel.abline(h=0.146,lty='dotted')}))


###################################################
### chunk number 31: 
###################################################
#sig <- seq(0.515,length=100,by=0.015)
#prl <- (0.98/sig)^25*exp(12.5*(1-(0.96/sig^2)))
#mrl <- (1/sig)^24*exp(12*(1-(1/sig^2)))
print(xyplot(0:1~0.6:2,type="n",
xlab=expression(sigma),ylab='relative likelihood',
panel=function(...){
panel.curve((0.98/x)^25*exp(12.5*(1-(0.96/x^2))))
panel.curve((1/x)^24*exp(12*(1-(1/x^2))),lty="dotted")
panel.abline(h=0.146,lty='dotted')}))


###################################################
### chunk number 32: sigma.likelihood eval=FALSE
###################################################
## print(xyplot(0:1~0.6:2,type="n",xlab=expression(sigma),
## ylab='relative likelihood',
## panel=function(...){
## panel.xyplot(...)
## panel.curve((1/x)^24*exp(12*(1-(1/x^2))))
## panel.curve(exp(-0.5*(((x^(-2/3))-1)*(3*sqrt(12)))^2),lty="dotted",lwd=3)
## panel.abline(h=0.146,lty="dotted")}))


###################################################
### chunk number 33: 
###################################################
print(xyplot(0:1~0.6:2,type="n",xlab=expression(sigma),
ylab='relative likelihood',
panel=function(...){
panel.xyplot(...)
panel.curve((1/x)^24*exp(12*(1-(1/x^2))))
panel.curve(exp(-0.5*(((x^(-2/3))-1)*(3*sqrt(12)))^2),lty="dotted",lwd=3)
panel.abline(h=0.146,lty="dotted")}))


###################################################
### chunk number 34:  eval=FALSE
###################################################
## theta <- seq(from=0,to=1,by=0.001)
## bl <- theta*(1-theta)^9
## plot(theta,bl,xlab=expression(theta),ylab='likelihood',type='l')


###################################################
### chunk number 35:  eval=FALSE
###################################################
## phi <- pbeta(theta,1/3,1/3)
## plot(theta,phi,xlab=expression(theta),
## ylab=expression(phi),type='l')


###################################################
### chunk number 36: binomial eval=FALSE
###################################################
## theta <- seq(from=0,to=1,by=0.001)
## bl <- theta*(1-theta)^9
## plot.binomial.a <- xyplot(bl~theta, xlab=expression(theta),ylab='likelihood',type='l',main="(a)")
## phi <- pbeta(theta,1/3,1/3)
## plot.binomial.b <- xyplot(phi~theta,xlab=expression(theta),ylab=expression(phi),type='l',main="(b)")
## s <- 0.1^(-1/6)*0.9^(-1/6)/(5.3*sqrt(10))
## nl <- 0.1*0.9^9*exp(-0.5*((phi-0.267)/s)^2)
## plot.binomial.c <- xyplot(bl~phi,type='l',xlab=expression(phi),ylab='likelihood',main="(c)",
## panel=function(x,y,...){
## panel.xyplot(x,y,...)
## s <- 0.1^(-1/6)*0.9^(-1/6)/(5.3*sqrt(10))
## panel.curve(0.1*0.9^9*exp(-0.5*((x-0.267)/s)^2),lty='dotted',cex=1.5)})
## #llines(phi,nl,lty='dotted')})
## print(plot.binomial.a,position=c(0,0.5,0.5,1),more=TRUE)
## print(plot.binomial.b,position=c(0.5,0.5,1,1),more=TRUE)
## print(plot.binomial.c,position=c(0,0,0.5,0.5))


###################################################
### chunk number 37: 
###################################################
theta <- seq(from=0,to=1,by=0.001)
bl <- theta*(1-theta)^9
plot.binomial.a <- xyplot(bl~theta, xlab=expression(theta),ylab='likelihood',type='l',main="(a)")
phi <- pbeta(theta,1/3,1/3)
plot.binomial.b <- xyplot(phi~theta,xlab=expression(theta),ylab=expression(phi),type='l',main="(b)")
s <- 0.1^(-1/6)*0.9^(-1/6)/(5.3*sqrt(10))
nl <- 0.1*0.9^9*exp(-0.5*((phi-0.267)/s)^2)
plot.binomial.c <- xyplot(bl~phi,type='l',xlab=expression(phi),ylab='likelihood',main="(c)",
panel=function(x,y,...){
panel.xyplot(x,y,...)
s <- 0.1^(-1/6)*0.9^(-1/6)/(5.3*sqrt(10))
panel.curve(0.1*0.9^9*exp(-0.5*((x-0.267)/s)^2),lty='dotted',cex=1.5)})
#llines(phi,nl,lty='dotted')})
print(plot.binomial.a,position=c(0,0.5,0.5,1),more=TRUE)
print(plot.binomial.b,position=c(0.5,0.5,1,1),more=TRUE)
print(plot.binomial.c,position=c(0,0,0.5,0.5))


###################################################
### chunk number 38: binomial.likelihood.function
###################################################
theta <- seq(0,1,length=1000)
L <- dbinom(1,size=10,prob=theta)
maxL <- max(L)
L <- L/maxL
plot(theta,L,type='l',lwd=2,
ylab="Relative likelihood, binomial(1,10)",
xlab=expression(theta))
abline(h=0.146,lty=2,lwd=2)
a <- theta[L>0.146]
limits <- range(a)
round(limits,3)


###################################################
### chunk number 39: 
###################################################
uniroot(function(theta) 0.025-(1-pbinom(0,10,theta)),c(0,1))$root
uniroot(function(theta) 0.025-(pbinom(1,10,theta)),c(0,1))$root


###################################################
### chunk number 40: binomial.confidence.interval
###################################################
bin.conf.int <- function(r,n,alpha=0.05){
theta <- r/n
lambda <- qnorm(1-alpha/2)
ci2 <- lambda*sqrt(theta*(1-theta)/n+lambda^2/4/n/n)
lower <- (theta+lambda^2/2/n-ci2)/(1+lambda^2/n)
upper <- (theta+lambda^2/2/n+ci2)/(1+lambda^2/n)
list(lower=lower,upper=upper)
}


###################################################
### chunk number 41: 
###################################################
library(pscl)


###################################################
### chunk number 42: median.likelihood eval=FALSE
###################################################
## inc <- c(26,35,38,39,42,46,47,47,47,52,53,55,55,56,
##     58,60,60,60,60,60,65,65,67,67,69,70,71,72,
##     75,77,80,81,85,93,96,104,104,107,119,120)
## psi <- seq(40,90)
## like <- NULL
## for (i in seq(along=psi)){
## like[i] <- dbinom(sum(psi[i]>inc),length(inc),0.5)}
## print(xyplot(like~psi,type='s',xlab='median',ylab='likelihood'))


###################################################
### chunk number 43: Q3.likelihood eval=FALSE
###################################################
##  psi <- seq(60,110)
##  like <- NULL
##  for (i in seq(along=psi)){
##  like[i] <- dbinom(sum(psi[i]>inc),length(inc),0.75)
##  }
##  print(xyplot(like~psi,type='s',xlab='75-th percentile',ylab='likelihood'))


###################################################
### chunk number 44: percentile90.likelihood eval=FALSE
###################################################
## psi <- seq(80,130)
## like <- NULL
## for (i in seq(along=psi)){
## like[i] <- dbinom(sum(psi[i]>inc),length(inc),0.90)
## }
## print(xyplot(like~psi,type='s',xlab='90-th percentile',ylab='likelihood'))


###################################################
### chunk number 45: likelihoodmedianq3p90
###################################################
inc <- c(26,35,38,39,42,46,47,47,47,52,53,55,55,56,
    58,60,60,60,60,60,65,65,67,67,69,70,71,72,
    75,77,80,81,85,93,96,104,104,107,119,120)
psi <- seq(40,90)
like <- NULL
for (i in seq(along=psi)){
like[i] <- dbinom(sum(psi[i]>inc),length(inc),0.5)}
q2.plot <- xyplot(like~psi,type='s',xlab='median',ylab='likelihood',main="(a)")
#
 psi <- seq(60,110)
 like <- NULL
 for (i in seq(along=psi)){
 like[i] <- dbinom(sum(psi[i]>inc),length(inc),0.75)
 }
q3.plot <- xyplot(like~psi,type='s',xlab='75-th percentile',ylab='likelihood',main="(b)")
#
psi <- seq(80,130)
like <- NULL
for (i in seq(along=psi)){
like[i] <- dbinom(sum(psi[i]>inc),length(inc),0.90)
}
p90.plot <- xyplot(like~psi,type='s',xlab='90-th percentile',ylab='likelihood',main="(c)")
#
print(q2.plot,more=TRUE,position=c(0,0.49,0.51,1))
print(q3.plot,more=TRUE,position=c(0.49,0.49,1,1))
print(p90.plot,position=c(0,0,0.5,0.5))


###################################################
### chunk number 46: mean_empirical_likelihood
###################################################
source('el.s')
l <- NULL
mu <- seq(55,80,length=100)
for (i in seq(along=mu))l[i] <- exp(elm(inc,mu[i])$logelr)
print(xyplot(l~mu,type='l',xlab='mean',ylab='relative likelihood'))


###################################################
### chunk number 47: empirical_parametric_like eval=FALSE
###################################################
## lmu <- log(mu)
## sdseq <- seq(0.25,0.45,len=100)
## llnorm <- matrix(NA,nrow=100,ncol=100)
## for (i in 1:100){
##   for (j in 1:100){
##     llnorm[i,j] <- sum(dlnorm(inc,lmu[i]-var(log(inc))/2,sdseq[j],log=TRUE))
##   }
## }
## Llnorm <- exp(llnorm-max(llnorm))
## # contour(mu,sdseq,Llnorm)
## pll <- apply(Llnorm,1,max)
## #plot(exp(lmu),pl1,type="l")
## #
## #mu <- seq(55,85,len=100)
## scaleseq <- seq(4,10,len=100)
## lgamma <- matrix(NA,nrow=100,ncol=100)
## for (i in 1:100){
##   for (j in 1:100){
##     lgamma[i,j] <- sum(dgamma(inc,mu[i]/scaleseq[j],scale=scaleseq[j],log=TRUE))
##   }
## }
## Lgamma <- exp(lgamma-max(lgamma))
## # contour(mu,scaleseq,Lgamma)
## plg <- apply(Lgamma,1,max)
## #
## print(xyplot(l+pll+plg~mu,type="l",lty=1:3,ylab="relative likelihood",xlab="mean"))


###################################################
### chunk number 48: 
###################################################
lmu <- log(mu)
sdseq <- seq(0.25,0.45,len=100)
llnorm <- matrix(NA,nrow=100,ncol=100)
for (i in 1:100){
  for (j in 1:100){
    llnorm[i,j] <- sum(dlnorm(inc,lmu[i]-var(log(inc))/2,sdseq[j],log=TRUE))
  }
}
Llnorm <- exp(llnorm-max(llnorm))
# contour(mu,sdseq,Llnorm)
pll <- apply(Llnorm,1,max)
#plot(exp(lmu),pl1,type="l")
#
#mu <- seq(55,85,len=100)
scaleseq <- seq(4,10,len=100)
lgamma <- matrix(NA,nrow=100,ncol=100)
for (i in 1:100){
  for (j in 1:100){
    lgamma[i,j] <- sum(dgamma(inc,mu[i]/scaleseq[j],scale=scaleseq[j],log=TRUE))
  }
}
Lgamma <- exp(lgamma-max(lgamma))
# contour(mu,scaleseq,Lgamma)
plg <- apply(Lgamma,1,max)
#
print(xyplot(l+pll+plg~mu,type="l",lty=1:3,ylab="relative likelihood",xlab="mean"))