#
#
# MLE AIC 

x<-softdrinks

n<-length(x[,1])

tab <- matrix( nrow=4, ncol=4 ) 

colnames(tab) <- c('deviance', 'dm', 'AIC', 'BIC')

# model 1
m<-lm(Time~1, data=x)
dm<- length(m$coef)+1

dev <- -2*logLik(m)[1]

AIC <- dev + dm*2
BIC <- dev + dm*log(n)

tab[1,] <- c( dev, dm, AIC, BIC )


# model 2
m<-lm(Time~Cases, data=x)
sigma <- summary(m)$sigma
dm<- length(m$coef)+1

dev <- -2*logLik(m)[1]

AIC <- dev + dm*2
BIC <- dev + dm*log(n)

tab[2,] <- c( dev, dm, AIC, BIC )

# model 3
m<-lm(Time~Distance, data=x)
dm<- length(m$coef)+1

dev <- -2*logLik(m)[1]

AIC <- dev + dm*2
BIC <- dev + dm*log(n)

tab[3,] <- c( dev, dm, AIC, BIC )

# model 4
m<-lm(Time~Cases+Distance, data=x)
dm<- length(m$coef)+1

dev <- -2*logLik(m)[1]

AIC <- dev + dm*2
BIC <- dev + dm*log(n)

tab[4,] <- c( dev, dm, AIC, BIC )

# --------------------------------------------------
# --------------------------------------------------
# --------------------------------------------------


# MINIMUM DEVIANCE 

tab2 <- matrix( nrow=4, ncol=6 ) 
colnames(tab2) <- c('minD5', 'minD20', 'minD100', 'Datmean', 'Datmed', 'Dpbf')


m.dev<-read.table('m1_dev.txt')

dev<-m.dev[,2]

tab2[1,1:3]<- c( min(dev[1:5000]), min(dev[1:20000]), min(dev[1:100000]))




# --------------------------------------------------
m.dev<-read.table('m2_dev.txt')

dev <- m.dev[,2]

tab2[2,1:3]<- c( min(dev[1:5000]), min(dev[1:20000]), min(dev[1:100000]))
# --------------------------------------------------
m.dev<-read.table('m3_dev.txt')

dev<-m.dev[,2]

tab2[3,1:3]<- c( min(dev[1:5000]), min(dev[1:20000]), min(dev[1:100000]))
# --------------------------------------------------
m.dev<-read.table('m4_dev.txt')

dev<-m.dev[,2]

tab2[4,1:3]<- c( min(dev[1:5000]), min(dev[1:20000]), min(dev[1:100000]))


# --------------------------------------------------
# --------------------------------------------------
# DEVIANCE AT MEAN

X<-as.matrix(cbind(1, x[,2:3]))

meantab <-matrix( nrow=4, ncol=4 ) 

meantab[1,]<-c( 22.45, 0,0, 0.004307)
meantab[2,]<-c( 	3.394, 2.169,0, 0.06188)
meantab[3,]<-c( 	4.937, 0,0.04267, 0.02105)
meantab[4,]<-c( 	2.357, 1.598,0.01475, 0.1066)

# ----------------------------------------------------

k<-1 
mu <- X %*% meantab[k,1:3] 
s  <- 1/sqrt(meantab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,4] <- dev
# ----------------------------------------------------

k<-2 
mu <- X %*% meantab[k,1:3] 
s  <- 1/sqrt(meantab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,4] <- dev
# ----------------------------------------------------
k<-3 
mu <- X %*% meantab[k,1:3] 
s  <- 1/sqrt(meantab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,4] <- dev
# ----------------------------------------------------
k<-4
mu <- X %*% meantab[k,1:3] 
s  <- 1/sqrt(meantab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,4] <- dev
# ----------------------------------------------------
# ----------------------------------------------------
# DEVIANCE AT MEDIAN
medtab <-matrix( nrow=4, ncol=4 ) 
medtab[1,]<-c( 22.4, 0,0, 0.004207)
medtab[2,]<-c( 	3.39, 2.169,0, 0.06049)
medtab[3,]<-c( 	4.921, 0,0.04265, 0.02055)
medtab[4,]<-c( 	2.352, 1.598,0.01474, 0.1039)


# ----------------------------------------------------

k<-1 
mu <- X %*% medtab[k,1:3] 
s  <- 1/sqrt(medtab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,5] <- dev
# ----------------------------------------------------

k<-2 
mu <- X %*% medtab[k,1:3] 
s  <- 1/sqrt(medtab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,5] <- dev
# ----------------------------------------------------
k<-3 
mu <- X %*% medtab[k,1:3] 
s  <- 1/sqrt(medtab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,5] <- dev
# ----------------------------------------------------
k<-4
mu <- X %*% medtab[k,1:3] 
s  <- 1/sqrt(medtab[k,4] )

dev <- -2*sum( dnorm( x$Time, mu, s, log=TRUE ) )

tab2[k,5] <- dev




# ----------------------------------------------------
# ----------------------------------------------------
# DEVIANCE BASED ON PBF
logPBF <- log( c(4.924E-46, 1.135E-31, 1.572E-37, 7.336E-29) )


tab2[1,6]<- -2*logPBF[1] - 2*log(2)
tab2[2,6]<- -2*logPBF[2] - 3*log(2)
tab2[3,6]<- -2*logPBF[3] - 3*log(2)
tab2[4,6]<- -2*logPBF[4] - 4*log(2)



# ----------------------------------------------------
# ----------------------------------------------------
# ----------------------------------------------------

m1.dev<-read.table('m1_dev.txt')
m2.dev<-read.table('m2_dev.txt')
m3.dev<-read.table('m3_dev.txt')
m4.dev<-read.table('m4_dev.txt')

tab3 <-matrix( nrow=4, ncol=4 ) 

colnames(tab3) <- c('min', 'p025', 'mean', 'p975')

dev<-m1.dev[1:5000,2]
tab3[1,] <- c( min(dev), quantile(dev, probs=0.025), mean(dev), quantile(dev, probs=0.975) )
dev<-m2.dev[1:5000,2]
tab3[2,] <- c( min(dev), quantile(dev, probs=0.025), mean(dev), quantile(dev, probs=0.975) )
dev<-m3.dev[1:5000,2]
tab3[3,] <- c( min(dev), quantile(dev, probs=0.025), mean(dev), quantile(dev, probs=0.975) )
dev<-m4.dev[1:5000,2]
tab3[4,] <- c( min(dev), quantile(dev, probs=0.025), mean(dev), quantile(dev, probs=0.975) )

dm<-matrix( c(2,3,3,4), nrow=4, ncol=4)
AIC<-tab3 + dm*2
BIC<-tab3 + dm*log(25)


postscript('c:/chap11_ex05_aic.ps', width = 10.0, height = 6.0, horizontal=FALSE)
#
# load package Hmisc
errbar( x=paste('model',1:4), y= AIC[,3], yminus= AIC[,2], yplus= AIC[,4],  ylab='AIC' , xaxs='r', yaxs='r', bty='l' , cex=1.2, ylim=c(120,240))


graphics.off()

postscript('c:/chap11_ex05_bic.ps', width = 10.0, height = 6.0, horizontal=FALSE)
#
# load package Hmisc
errbar( x=paste('model',1:4), y= BIC[,3], yminus= BIC[,2], yplus= BIC[,4],  ylab='BIC',xaxs='r', yaxs='r', bty='l' , cex=1.2, ylim=c(120,240))

graphics.off()




dev <- round(cbind( tab[,2],tab[,1],tab2),3)
dev2<- round(tab3,2)

AIC2<- dev[,-c(1,4,5)] +2*matrix( c(2,3,3,4), nrow=4, ncol=5 )
BIC2<- dev[,-c(1,4,5)] +  matrix( c(2,3,3,4), nrow=4, ncol=5 )*log(25)


round(AIC2,2)
round(BIC2,2)


dm<-c(2,3,3,4)
n<-length(x[,1])

postscript('c:/chap11_ex05_aic.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par(xaxs='r', yaxs='r', bty='l' , cex=0.8, cex.axis=1.5)

 boxplot(m1.dev[1:5000,2]+dm[1]*2,m2.dev[1:5000,2]+dm[2]*2,m3.dev[1:5000,2]+dm[3]*2,m4.dev[1:5000,2]+dm[4]*2,horizontal=F, names=paste('model',1:4))

graphics.off()



postscript('c:/chap11_ex05_bic.ps', width = 10.0, height = 6.0, horizontal=FALSE)
par(xaxs='r', yaxs='r', bty='l' , cex=0.8, cex.axis=1.5)
#

 boxplot(m1.dev[1:5000,2]+dm[1]*log(n),m2.dev[1:5000,2]+dm[2]*log(n),m3.dev[1:5000,2]+dm[3]*log(n),m4.dev[1:5000,2]+dm[4]*log(n),horizontal=F, names=paste('model',1:4))

graphics.off()



#
#
# calculation of AIC model weights

paic<-exp(-0.5*AIC)
temp<-apply(paic,2,sum)
paic/matrix(temp,nrow=4,ncol=4,byrow=TRUE)
round(paic/matrix(temp,nrow=4,ncol=4,byrow=TRUE),4)


paic2<-exp(-0.5*AIC2)
temp<-apply(paic2,2,sum)
paic2/matrix(temp,nrow=4,ncol=5,byrow=TRUE)
round(paic2/matrix(temp,nrow=4,ncol=5,byrow=TRUE),4)
#
#
#
# calculation of AIC model weights
pbic<-exp(-0.5*BIC)
temp<-apply(pbic,2,sum)
pbic/matrix(temp,nrow=4,ncol=4,byrow=TRUE)
round(pbic/matrix(temp,nrow=4,ncol=4,byrow=TRUE),4)


pbic2<-exp(-0.5*BIC2)
temp<-apply(pbic2,2,sum)
pbic2/matrix(temp,nrow=4,ncol=5,byrow=TRUE)
round(pbic2/matrix(temp,nrow=4,ncol=5,byrow=TRUE),4)