###################################################
### chunk number 1: 
###################################################
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
#load(".Rchapter7")
}))


###################################################
### chunk number 2: 
###################################################
statlab <- read.csv('Data/statlab.csv')
#data(statlab,package="SMIR")
girls <- subset(statlab,sex=="girl")
n <- with(girls,table(c.b.wgt))
y <- as.numeric(names(n))
cn <- cumsum(n)
cn <- cn/(sum(n)+1)
print(xyplot(qnorm(cn)~y,type='p',xlab='birthweight',
ylab='normal deviate',col="black",
ylim=c(-4.5,3.5),
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.curve(qnorm(pnorm(x,7.22,1.14)))}))


###################################################
### chunk number 3:  eval=FALSE
###################################################
## data(statlab,package="SMIR")


###################################################
### chunk number 4: 
###################################################
statlab <- read.csv("Data/statlab.csv",header=TRUE)


###################################################
### chunk number 5: girls.k1
###################################################
girls <- subset(statlab,sex=="girl")
library(npmlreg)
girls.k1 <- alldist(c.b.wgt~1,data=girls,k=1)
summary(girls.k1)
round(girls.k1$disparity,2)


###################################################
### chunk number 6: girls.np eval=FALSE
###################################################
## girls.k2 <- update(girls.k1,k=2,plot.opt=0,verbose=FALSE)
## girls.k3 <- update(girls.k2,k=3)
## girls.k4 <- update(girls.k3,k=4)
## girls.k5 <- update(girls.k4,k=5)
## girls.k6 <- update(girls.k5,k=6)


###################################################
### chunk number 7: girls.np.summaries
###################################################
summary(girls.k2)
round(girls.k2$disparity,2)
summary(girls.k3)
round(girls.k3$disparity,2)
summary(girls.k4)
round(girls.k4$disparity,2)
#summary(girls.k5)
#summary(girls.k6)


###################################################
### chunk number 8: lrts4 eval=FALSE
###################################################
## lrts <- function(mu,sd,p,N,k1){
## k <- length(p)
## sigma <- ifelse(length(sd)==1,rep(sd,k),sd)
## n1 <- round(N*p[-k])
## nk <- N-sum(n1)
## n <- c(n1,nk)
## y <- NULL
## for(i in 1:k) y <- c(y,rnorm(n[i],mu[i],sigma[i]))
## data <- data.frame(yy=y)
## fit1 <- alldist(yy~1,k=k,data=data,plot.opt=0,verbose=FALSE)
## fit2 <- update(fit1,k=k1)
## fit1$disparity-fit2$disparity
## }
## lrts3vs4 <- NULL
## for (R in 1:19){
## lrts3vs4 <- c(lrts3vs4,lrts(mu=c(3.892,7.078,9.132),sd=0.8645,p=c(0.017,0.885,0.0974),N=648,k1=4))
## }


###################################################
### chunk number 9: 
###################################################
round(sort(lrts3vs4),4)


###################################################
### chunk number 10: single_vs_three_mixture eval=FALSE
###################################################
## lrts1vs3 <- NULL
## for (R in 1:19){
## lrts1vs3 <- c(lrts1vs3,lrts(mu=7.223,sd=1.146,p=1,N=648,k1=3))
## }


###################################################
### chunk number 11: 
###################################################
sort(round(lrts1vs3,2))


###################################################
### chunk number 12: girl.regk1
###################################################
girl.regk1 <- alldist(c.b.wgt~m.b.ag+m.b.wgt,data=girls,k=1,plot.opt=0,verbose=FALSE)
summary(girl.regk1)
round(girl.regk1$disparity,2)


###################################################
### chunk number 13: girl.regk4
###################################################
girl.regk1a <- update(girl.regk1,.~.-m.b.ag)
summary(girl.regk1a)
round(girl.regk1a$disparity,2)
girl.regk3a <- update(girl.regk1a,k=3)
summary(girl.regk3a)
round(girl.regk3a$disparity,2)


###################################################
### chunk number 14: galaxy eval=FALSE
###################################################
## data(galaxies,data="SMIR")


###################################################
### chunk number 15: 
###################################################
data(galaxies,package="MASS")


###################################################
### chunk number 16: galaxy.plot eval=FALSE
###################################################
## galaxies[78] <- 26960
## velocity <- galaxies/1000
## #galaxies.df <- data.frame(velocity=galaxies/1000)
## print(xyplot(lower+upper~x,data=NPL.bands(velocity),pch=20, col=1,
## xlab="velocity",ylab="probability",
## panel=function(...){
## panel.xyplot(...)
## mgalaxy <- mean(velocity)
## sdgalaxy <- sd(velocity)
## panel.curve(pnorm(x,mgalaxy,sdgalaxy),lwd=2)
## }))


###################################################
### chunk number 17: 
###################################################
galaxies[78] <- 26960
velocity <- galaxies/1000
#galaxies.df <- data.frame(velocity=galaxies/1000)
print(xyplot(lower+upper~x,data=NPL.bands(velocity),pch=20, col=1,
xlab="velocity",ylab="probability",
panel=function(...){
panel.xyplot(...)
mgalaxy <- mean(velocity)
sdgalaxy <- sd(velocity)
panel.curve(pnorm(x,mgalaxy,sdgalaxy),lwd=2)
}))


###################################################
### chunk number 18: galaxy_mixture eval=FALSE
###################################################
## #<<galaxy.np,echo=FALSE>>=
## library(npmlreg)
## galaxies.df <- data.frame(velocity=velocity)
## galaxy.k1 <- alldist(velocity~1,data=galaxies.df,k=1,plot.opt=0,lambda=0,verbose=FALSE)
## galaxy.k2 <- update(galaxy.k1,k=2)


###################################################
### chunk number 19:  eval=FALSE
###################################################
## galaxy.k8 <- update(galaxy.k1,k=8)#,tol=0.2)
## #galaxy.k9 <- update(galaxy.k1,k=9,tol=0.06)
## #galaxy.k10 <- update(galaxy.k1,k=10)
## #
## galaxy.k2u <- update(galaxy.k2,lambda=1,spike.protect=1)
## galaxy.k3u <- update(galaxy.k2u,k=3)


###################################################
### chunk number 20:  eval=FALSE
###################################################
## galaxy.k8u <- update(galaxy.k2u,k=8,lambda=0.99)#tol=0.5,spike.protect=TRUE)
## #galaxy.k9u <- update(galaxy.k2u,k=9)


###################################################
### chunk number 21: galaxy.np eval=FALSE
###################################################
## library(npmlreg)
## galaxy.k1 <- alldist(velocity~1,data=galaxies.df,k=1,plot.opt=0,lambda=0,verbose=FALSE)
## galaxy.k2 <- update(galaxy.k1,k=2)
## galaxy.k3 <- update(galaxy.k1,k=3)
## galaxy.k4 <- update(galaxy.k1,k=4)
## galaxy.k5 <- update(galaxy.k1,k=5)
## galaxy.k6 <- update(galaxy.k1,k=6,tol=0.2)
## galaxy.k7 <- update(galaxy.k1,k=7,tol=0.12)
## galaxy.k8 <- update(galaxy.k1,k=8,tol=0.2)
## #galaxy.k9 <- update(galaxy.k1,k=9,tol=0.06)
## #galaxy.k10 <- update(galaxy.k1,k=10)
## #
## galaxy.k2u <- update(galaxy.k2,lambda=1,spike.protect=1)
## galaxy.k3u <- update(galaxy.k2u,k=3,tol=0.5,lambda=0.99)
## galaxy.k4u <- update(galaxy.k2u,k=4,tol=0.5,lambda=0.99)
## galaxy.k5u <- update(galaxy.k2u,k=5,tol=0.5,lambda=0.99)
## galaxy.k6u <- update(galaxy.k5u,k=6)
## galaxy.k7u <- update(galaxy.k2u,k=7)
## galaxy.k8u <- update(galaxy.k2u,k=8,lambda=0.999,tol=0.5,spike.protect=TRUE)
## #galaxy.k9u <- update(galaxy.k2u,k=9)
## #galaxy.k10u <- update(galaxy.k2u,k=10)


###################################################
### chunk number 22: 
###################################################
model.list <- list(galaxy.k1,galaxy.k2,galaxy.k3,galaxy.k4,galaxy.k5,
galaxy.k6,galaxy.k7,galaxy.k8)#,galaxy.k9,galaxy.k10)
modelu.list <- list(galaxy.k1,galaxy.k2u,galaxy.k3u,galaxy.k4u,galaxy.k5u,
galaxy.k6u,galaxy.k7u,galaxy.k8u)#,galaxy.k9a,galaxy.k10a)
#
#table.matrix <- matrix(NULL,nrow=36,ncol=10)
K <- NULL;for(i in 1:8) K <- c(K,rep(i,i))
k <- NULL;for ( i in 1:8) k <- c(k,(seq(1,i)))
meane <- NULL;for (i in 1:8) meane <- c(meane,model.list[[i]]$mass.points)
prope <- NULL;for (i in 1:8) prope <- c(prope,model.list[[i]]$masses)
sde <- NULL; for (i in 1:8)  sde <- c(sde,model.list[[i]]$sdev$sdevk)
disparitye <- NULL; for (i in 1:8) disparitye <- c(disparitye,
                                                   rep(model.list[[i]]$disparity,i))
#
meanu <- NULL;for (i in 1:8) meanu <- c(meanu,modelu.list[[i]]$mass.points)
propu <- NULL;for (i in 1:8) propu <- c(propu,modelu.list[[i]]$masses)
sdu <- NULL; for (i in 1:8)  sdu <- c(sdu,modelu.list[[i]]$sdev$sdevk)
disparityu <- NULL; for (i in 1:8) disparityu <- c(disparityu,
                                                   rep(modelu.list[[i]]$disparity,i))


###################################################
### chunk number 23:  eval=FALSE
###################################################
## galaxy.k3$disparity-galaxy.k3u$disparity
## galaxy.k4$disparity-galaxy.k4u$disparity
## galaxy.k5$disparity-galaxy.k5u$disparity
## galaxy.k6$disparity-galaxy.k6u$disparity


###################################################
### chunk number 24:  eval=FALSE
###################################################
## simulate.alldist <- function(object){
## n <- dim(object$data)[[1]]
## #nz <- length(object$masses)
## massk <- object$mass.points
## cump <- c(0,cumsum(object$masses))
## sdevk <- object$sdev$sdevk
## z <- runif(n)
## compon <- cut(sort(z),cump)
## error <- rnorm(n)*sdevk[compon]
## fval <- massk[compon]
## data <- data.frame(y=fval+error)
## data
## }


###################################################
### chunk number 25: lrt2u eval=FALSE
###################################################
## lrt2u <- NULL
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k2)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=2,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, lambda=1)
## lrt2u <- c(lrt2u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 26: 
###################################################
round(sort(lrt2u),2)


###################################################
### chunk number 27: lrt3u eval=FALSE
###################################################
## lrt3u <- NULL
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k3)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=3,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, lambda=0.999)
## lrt3u <- c(lrt3u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 28:  eval=FALSE
###################################################
## round(sort(lrt3u),2)


###################################################
### chunk number 29: lrt4u eval=FALSE
###################################################
## lrt4u <- NULL
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k4)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=4,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, lambda=0.999)
## lrt4u <- c(lrt4u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 30: 
###################################################
round(sort(lrt4u),2)


###################################################
### chunk number 31: lrt5u eval=FALSE
###################################################
## lrt5u <- NULL
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k5)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=5,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, lambda=0.999)
## lrt5u <- c(lrt5u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 32: 
###################################################
round(sort(lrt5u),2)


###################################################
### chunk number 33: lrt6u eval=FALSE
###################################################
## lrt6u <- NULL
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k6)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=6,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, lambda=0.999)
## lrt6u <- c(lrt6u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 34: 
###################################################
round(sort(lrt6u),2)


###################################################
### chunk number 35: galaxy.k1.vs.k3 eval=FALSE
###################################################
## lrts13 <- NULL
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k1)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=1,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, k=3)
## lrts13 <- c(lrts13, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 36: 
###################################################
round(sort(lrts13),4)


###################################################
### chunk number 37: simulate.alldist eval=FALSE
###################################################
## simulate.alldist <- function(object){
## n <- dim(object$data)[[1]]
## massk <- object$mass.points
## cump <- c(0,cumsum(object$masses))
## sdevk <- object$sdev$sdevk
## z <- runif(n)
## compon <- cut(sort(z),cump)
## error <- rnorm(n)*sdevk[compon]
## fval <- massk[compon]
## data <- data.frame(y=fval+error)
## data
## }


###################################################
### chunk number 38: lrt12 eval=FALSE
###################################################
## lrt12 <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k1)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=1,plot.opt=0,verbose=FALSE)
## gal.k2 <- update(gal.k1, k=2)
## lrt12 <- c(lrt12, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 39: 
###################################################
round(sort(lrt12),2)


###################################################
### chunk number 40: lrt23 eval=FALSE
###################################################
## lrt23 <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k2)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=2,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, k=3)
## lrt23 <- c(lrt23, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 41:  eval=FALSE
###################################################
## round(sort(lrt23),2)


###################################################
### chunk number 42:  eval=FALSE
###################################################
## lrt34 <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k3)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=3,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, k=4)
## lrt34 <- c(lrt34, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 43:  eval=FALSE
###################################################
## round(sort(lrt34),2)


###################################################
### chunk number 44: lrt45 eval=FALSE
###################################################
## lrt45 <- NULL
## #library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k4)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=4,plot.opt=0,verbose=FALSE,spike.protect=1)
## gal.k2 <- update(gal.k1, k=5)
## lrt45 <- c(lrt45, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 45:  eval=FALSE
###################################################
## round(sort(lrt45),2)


###################################################
### chunk number 46: lrt56 eval=FALSE
###################################################
## #library(npmlreg)
## lrt56 <- NULL
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k5)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=5,plot.opt=0,verbose=FALSE,spike.protect=1)
## gal.k2 <- update(gal.k1, k=6)
## lrt56 <- c(lrt56, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 47:  eval=FALSE
###################################################
## round(sort(lrt56),2)


###################################################
### chunk number 48: lrt67 eval=FALSE
###################################################
## lrt67 <- NULL
## #library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k6)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=6,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, k=7)
## lrt67 <- c(lrt67, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 49:  eval=FALSE
###################################################
## round(sort(lrt67),2)


###################################################
### chunk number 50: lrt36 eval=FALSE
###################################################
## lrt36 <- NULL
## #library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k3)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=3,plot.opt=0,verbose=FALSE,
## spike.protect=1)
## gal.k2 <- update(gal.k1, k=6,tol=1)
## lrt36 <- c(lrt36, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 51:  eval=FALSE
###################################################
## round(sort(lrt36),2)


###################################################
### chunk number 52:  eval=FALSE
###################################################
## round((100-sum((galaxy.k3$disparity-galaxy.k4$disparity)>round(sort(lrt34),2)))/100,2)
## round((100-sum((galaxy.k4$disparity-galaxy.k5$disparity)>round(sort(lrt45),2)))/100,2)
## round((100-sum((galaxy.k5$disparity-galaxy.k6$disparity)>round(sort(lrt56),2)))/100,2)
## round((100-sum((galaxy.k6$disparity-galaxy.k7$disparity)>round(sort(lrt67),2)))/100,2)


###################################################
### chunk number 53: lrt12u eval=FALSE
###################################################
## lrt12u <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k1)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=1,plot.opt=0,verbose=FALSE)
## gal.k2 <- update(gal.k1, k=2, lambda=1,spike.protect=1)
## lrt12u <- c(lrt12u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 54:  eval=FALSE
###################################################
## round(sort(lrt12u),2)


###################################################
### chunk number 55: lrt23u eval=FALSE
###################################################
## lrt23u <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k2u)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=2,plot.opt=0,verbose=FALSE,
## spike.protect=1,lambda=1)
## gal.k2 <- update(gal.k1, k=3)
## lrt23u <- c(lrt23u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 56:  eval=FALSE
###################################################
## round(sort(lrt23u),2)


###################################################
### chunk number 57: lrt34u eval=FALSE
###################################################
## lrt34u <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k3u)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=3,plot.opt=0,verbose=FALSE,
##                   spike.protect=1,lambda=1)
## gal.k2 <- update(gal.k1, k=4,lambda=1,spike.protect=1)
## lrt34u <- c(lrt34u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 58:  eval=FALSE
###################################################
## round(sort(lrt34u),2)


###################################################
### chunk number 59: lrt45u eval=FALSE
###################################################
## lrt45u <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k4u)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=4,plot.opt=0,verbose=FALSE,
## spike.protect=1,lambda=1)
## gal.k2 <- update(gal.k1, k=5)
## lrt45u <- c(lrt45u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 60:  eval=FALSE
###################################################
## round(sort(lrt45u),2)


###################################################
### chunk number 61: lrt56u eval=FALSE
###################################################
## lrt56u <- NULL
## library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k5u)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=5,plot.opt=0,verbose=FALSE,
## spike.protect=1,lambda=1)
## gal.k2 <- update(gal.k1, k=6)
## lrt56u <- c(lrt56u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 62:  eval=FALSE
###################################################
## round(sort(lrt56u),2)


###################################################
### chunk number 63: lrt67u eval=FALSE
###################################################
## lrt67u <- NULL
## #library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k6u)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=6,plot.opt=0,verbose=FALSE,
## spike.protect=1,lambda=0.999)
## gal.k2 <- update(gal.k1, k=7)
## lrt67u <- c(lrt67u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 64:  eval=FALSE
###################################################
## round(sort(lrt67u),2)


###################################################
### chunk number 65: lrt36u eval=FALSE
###################################################
## lrt36u <- NULL
## #library(npmlreg)
## for ( i in 1:99) {
## data <- simulate.alldist(galaxy.k3u)
## gal.k1 <- alldist(y~1,data=data,random.distribution="np",k=3,plot.opt=0,verbose=FALSE,
## spike.protect=1,lambda=1)
## gal.k2 <- update(gal.k1, k=6,lambda=0)
## lrt36u <- c(lrt36u, gal.k1$disparity - gal.k2$disparity)
## }


###################################################
### chunk number 66:  eval=FALSE
###################################################
## round(sort(-lrt36u),2)


###################################################
### chunk number 67: galaxy.cdfs eval=FALSE
###################################################
## #tmp <- NPL.bands(galaxies.df$velocity)
## print(xyplot(upper+lower~x,data=NPL.bands(galaxies.df$velocity),
## xlab="velocity",ylab="probability",
## panel=function(...){
## panel.xyplot(...,pch=20)
## panel.curve(0.085*pnorm(x,9.71,0.81)+
##             0.025*pnorm(x,16.14,0.81)+
##             0.453*pnorm(x,19.93,0.81)+
##             0.357*pnorm(x,23.05,0.81)+
##             0.044*pnorm(x,26.24,0.81)+
##             0.037*pnorm(x,33.04,0.81)
##             , lwd=2)
## panel.curve(0.08536*pnorm(x,9.71,0.64)+
##             0.878*pnorm(x,21.40,2.2038)+
##             0.037*pnorm(x,33.04,1.13)
##             ,lty=2,lwd=2)
## }))


###################################################
### chunk number 68: galaxy.6cdfs eval=FALSE
###################################################
## tmp <- NPL.bands(galaxies.df$velocity)
## print(xyplot(upper+lower~x,data=NPL.bands(galaxies.df$velocity),
## xlab="velocity",ylab="probability",
## panel=function(...){
## panel.xyplot(...,pch=20)
## masses = galaxy.k6$masses
## masspoints = galaxy.k6$mass.points
## sdevs = galaxy.k6$sdev$sdevk
## panel.curve(
##             0.085*pnorm(x,9.71,0.81)+
##             0.025*pnorm(x,16.14,0.81)+
##             0.453*pnorm(x,19.93,0.81)+
##             0.357*pnorm(x,23.05,0.81)+
##             0.044*pnorm(x,26.24,0.81)+
##             0.037*pnorm(x,33.04,0.81), lwd=2,lty=2)
## panel.curve(
##             0.08536*pnorm(x,9.71,0.64)+
##             0.878*pnorm(x,21.40,2.2038)+
##             0.037*pnorm(x,33.04,1.13),
##             lty=1,lwd=2)
## }))


###################################################
### chunk number 69: 
###################################################
#<<galaxy.cdfs>>
#tmp <- NPL.bands(galaxies.df$velocity)
print(xyplot(upper+lower~x,data=NPL.bands(galaxies.df$velocity),
xlab="velocity",ylab="probability",
panel=function(...){
panel.xyplot(...,pch=20)
panel.curve(0.085*pnorm(x,9.71,0.81)+
            0.025*pnorm(x,16.14,0.81)+
            0.453*pnorm(x,19.93,0.81)+
            0.357*pnorm(x,23.05,0.81)+
            0.044*pnorm(x,26.24,0.81)+
            0.037*pnorm(x,33.04,0.81)
            , lwd=2)
panel.curve(0.08536*pnorm(x,9.71,0.64)+
            0.878*pnorm(x,21.40,2.2038)+
            0.037*pnorm(x,33.04,1.13)
            ,lty=2,lwd=2)
}))


###################################################
### chunk number 70: galaxy.mixture.density.plot eval=FALSE
###################################################
## print(xyplot(jitter(rep(0,82),0.1)~velocity,data=galaxies.df,plot.points=T,plot=FALSE,
## xlab="velocity",ylab="density",xlim=c(5,40),ylim=c(-0.01,0.25),
## panel=function(...){
## panel.xyplot(...)
## panel.curve(
##             0.085*dnorm(x,9.71,0.81)+
##             0.025*dnorm(x,16.14,0.81)+
##             0.453*dnorm(x,19.93,0.81)+
##             0.357*dnorm(x,23.05,0.81)+
##             0.044*dnorm(x,26.04,0.81)+
##             0.037*dnorm(x,33.04,0.81),
##             lwd=1.5,lty=2)
## panel.curve(
##             0.08536*dnorm(x,9.71,0.4225)+
##             0.878*dnorm(x,21.40,2.2038)+
##             0.037*dnorm(x,33.04,0.9217171),
##             lty=1,lwd=1.5)
## }))


###################################################
### chunk number 71: 
###################################################
print(xyplot(jitter(rep(0,82),0.1)~velocity,data=galaxies.df,plot.points=T,plot=FALSE,
xlab="velocity",ylab="density",xlim=c(5,40),ylim=c(-0.01,0.25),
panel=function(...){
panel.xyplot(...)
panel.curve(
            0.085*dnorm(x,9.71,0.81)+
            0.025*dnorm(x,16.14,0.81)+
            0.453*dnorm(x,19.93,0.81)+
            0.357*dnorm(x,23.05,0.81)+
            0.044*dnorm(x,26.04,0.81)+
            0.037*dnorm(x,33.04,0.81),
            lwd=1.5,lty=2)
panel.curve(
            0.08536*dnorm(x,9.71,0.4225)+
            0.878*dnorm(x,21.40,2.2038)+
            0.037*dnorm(x,33.04,0.9217171),
            lty=1,lwd=1.5)
}))


###################################################
### chunk number 72: galaxy.density eval=FALSE
###################################################
## library(lattice)
## print(densityplot(~velocity,data=galaxies.df,
##                   xlab="velocity",ylim=c(-0.01,0.23),
## panel=function(x,...){
## panel.densityplot(x,darg=list(bw=0.5))
## panel.densityplot(x,darg=list(bw=1),plot.points=FALSE,lty=2)
## panel.densityplot(x,darg=list(bw=2),plot.points=FALSE,lty=3)
## panel.densityplot(x,darg=list(bw=4),plot.points=FALSE,lty=4)
## panel.text(c(22,20,20,20),c(0.19,0.16,0.12,0.08),c(0.5,1,2,4))
## }))


###################################################
### chunk number 73: 
###################################################
library(lattice)
print(densityplot(~velocity,data=galaxies.df,
                  xlab="velocity",ylim=c(-0.01,0.23),
panel=function(x,...){
panel.densityplot(x,darg=list(bw=0.5))
panel.densityplot(x,darg=list(bw=1),plot.points=FALSE,lty=2)
panel.densityplot(x,darg=list(bw=2),plot.points=FALSE,lty=3)
panel.densityplot(x,darg=list(bw=4),plot.points=FALSE,lty=4)
panel.text(c(22,20,20,20),c(0.19,0.16,0.12,0.08),c(0.5,1,2,4))
}))