R/ZIPInt.R

In ZIPBayes: Bayesian Methods in the Analysis of Zero-Inflated Poisson Model

ZIPInt <- function(Ystar, Covarmainphi, Covarmainmu, Covarplus, Covarminus,
                                Ystarval, Yval, Covarvalmainphi, Covarvalmainmu, Covarvalplus, Covarvalminus,
                                betaphi, betamu, alphaplus, alphaminus,
                                Uibound = c(7,11),
                                priorgamma, priormu, priorSigma, 
                                propsigmaphi,  propsigmamu = propsigmaphi, propsigmaplus = propsigmaphi,  propsigmaminus = propsigmaphi,
                                seed = 1, nmcmc = 500){
  
  set.seed(seed)
  nmain <- length(Ystar)
  nval <- length(Ystarval)
  
  nsample <- nmain + nval
  
  Zplus <- rep(0,nsample)
  Zminus<- rep(0,nsample)
  
  Covarmain1all <- rbind(Covarmainphi, Covarvalmainphi)
  Covarmain2all <- rbind(Covarmainmu, Covarvalmainmu)
  
  mu1 <-  as.matrix(rbind(Covarmain1all)) %*% t(t(betaphi))
  mu2 <- as.matrix(rbind(Covarmain2all)) %*% t(t(betamu))
  
  if (is.vector(Covarplus)) {Covarplusall <- c(Covarplus, Covarvalplus)} else {
    names(Covarvalplus) <- names(Covarplus) 
    Covarplusall <- rbind(Covarplus, Covarvalplus)
  }
  
  if (is.vector(Covarminus)) {Covarminusall <- c(Covarminus, Covarvalminus)} else {
    names(Covarvalminus) <- names(Covarminus) 
    Covarminusall <- rbind(Covarminus, Covarvalminus)
  }
  
  muZplus <- as.matrix(Covarplusall) %*% t(t(alphaplus)) 
  muZminus <- as.matrix(Covarminusall) %*% t(t(alphaminus))
  U1 <- rep(0,nsample)
  U2 <- rep(0,nsample)
  
  Y <- Ystar - Zplus[1:nmain] + Zminus[1:nmain]
  Y <- c(Y,Yval)
  
  betaphirecord <- betaphi
  betamurecord <- betamu
  alphaplusrecord <- alphaplus
  alphaminusrecord <- alphaminus
  for (repk in 1:nmcmc){
    
    ## Adjust the proposal standard error 
    if (repk %% 500 ==0){
      cat(repk," ")
    }
    
    ### Step 1: Measurement Error: Positive Process
    
    for (t in 1:nmain){
      NewGen <-  ZI_GenerateBigJoint(Ystar[t], Uibound[1], Uibound[2],  mu1[t], mu2[t], muZplus[t],  muZminus[t])
      Zplus[t] <- NewGen[1]
      Zminus[t] <- NewGen[2]
      Y[t] <- NewGen[3]
      U1[t] <- NewGen[4]
      U2[t] <- NewGen[5]
    }
    
    for (t in 1:nval){
      NewGen2 <- ZI_GenerateZpZmJoint(Ystarval[t], Yval[t], muZplus[t+nmain], muZminus[t+nmain])
      Zminus[t+nmain] <- NewGen2[1]
      Zplus[t+nmain] <- NewGen2[2]
      
      U1[t+nmain] <- GenerateU1i(Yval[t], U2[t+nmain],  mu1[t+nmain])
      U2[t+nmain] <- GenerateU2i(Yval[t], U1[t+nmain],  mu2[t+nmain])
    }
    
    alphaplussafe <- alphaplus
    
    alphaplus <- ZI_GenerateBetaMetro(alphaplus, Zplus, as.matrix(Covarplusall), propsigmaplus,  priorgamma)
    
    muZplus <- as.matrix(Covarplusall) %*% t(t(alphaplus))
    
    
    # ### Step 2: Measurement Error: Negative Process
    
    alphaminussafe <- alphaminus
    
    alphaminus <- ZI_GenerateAlphaMNMetro(alphaminus, Y, Zminus, as.matrix(Covarminusall),  propsigma= propsigmaminus,  priormu,  priorSigmas=priorSigma)
    
    muZminus <- as.matrix(Covarminusall) %*% t(t(alphaminus))
    
    # ### Step 3: Main Model
    
    ### Step 3.1 Update betaphi and beta 2
    
    betaphi <- ZI_GenerateBetaMetro(Par=betaphi, U=U1, Covar=as.matrix(Covarmain1all), propsigma=propsigmamu, priorgamma)
    mu1 <-  as.matrix(Covarmain1all) %*% t(t(betaphi))
    
    betamu <- ZI_GenerateBetaMetro(Par=betamu, U=U2, Covar=as.matrix(Covarmain2all), propsigma=propsigmaphi, priorgamma)
    
    mu2 <- as.matrix(Covarmain2all) %*% t(t(betamu))
    
    
    betaphirecord <- rbind(betaphirecord,betaphi)
    betamurecord <- rbind(betamurecord,betamu)
    if (is.vector(Covarplus)) {alphaplusrecord <- c(alphaplusrecord,alphaplus)} else {alphaplusrecord <- rbind(alphaplusrecord,alphaplus)}
    if (is.vector(Covarminus)) {alphaminusrecord <- c(alphaminusrecord,alphaminus)} else {alphaminusrecord <- rbind(alphaminusrecord,alphaminus)}
  }
  
  BayesResults <- list(betaphi_trace = betaphirecord, betamu_trace = betamurecord,  
                       alphaplus_trace = alphaplusrecord, alphaminus_trace = alphaminusrecord)
  
  class(BayesResults) <- "ZIPBayes"
  
  
  return(BayesResults)
}