# dic calculation 
x<-gss1990.original


DIC<-1:3*0
AIC<-1:3*0

dm<-c( 4, 6, 6 )
pm<-1:3*0

# zip 


p0    <- c( 0.2716, 0.4133)
lambda<- c( 8.045, 	7.324)

mean.deviance <- 3663.0	

# -----------------------------------------------
# probability function of Zero Inflated Poisson
# -----------------------------------------------

dzip<-function( x, lambda, p0, log=FALSE ){
     y<-q
     ld <- dpois( x, lambda, log=TRUE )
     
     zeros <- x==0
     log.density <- 1:length(x)*0
     log.density[zeros]  <- log( p0 + (1-p0)*exp(ld[zeros]) )
     log.density[!zeros] <- ld[!zeros] + log(1-p0) 
     if (!log){ res <- exp(log.density) } else {res<-log.density}  
     return(res)
}

ly1 <- dzip( x[,1], lambda[1], p0[1], log = TRUE)
ly2 <- dzip( x[,1], lambda[2], p0[2], log = TRUE)
ll1 <- sum( x[,2]*ly1 )
ll2 <- sum( x[,3]*ly2 )

ll<- -2*(ll1+ll2)

DIC[1] <- 2*mean.deviance - ll
AIC[1] <- ll + 2*dm[1]
pm[1] <- mean.deviance - ll


# zinb

p0    <- c( 0.2239, 0.3542)
r <- c( 1.655, 1.425)
lambda<- c( 7.578, 6.693)

mean.deviance <- 2815.0	

# -----------------------------------------------
# probability function of zero inflated binomial
# -----------------------------------------------
dzinb<-function( x, lambda, size, p0, log=FALSE ){
     y<-q
     ld <- dnbinom( x, mu=lambda, size=size, log=TRUE )
     
     zeros <- x==0
     log.density <- 1:length(x)*0
     log.density[zeros]  <- log( p0 + (1-p0)*exp(ld[zeros]) )
     log.density[!zeros] <- ld[!zeros] + log(1-p0) 

     if (!log){ res <- exp(log.density) } else {res<-log.density}  
     return(res)
}



ly1 <- dzinb( x[,1], lambda[1], r[1], p0[1], log = TRUE)
ly2 <- dzinb( x[,1], lambda[2], r[2], p0[2], log = TRUE)
ll1 <- sum( x[,2]*ly1 )
ll2 <- sum( x[,3]*ly2 )

ll<- -2*(ll1+ll2)

DIC[2] <- 2*mean.deviance - ll
AIC[2] <- ll + 2*dm[2]
pm[2] <- mean.deviance - ll


# zigp

p0    <- c(0.2423, 0.3757)
omega <- c(0.5807, 0.5888)
lambda<-c(7.728, 6.928)

mean.deviance <- 2809.0	

# -----------------------------------------------
# probability function of generalized Poisson (GP)
# -----------------------------------------------
dgenpois<-function( x, lambda, omega, log=FALSE ){
     y<-x
     lambda.star <- (1-omega)*lambda + omega*y
     log.density <- log(1-omega) + log(lambda) + (y-1)*log(lambda.star) - lgamma(y+1) - lambda.star 
     if (!log){ res <- exp(log.density) } else {res<-log.density}  
     return(res)
}


# -----------------------------------------------
# probability function of zero inflated GP 
# -----------------------------------------------
dzigp<-function( x, lambda, omega, p0, log=FALSE ){

     ld <- dgenpois( x, lambda, omega, log=TRUE )
     
     zeros <- x==0
     log.density <- 1:length(x)*0
     log.density[zeros]  <- log( p0 + (1-p0)*exp(ld[zeros]) )
     log.density[!zeros] <- ld[!zeros] + log(1-p0) 
     if (!log){ res <- exp(log.density) } else {res<-log.density}  
     return(res)
}

ly1 <- dzigp( x[,1], lambda[1], omega[1], p0[1], log = TRUE)
ly2 <- dzigp( x[,1], lambda[2], omega[2], p0[2], log = TRUE)
ll1 <- sum( x[,2]*ly1 )
ll2 <- sum( x[,3]*ly2 )

ll<- -2*(ll1+ll2)

DIC[3] <- 2*mean.deviance - ll
AIC[3] <- ll + 2*dm[3]
pm[3] <- mean.deviance - ll

pm