ltacov.mcmc.multi <-
  function(str, ci, si, cov.gamma, cov.delta, cov.tau,
           nrc, nrs, nc, ns, niter, burn.in, t.df,
           scale, scale.delta, scale.tau,
           is.big.rho.equal.across.time, is.big.rho.equal.across.cls,
           tau.const, big.rho.const,
           rho, big.rho, beta, delta.beta, tau.beta,
           prior.rho, prior.big.rho,
           prior.beta, prior.delta.beta, prior.tau.beta,
           cov.gam.nam, cov.del.nam, cov.tau.nam,
           cov.beta, cov.delta.beta, cov.tau.beta,
           cls.nam, sta.nam, ci.nam, si.nam, time.nam, cov.nam, str.nam,
           update.cov.tau.beta, update.every, assign.sub, n.fixed)
  {
    is.assign <- T
    if( missing(assign.sub) ) { is.assign  <- F }
    if( missing(str) ) {
      str <- rep(1, dim(si)[1])
    }
    if( missing(ci) ) {
      is.cls <- F
      nci <- 0
      nc <- 1
      is.big.rho.equal.across.cls <- F
    }
    else {
      is.cls <- T
      nci <- ncol(ci)
      ncol.gam <- dim(cov.gamma)[2]
    }
    nkase <- dim(si)[1]
    nsi <- dim(si)[2]
    nt <- dim(si)[3]
    ncol.del <- dim(cov.delta)[2]
    ncol.tau <- dim(cov.tau)[2]
    npatt <- ns**nt
    ngrp <- length(table(str))
    if( update.cov.tau.beta == T ){
      update.cov <- seq(update.every, niter, update.every)
    }
    ## pattern locator
    patt.loc <- array(0, c(ns, ns, nt))
    patt.loc[,,1] <- matrix(c(1:npatt), ns, ns)
    patt.loc[,,2] <- matrix(c(1:npatt), ns, ns, byrow=T)
    ## status patterns over time
    patt.s <- matrix(0, npatt, nt)
    patt.s[1,] <- rep(1, nt)
    for(patt in 2:npatt){
      patt.s[patt,] <- patt.s[patt-1,]
      for(t in 1:nt){
        if( patt.s[patt,t] < ns) {patt.s[patt,t] <- patt.s[patt,t] + 1
                                  break}
        else {patt.s[patt,t] <- 1}
      }
    }
    ## covariates patterns
    if( ncol.tau > 1 ){
      npatt.cov <- 1
      for(col in 2:(ncol.tau)){
        npatt.cov <- npatt.cov * length(table(cov.tau[,col,1]))
      }
      patt.cov <- matrix(0, npatt.cov, (ncol.tau-1))
      for(col in 1:(ncol.tau-1)){
        patt.cov[,col] <- min(cov.tau[,(col+1),1])
      }
      for(patt in 2:npatt.cov){
        patt.cov[patt,] <- patt.cov[patt-1,]
        for(col in 1:(ncol.tau-1)){
          if( patt.cov[patt,col] < max(cov.tau[,(col+1),1]) ) {
            tmp <- names(table(cov.tau[,col+1,1])) == patt.cov[patt,col]
            tmp.1 <- (c(1:length(table(cov.tau[,col+1,1])))[tmp]) + 1
            patt.cov[patt,col] <- names(table(cov.tau[,col+1,1]))[tmp.1]
            break
          }
          else {patt.cov[patt,col] <- min(cov.tau[,(col+1),1])}
        }
      }
    }
    ## Arrays for parameters generated from MCMC
    nkase.grp <- rep(0, ngrp)
    for(grp in 1:ngrp){
      nkase.grp[grp] <- table(str)[grp]
    }
    if( is.cls == T ){
      rho.sam <- array( 0, c(niter, nci, max(nrc), nc) )
      beta.sam <- array( 0, c(niter, ncol.gam ,nc, ngrp) )
      gamma.sam <- array( 0, c(niter, nc, ngrp) )
      gamma <- array(0, c(nkase, nc))
      gam.num <- array(0, c(nkase, nc))
    }
    big.rho.sam <- array( 0, c(niter, nsi, max(nrs), ns, nc, nt) )
    delta.beta.sam <- array( 0, c(niter, ncol.del, ns, nc, ngrp) )
    tau.beta.sam <- array( 0, c(niter, ncol.tau, ns, ns, nc, nt-1, ngrp) )
    delta.sam <- array( 0, c(niter, ns, nt, nc, ngrp) )
    delta <- array( 0, c(nkase, ns, nt, nc) )
    delta.num <- array( 0, c(nkase, ns) )
    tau.sam <- array( 0, c(niter, ns, ns, nc, nt-1, ngrp) )
    tau <- array(0, c(nkase, ns, ns, nc, nt-1))
    tau.num <- array( 0, c(nkase, ns) )
    ## Arrays for confidence interval
    if( is.cls == T ){
      rho.LCI <- rho
      rho.UCI <- rho
      dimnames(rho.LCI) <- list(ci.nam, NULL, cls.nam)
      dimnames(rho.UCI) <- list(ci.nam, NULL, cls.nam)
      beta.LCI <- beta
      beta.UCI <- beta
      dimnames(beta.LCI) <- list(cov.gam.nam, cls.nam, str.nam)
      dimnames(beta.UCI) <- list(cov.gam.nam, cls.nam, str.nam)
      gamma.LCI <- matrix( 0, nc, ngrp )
      gamma.UCI <- matrix( 0, nc, ngrp )
      dimnames(gamma.LCI) <- list(cls.nam, str.nam)
      dimnames(gamma.UCI) <- list(cls.nam, str.nam)
    }
    big.rho.LCI <- big.rho
    big.rho.UCI <- big.rho
    dimnames(big.rho.LCI) <- list(si.nam, NULL, sta.nam, cls.nam, time.nam)
    dimnames(big.rho.UCI) <- list(si.nam, NULL, sta.nam, cls.nam, time.nam)
    delta.beta.LCI <- delta.beta
    delta.beta.UCI <- delta.beta
    dimnames(delta.beta.LCI) <- list(cov.del.nam, sta.nam, cls.nam, str.nam)
    dimnames(delta.beta.UCI) <- list(cov.del.nam, sta.nam, cls.nam, str.nam)
    tau.beta.LCI <- tau.beta
    tau.beta.UCI <- tau.beta
    dimnames(tau.beta.LCI) <-
      list(cov.tau.nam, sta.nam, sta.nam, cls.nam, time.nam[2:nt], str.nam)
    dimnames(tau.beta.UCI) <- 
      list(cov.tau.nam, sta.nam, sta.nam, cls.nam, time.nam[2:nt], str.nam)
    delta.LCI <- array( 0, c(ns, nt, nc, ngrp) )
    delta.UCI <- array( 0, c(ns, nt, nc, ngrp) )
    dimnames(delta.LCI) <- list(sta.nam, time.nam, cls.nam, str.nam)
    dimnames(delta.UCI) <- list(sta.nam, time.nam, cls.nam, str.nam)
    tau.LCI <- array( 0, c(ns, ns, nc, nt-1, ngrp) )    
    tau.UCI <- array( 0, c(ns, ns, nc, nt-1, ngrp) )
    dimnames(tau.LCI) <-
      list(sta.nam, sta.nam, cls.nam, time.nam[2:nt], str.nam)
    dimnames(tau.UCI) <-
      list(sta.nam, sta.nam, cls.nam, time.nam[2:nt], str.nam)
    ## Arrays needed in calculation
    comlik <- array(0, c(nkase, npatt, nc))
    post <- array(0, c(nkase, npatt, nc))
    ## Array for acceptance rate for Metropolis
    if( is.cls == T ){
      acp.beta <- matrix(0, niter, ngrp)
    }
    acp.delta.beta <- array(0, c(niter, nc, ngrp))
    acp.tau.beta <- array(0, c(niter, ns, nc, nt-1, ngrp))
    ###################
    ## starting MCMC ##
    ###################
    iter <- 0
    update.seq <- 1
    while( iter < niter ){
      iter <- iter + 1
      cat(iter, ".")
      ## update parameters
      if( is.cls == T ){
        rho.sam[iter,,,] <- rho
        beta.sam[iter,,,] <- beta
      }
      big.rho.sam[iter,,,,,] <- big.rho
      delta.beta.sam[iter,,,,] <- delta.beta
      tau.beta.sam[iter,,,,,,] <- tau.beta
      ############
      ## I-STEP ##
      ############      
      ## calculate gamma and posterior prob for the static class
      if( is.cls == T ) {
        tmp.rho <- array(1, c(nkase, nci, nc))
        for( grp in 1:ngrp ){
          w <- str == grp
          cov.tmp <- matrix(cov.gamma[w,], nkase.grp[grp], ncol.gam)
          beta.tmp <- matrix(beta[,,grp], ncol.gam, nc)
          gam.num[w,] <- exp( cov.tmp %*% beta.tmp )
          gamma[w,] <- gam.num[w,] / apply( gam.num[w,], 1, sum )          
          gamma.sam[iter,,grp] <- apply( gamma[w,], 2, mean )
          ## product of rho
          for(item in 1:nci){
            for(k in 1:nrc[item]){
              tmp <- (ci[,item] == k) & w
              for(cls in 1:nc){
                tmp.rho[tmp,item,cls] <- rho[item,k,cls]
              }
            }
          }
        }
        rho.prod <- matrix(0, nkase, nc)
        for(cls in 1:nc){
          rho.prod[,cls] <- apply( tmp.rho[,,cls], 1, prod )
        }
        num.cls <- gamma * rho.prod
      }
      if( is.cls == F ) {
        num.cls <- matrix(1, nkase, nc)
      }
      ## calculate delta at time 1
      for( grp in 1:ngrp ){
        w <- str == grp
        cov.tmp <- matrix(cov.delta[w,], nkase.grp[grp], ncol.del)
        for(cls in 1:nc){
          delta.beta.tmp <- matrix(delta.beta[,,cls,grp], ncol.del, ns)
          delta.num[w,] <- exp( cov.tmp %*% delta.beta.tmp )
          delta[w,,1,cls] <- delta.num[w,] / apply( delta.num[w,], 1, sum )
        }     
        delta.sam[iter,,1,,grp] <- apply(array(delta[w,,1,], c(sum(w),ns,nc)),
                                         c(2,3), mean)
      }
      ## calculate tau
      tmp.rho <- array(1, c(nkase, nsi, ns, nc, nt))
      for( grp in 1:ngrp ){
        w <- str == grp
        for(t in 1:(nt-1)){
          cov.tmp <- matrix(cov.tau[w,,t+1], nkase.grp[grp], ncol.tau)
          for(cls in 1:nc){
            for(prev.s in 1:ns){
              tmp <- tau.const[,prev.s,cls,t,grp] == 1
              tau.beta.tmp <-
                matrix(tau.beta[,tmp,prev.s,cls,t,grp], ncol.tau, sum(tmp))
              tau.num[w,tmp] <- exp( cov.tmp %*% tau.beta.tmp )
              tau.num.tmp <- matrix(tau.num[w,tmp], nkase.grp[grp], sum(tmp))
              tau[w,tmp,prev.s,cls,t] <-
                tau.num[w,tmp] / apply( tau.num.tmp, 1, sum )
            }
          }
        }
        tau.sam[iter,,,,,grp] <-
          apply(array(tau[w,,,,], c(sum(w),ns,ns,nc,nt-1)), c(2,3,4,5), mean)
        ## product of big rho
        for(time in 1:nt){
          for(item in 1:nsi){
            for(k in 1:nrs[item]){
              tmp <- (si[,item,time] == k) & w
              for(cls in 1:nc){
                for(sta in 1:ns){
                  tmp.rho[tmp,item,sta,cls,time] <-
                    big.rho[item,k,sta,cls,time]
                }
              }
            }
          }
        }
      }
      rho.prod <- array(0, c(nkase, ns, nc, nt))
      for(time in 1:nt){
        for(cls in 1:nc){
          for(sta in 1:ns){
            rho.prod[,sta,cls,time] <-
              apply( tmp.rho[,,sta,cls,time], 1, prod )
          }
        }
      }
      ## calculate delta at time t
      for( grp in 1:ngrp ){
        w <- str == grp
        for(cls in 1:nc){          
          for(t in 1:(nt-1)){
            for(sta in 1:ns){
              tmp <- rep(0, sum(w))
              for(pre.s in 1:ns){
                tmp <- tmp + delta[w,pre.s,t,cls]*tau[w,sta,pre.s,cls,t]
              }
              delta[w,sta,t+1,cls] <- tmp
              delta.sam[iter,sta,t+1,cls,grp] <- mean(tmp)
            }
          }
        }   
      }
      ## calculate complete data likelihood given cls
      for( grp in 1:ngrp ){
        w <- str == grp
        for(cls in 1:nc){
          tmp <- delta[w,,1,cls] * rho.prod[w,,cls,1]
          for(patt in 1:npatt){
            tmp.patt <- patt.s[patt,]
            comlik[w,patt,cls] <- num.cls[w,cls] * 
              tmp[,tmp.patt[1]] * tau[w,tmp.patt[2],tmp.patt[1],cls,1] *
                rho.prod[w,tmp.patt[2],cls,2]
          }
        }
      }
      ## calculate joint posterior probs
      post <- comlik / apply(comlik, 1, sum)
      ## Impute class and status membership
      imp <- array(0, c(nkase, npatt, nc))
      ran <- runif(nkase)
      for(kase in 1:nkase){
        cls <- 1
        patt <- 1
        repeat{
          if( ran[kase] < post[kase,patt,cls] ) break
          ran[kase] <- ran[kase] - post[kase,patt,cls]
          if( patt == npatt ) {
            cls <- cls + 1
            patt <- 0
          }
          patt <- patt + 1
        }
        imp[kase,patt,cls] <- 1
      }
      ## assign individual to status if any
      if( is.assign == T ){
        for(cls in 1:nc){
          for(patt in 1:npatt){
            w.1 <- patt.s[patt,1]
            w.2 <- patt.s[patt,2]
            if( tau.const[w.2,w.1,cls,1,1] == 1 ){
              for(ii in (nt+1):dim(assign.sub)[2]){
                for(jj in 1:(n.fixed[w.2,w.1,cls,1,1])){
                  imp[assign.sub[patt,ii,cls,jj],,] <- 0
                  imp[assign.sub[patt,ii,cls,jj],patt,cls] <- 1
                }
              }
            }
          }
        }
      }
      ## imputed class membership
      if( is.cls == T ){
        imp.cls <- matrix(0, nkase, nc)
        for(cls in 1:nc){
          imp.cls[,cls] <- apply(imp[,,cls], 1, sum)
        }
      }
      if( is.cls == F ){
        imp.cls <- matrix(1, nkase, nc)
      }
      ## imputed class and status membership at time t
      imp.cls.sta.at.t <- array(0, c(nkase, ns, nc, nt))
      for(kase in 1:nkase){
        for(cls in 1:nc){
          for(time in 1:nt){
            for(sta in 1:ns){
              imp.cls.sta.at.t[kase,sta,cls,time] <-
                 sum(imp[kase,patt.loc[sta,,time],cls])
            }
          }
        }
      }
      ## impute adjacent status membership
      imp.adj.sta <- array(0, c(nkase, ns, ns, nc, nt-1))
      for(time in 1:(nt-1)){
        for(cls in 1:nc){
          for(patt in 1:npatt){
            tmp <- patt.s[patt,]
            imp.adj.sta[,tmp[2],tmp[1],cls,time] <- imp[,patt,cls]
          }
        }
      }      
      ## Impute individual's missing item for static class
      if( is.cls == T ){
        ci.imp <- ci
        for(kase in 1:nkase){
          for(item in 1:nci){
            if(ci[kase,item] == 0){
              for(cls in 1:nc){
                if(imp.cls[kase,cls] == 1){
                  ran <- runif(1)
                  for(k in 1:nrc[item]){
                    if( ran <= rho[item,k,cls] ){
                      ci.imp[kase,item] <- k
                      break
                    }
                    else{ ran <- ran - rho[item,k,cls] }
                  }
                }
              }
            }
          }
        }
      }
      ## Impute individual's missing item for latent status
      si.imp <- si
      for(t in 1:nt){
        for(kase in 1:nkase){
          for(item in 1:nsi){
            if(si[kase,item,t] == 0){
              for(cls in 1:nc){
                for(sta in 1:ns){
                  if(imp.cls.sta.at.t[kase,sta,cls,t] == 1){
                    ran <- runif(1)
                    for(k in 1:nrs[item]){
                      if( ran <= big.rho[item,k,sta,cls,t] ){
                        si.imp[kase,item,t] <- k
                        break
                      }
                      else{ ran <- ran - big.rho[item,k,sta,cls,t] }
                    }
                  }
                }
              }
            }
          }
        }
      }   
      ############
      ## P-STEP ##
      ############
      ## generate rho parameters
      if( is.cls == T ){
        for(cls in 1:nc){
          for(item in 1:nci){
            gam <- rep(0, nrc[item])
            for(k in 1:nrc[item]){
              tmp <- sum( imp.cls[,cls] * (ci.imp[,item]==k) )
              gam[k] <- rgamma(1, shape=tmp+prior.rho[item,k,cls])
            }
            for(k in 1:nrc[item]){
              rho[item,k,cls] <- gam[k] / sum(gam)
            }
          }
        }
      }
      ## generate big.rho parameters
      for(sta in 1:ns){
        for(item in 1:nsi){
          gam <- rep(0, nrs[item])
          if( is.big.rho.equal.across.cls == F ){          
            for(cls in 1:nc){            
              if( is.big.rho.equal.across.time == T ){
                for(k in 1:nrs[item]){
                  if( big.rho.const[item,k,sta,cls,1] == 1 ){
                    tmp <-
                      sum( imp.cls.sta.at.t[,sta,cls,] * (si.imp[,item,]==k) )
                    gam[k] <-
                      rgamma(1, shape=tmp+prior.big.rho[item,k,sta,cls,1])
                  }
                  if( big.rho.const[item,k,sta,cls,1] == 0 ){
                    gam[k] <- 0
                  }
                }
                for(k in 1:nrs[item]){
                  if( big.rho.const[item,k,sta,cls,1] == 1 ){
                    big.rho[item,k,sta,cls,1:nt] <- gam[k] / sum(gam)
                  }
                }
              }
              if( is.big.rho.equal.across.time == F ){
                for(t in 1:nt){
                  for(k in 1:nrs[item]){
                    if( big.rho.const[item,k,sta,cls,t] == 1 ){
                      tmp <- sum( imp.cls.sta.at.t[,sta,cls,t] *
                                 (si.imp[,item,t]==k) )
                      gam[k] <-
                        rgamma(1, shape=tmp+prior.big.rho[item,k,sta,cls,t])
                    }
                    if( big.rho.const[item,k,sta,cls,t] == 0 ){
                      gam[k] <- 0
                    }
                  }
                  for(k in 1:nrs[item]){
                    if( big.rho.const[item,k,sta,cls,t] == 1 ){
                      big.rho[item,k,sta,cls,t] <- gam[k] / sum(gam)
                    }
                  }
                }
              }          
            }
          }
          if( is.big.rho.equal.across.cls == T ){
            if( is.big.rho.equal.across.time == F ){
              for(t in 1:nt){
                for(k in 1:nrs[item]){
                  if( big.rho.const[item,k,sta,1,t] == 1 ){
                    tmp <- sum( imp.cls.sta.at.t[,sta,,t] *
                               (si.imp[,item,t]==k) )
                    gam[k] <-
                      rgamma(1, shape=tmp+prior.big.rho[item,k,sta,1,t])
                  }
                  if( big.rho.const[item,k,sta,1,t] == 0 ){
                    gam[k] <- 0
                  }
                }
                for(k in 1:nrs[item]){
                  if( big.rho.const[item,k,sta,1,t] == 1 ){
                    big.rho[item,k,sta,1:nc,t] <- gam[k] / sum(gam)
                  }
                }
              }
            }
            if( is.big.rho.equal.across.time == T ){
              tmp <- matrix(0, nt, nrs[item])
              for(t in 1:nt){
                for(k in 1:nrs[item]){
                  tmp[t,k] <-
                    sum( imp.cls.sta.at.t[,sta,,t] * (si.imp[,item,t]==k) )
                }
              }
              tmp <- apply(tmp, 2, sum)
              for(k in 1:nrs[item]){
                if( big.rho.const[item,k,sta,1,1] == 1 ){
                  gam[k] <-
                    rgamma(1, shape=tmp[k]+prior.big.rho[item,k,sta,1,1])
                }
                if( big.rho.const[item,k,sta,1,1] == 0 ){
                  gam[k] <- 0
                }
              }
              for(k in 1:nrs[item]){
                if( big.rho.const[item,k,sta,1,1] == 1 ){
                  big.rho[item,k,sta,1:nc,1:nt] <- gam[k] / sum(gam)
                }
              }
            }
          }
        }
      }
      ## generate beta parameters
      if( is.cls == T ){
        for( grp in 1:ngrp ){
          w <- str == grp
          ###############################
          if( update.cov.tau.beta == T ){
            if(iter == update.cov[update.seq]){
              cov.beta[,,grp] <-
                cov(matrix(beta.sam[1:update.cov[update.seq],
                                    1:ncol.gam,1:(nc-1),grp],
                           update.cov[update.seq], (nc-1)*ncol.gam, byrow=F))
            }
          }
          ###############################
          if(t.df >= 20){
            tmp <- matrix(rnorm(ncol.gam*(nc-1)), ncol=ncol.gam*(nc-1)) %*%
              chol(cov.beta[,,grp])
          }
          if(t.df < 20){
            tmp <-
              matrix(rt(ncol.gam*(nc-1), df=t.df), ncol=ncol.gam*(nc-1)) %*%
              chol(cov.beta[,,grp])
          }
          tmp <- matrix(c(tmp, rep(0, ncol.gam)) , ncol.gam, nc)
          beta.tmp <- (scale[grp]/sqrt(ncol.gam*(nc-1))) * tmp +
            matrix(beta[,,grp], ncol.gam, nc)
          cov.tmp <- matrix(cov.gamma[w,], nkase.grp[grp], ncol.gam)
          gam.num[w,] <- exp(cov.tmp %*% beta.tmp)
          gamma.tmp <- gam.num[w,] / apply( gam.num[w,], 1, sum )
          ##############################################################
          ## put prior in prior
          prior <- matrix(0, nkase, nc)
          for(cls in 1:nc){
            imp.tmp <- imp.cls[w,cls]
            prior.tmp <- prior[w,cls]            
#            if( ncol.gam == 1 ){
              if( sum(imp.tmp) == 0 ){
                prior.tmp <- prior.tmp + prior.beta[cls,grp] / sum(w)
              }
#            }
#            if( ncol.gam >= 2 ){
#              for( ii in 1:npatt.cov ){
#                w.1 <- apply(t(cov.gamma[w,2:ncol.gam])==patt.cov[ii,],
#                             2, prod) == 1
#                if( sum(imp.tmp[w.1]) == 0 ){
#                  prior.tmp[w.1] <- prior.tmp[w.1] +
#                    prior.beta[cls,grp] / sum(w.1)
#                }
#              }
#            }
            prior[w,cls] <- prior.tmp
          }
          ##############################################################
          tmp <-
            apply( (gamma.tmp / gamma[w,]) ** (imp.cls[w,] + prior[w,]),
                  1, prod )
          acp.beta[iter,grp] <- min( 1, prod(tmp) )
          ran <- runif(1)
          if( ran <= acp.beta[iter,grp] ) { beta[,,grp] <- beta.tmp }
        }
      }
      ## generate delta.beta parameters
      for( grp in 1:ngrp ){
        w <- str == grp
        for(cls in 1:nc){
          ###############################
          if( update.cov.tau.beta == T ){
            if(iter == update.cov[update.seq]){
              cov.delta.beta[,,cls,grp] <-
                cov(matrix(delta.beta.sam[1:update.cov[update.seq],
                                          1:ncol.del,1:(ns-1),cls,grp],
                           update.cov[update.seq], (ns-1)*ncol.del, byrow=F))
            }
          }
          ###############################
          if(t.df >= 20){
            tmp <- matrix(rnorm(ncol.del*(ns-1)), ncol=ncol.del*(ns-1)) %*%
              chol(cov.delta.beta[,,cls,grp])
          }
          if(t.df < 20){
            tmp <- matrix(rt(ncol.del*(ns-1), df=t.df),
                          ncol=ncol.del*(ns-1)) %*%
              chol(cov.delta.beta[,,cls,grp])
          }          
          tmp <- matrix(c(tmp, rep(0, ncol.del)), ncol.del, ns)
          delta.beta.tmp <- (scale.delta[cls,grp]/sqrt(ncol.del*(ns-1))) *
            tmp + matrix(delta.beta[,,cls,grp], ncol.del, ns)
          cov.tmp <- matrix(cov.delta[w,], nkase.grp[grp], ncol.del)
          delta.num[w,] <- exp(cov.tmp %*% delta.beta.tmp)
          delta.tmp <- delta.num[w,] / apply( delta.num[w,], 1, sum )
          ############################################################
          ## put prior in prior
          prior <- matrix(0, nkase, ns)
          ww <- (imp.cls[,cls]==1) & w
          for(sta in 1:ns){
            imp.tmp <- imp.cls.sta.at.t[ww,sta,cls,1]
            prior.tmp <- prior[ww,sta]           
#            if( ncol.del == 1 ){
              if( sum(imp.tmp) == 0 ){
                prior.tmp <- prior.tmp +
                  prior.delta.beta[sta,cls,grp] / sum(ww)
              }
#            }
#            if( ncol.del >= 2 ){
#              for( ii in 1:npatt.cov ){
#                w.1 <- apply(t(cov.del[ww,2:ncol.del])==patt.cov[ii,],
#                             2, prod) == 1
#                if( sum(imp.tmp[w.1]) == 0 ){
#                  prior.tmp[w.1] <- prior.tmp[w.1] +
#                    prior.delta.beta[sta,cls,grp] / sum(w.1)
#                }
#              }
#            }
            prior[ww,sta] <- prior.tmp
          }
          ############################################################
          tmp <- apply( (delta.tmp / delta[w,,1,cls]) **
                       (imp.cls.sta.at.t[w,,cls,1] + prior[w,]), 1, prod )
          if( is.cls == T ){
            acp.delta.beta[iter,cls,grp] <- min(1, prod(tmp[ww]))
          }
          else{
            acp.delta.beta[iter,cls,grp] <- min(1, prod(tmp))
          }
          ran <- runif(1)
          if( ran <= acp.delta.beta[iter,cls,grp] ) {
            delta.beta[,,cls,grp] <- delta.beta.tmp
          }
        }
      }
      ## generate tau.beta.parameters     
      for( grp in 1:ngrp ) {
        w.1 <- str == grp
        kk <- 0
        for(cls in 1:nc){
          for(t in 1:(nt-1)){
            itmp <- apply(tau.const[,,cls,t,grp], 2, sum) - 1
            w <- itmp != 0
            w <- (1:ns)[w]
            for(pre.s in w){
              kk <- kk + 1     
              n.fr.beta <- ncol.tau*itmp[pre.s]
              
              a.tmp <- tau.const[,pre.s,cls,t,grp] == 1
              b.tmp <- (1:ns)[a.tmp]
              a.tmp[pre.s] <- FALSE
              ###############################
              if( update.cov.tau.beta == T ){
                if(iter == update.cov[update.seq]){
                  cov.tau.beta[1:n.fr.beta,1:n.fr.beta,pre.s,cls,t,grp] <-
                    cov(matrix(tau.beta.sam[1:update.cov[update.seq],
                                            1:ncol.tau,a.tmp,pre.s,cls,t,grp],
                               update.cov[update.seq], n.fr.beta, byrow=F))
                }
              }
              ###############################
              if( t.df >= 20 ){
                tmp <- matrix(rnorm(n.fr.beta), ncol=n.fr.beta) %*%
                  chol(cov.tau.beta[1:n.fr.beta,1:n.fr.beta,pre.s,cls,t,grp])
              }
              if( t.df < 20 ){
                tmp <- matrix(rt(n.fr.beta, df=t.df), ncol=n.fr.beta) %*%
                  chol(cov.tau.beta[1:n.fr.beta,1:n.fr.beta,pre.s,cls,t,grp])
              }
              tmp <- matrix(tmp, ncol.tau, itmp[pre.s])            
              tau.beta.tmp <- matrix(0, ncol.tau, ns)
              tau.beta.tmp[,a.tmp] <- tmp
              tau.beta.tmp <- (scale.tau[cls,pre.s,grp]/sqrt(n.fr.beta)) *
                tau.beta.tmp +
                  matrix(tau.beta[,,pre.s,cls,t,grp], ncol.tau, ns)
              a.tmp <- tau.const[,pre.s,cls,t,grp] == 1
              cov.tmp <- matrix(cov.tau[w.1,,t+1], nkase.grp[grp], ncol.tau)
              tau.num[w.1,a.tmp] <- exp( cov.tmp %*% tau.beta.tmp[,a.tmp])
              tau.tmp <-
                tau.num[w.1,a.tmp] / apply( tau.num[w.1,a.tmp], 1, sum )
              ##############################################################
              ## put prior in prior
              prior <- matrix( 0, nkase, sum(a.tmp) )
              ww <- (imp.cls.sta.at.t[,pre.s,cls,t] == 1) & w.1
              jj <- 0
              for(sta in b.tmp){
                jj <- jj + 1
                imp.tmp <- imp.adj.sta[ww,sta,pre.s,cls,t]
                prior.tmp <- prior[ww,jj]
                if( ncol.tau == 1 ){
                  if( sum(imp.tmp) == 0 ){
                    prior.tmp <- prior.tmp +
                      prior.tau.beta[sta,pre.s,cls,grp] / sum(ww)
                  }
                }
                if( ncol.tau >= 2 ){
                  for( ii in 1:npatt.cov ){
                    w <- apply(t(cov.tau[ww,2:ncol.tau,1])==patt.cov[ii,],
                               2, prod) == 1
                    if( sum(imp.tmp[w]) == 0 ){
                      prior.tmp[w] <- prior.tmp[w] +
                        prior.tau.beta[sta,pre.s,cls,grp] / sum(w)
                    }
                  }
                }
                prior[ww,jj] <- prior.tmp
              }
              ##############################################################
              tmp <- apply( (tau.tmp / tau[w.1,a.tmp,pre.s,cls,t]) **
                           (imp.adj.sta[w.1,a.tmp,pre.s,cls,t] +
                            prior[w.1,] ), 1, prod )
              acp.tau.beta[iter,pre.s,cls,t,grp] <- min(1, prod(tmp[ww]))
              ran <- runif(1)
              if( ran <= acp.tau.beta[iter,pre.s,cls,t,grp] ) {
                tau.beta[,,pre.s,cls,t,grp] <- tau.beta.tmp
              }
            }
          }
        }
      }
      if( update.cov.tau.beta == T ){
        if(iter == update.cov[update.seq]){
          update.seq <- update.seq + 1
        }
      }
    }
    ## calculate CI for parameters
    n <- niter - burn.in
    int.loc <- round( n*0.025 )
    if( is.cls == T ){
      for(cls in 1:nc){
        for( grp in 1:ngrp ) {
          tmp <- sort(gamma.sam[(burn.in+1):niter,cls,grp])
          gamma.LCI[cls,grp] <- tmp[int.loc]
          gamma.UCI[cls,grp] <- tmp[(n-int.loc)]
          for(col in 1:ncol.gam){
            tmp <- sort(beta.sam[(burn.in+1):niter,col,cls,grp])
            beta.LCI[col,cls,grp] <- tmp[int.loc]
            beta.UCI[col,cls,grp] <- tmp[(n-int.loc)]
          }
        }
        for(item in 1:nci){
          for(k in 1:nrc[item]){
            tmp <- sort(rho.sam[(burn.in+1):niter,item,k,cls])
            rho.LCI[item,k,cls] <- tmp[int.loc]
            rho.UCI[item,k,cls] <- tmp[(n-int.loc)]
          }
        }
      }
    }
    for(cls in 1:nc){
      for(sta in 1:ns){
        for( grp in 1:ngrp ) {
          for(t in 1:nt){
            tmp <- sort(delta.sam[(burn.in+1):niter,sta,t,cls,grp])
            delta.LCI[sta,t,cls,grp] <- tmp[int.loc]
            delta.UCI[sta,t,cls,grp] <- tmp[(n-int.loc)]
          }
          for(col in 1:ncol.del){
            tmp <- sort(delta.beta.sam[(burn.in+1):niter,col,sta,cls,grp])
            delta.beta.LCI[col,sta,cls,grp] <- tmp[int.loc]
            delta.beta.UCI[col,sta,cls,grp] <- tmp[(n-int.loc)]
          }
          for(t in 1:(nt-1)){
            for(sta.p in 1:ns){
              tmp <- sort(tau.sam[(burn.in+1):niter,sta.p,sta,cls,t,grp])
              tau.LCI[sta.p,sta,cls,t,grp] <- tmp[int.loc]
              tau.UCI[sta.p,sta,cls,t,grp] <- tmp[(n-int.loc)]
              for(col in 1:ncol.tau){
                tmp <- sort(tau.beta.sam[(burn.in+1):niter,
                                         col,sta.p,sta,cls,t,grp])
                tau.beta.LCI[col,sta.p,sta,cls,t,grp] <- tmp[int.loc]
                tau.beta.UCI[col,sta.p,sta,cls,t,grp] <- tmp[(n-int.loc)]
              }
            }
          }
        }
        for(t in 1:nt){
          for(item in 1:nsi){
            for(k in 1:nrs[item]){
              tmp <- sort(big.rho.sam[(burn.in+1):niter,item,k,sta,cls,t])
              big.rho.LCI[item,k,sta,cls,t] <- tmp[int.loc]
              big.rho.UCI[item,k,sta,cls,t] <- tmp[(n-int.loc)]
            }
          }
        }
      }
    }
    ## save param
    param <- list()   
    if( is.cls == T ){
      param$rho <- apply(rho.sam[(burn.in+1):niter,,,], c(2,3,4), mean)    
      if(ngrp >= 2){
        param$gamma <- apply(gamma.sam[(burn.in+1):niter,,], c(2,3), mean)
        param$delta <-
          apply(delta.sam[(burn.in+1):niter,,,,], c(2,3,4,5), mean)
        param$tau <- apply(tau.sam[(burn.in+1):niter,,,,,], c(2,3,4,5), mean)
      }
      if(ngrp == 1){
        param$gamma <- apply(gamma.sam[(burn.in+1):niter,,], 2, mean)
        param$delta <- apply(delta.sam[(burn.in+1):niter,,,,], c(2,3,4), mean)
        param$tau <- apply(tau.sam[(burn.in+1):niter,,,,,], c(2,3,4), mean)
      }
      if(ncol.gam >= 2){
        if(ngrp >= 2){
          param$beta <- apply(beta.sam[(burn.in+1):niter,,,], c(2,3,4), mean)
        }
        if(ngrp == 1){
          param$beta <- apply(beta.sam[(burn.in+1):niter,,,], c(2,3), mean)
        }
      }
      if(ncol.del >= 2){
        if(ngrp >= 2){
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], c(2,3,4,5), mean)
        }
        if(ngrp == 1){
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], c(2,3,4), mean)
        }
      }
      if(ncol.tau >= 2){
        if(ngrp >= 2){
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3,4,5,6), mean)
        }
        if(ngrp == 1){
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3,4,5), mean)
        }
      }
      if(ncol.gam == 1){
        if(ngrp >= 2){
          param$beta <- apply(beta.sam[(burn.in+1):niter,,,], c(2,3), mean)
        }
        if(ngrp == 1){ 
          param$beta <- apply(beta.sam[(burn.in+1):niter,,,], 2, mean)
        }
      }
      if(ncol.del == 1){
        if(ngrp >= 2){
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], c(2,3,4), mean)
        }
        if(ngrp == 1){ 
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], c(2,3), mean)
        }
      }
      if(ncol.tau == 1){
        if(ngrp >= 2){
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3,4,5), mean)
        }
        if(ngrp == 1){ 
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3,4), mean)
        }
      }
      param$big.rho <-
        apply(big.rho.sam[(burn.in+1):niter,,,,,], c(2,3,4,5,6), mean)
      dimnames(param$rho) <- list(ci.nam, NULL, cls.nam)
      if(ngrp >= 2){
        dimnames(param$gamma) <- list(cls.nam, str.nam)
        dimnames(param$delta) <- list(sta.nam, time.nam, cls.nam, str.nam)
        dimnames(param$tau) <- list(sta.nam, sta.nam, cls.nam, str.nam)
      }
      if(ngrp == 1){
        names(param$gamma) <- cls.nam
        dimnames(param$delta) <- list(sta.nam, time.nam, cls.nam)
        dimnames(param$tau) <- list(sta.nam, sta.nam, cls.nam)
      }
      if(ncol.gam >= 2){
        if(ngrp >= 2){
          dimnames(param$beta) <- list(cov.gam.nam, cls.nam, str.nam)
        }
        if(ngrp == 1){
          dimnames(param$beta) <- list(cov.gam.nam, cls.nam)
        }
      }
      if(ncol.del >= 2){
        if(ngrp >= 2){
          dimnames(param$delta.beta) <-
            list(cov.del.nam, sta.nam, cls.nam, str.nam)
        }
        if(ngrp == 1){
          dimnames(param$delta.beta) <- list(cov.del.nam, sta.nam, cls.nam)
        }
      }
      if(ncol.tau >= 2){
        if(ngrp >= 2){
          dimnames(param$tau.beta) <-
            list(cov.tau.nam, sta.nam, sta.nam, cls.nam, str.nam)
        }
        if(ngrp == 1){
          dimnames(param$tau.beta) <-
            list(cov.tau.nam, sta.nam, sta.nam, cls.nam)
        }
      }
      if(ncol.gam == 1){
        if(ngrp >= 2){
          dimnames(param$beta) <- list(cls.nam, str.nam)
        }
        if(ngrp == 1){
          names(param$beta) <- cls.nam
        }
      }
      if(ncol.del == 1){
        if(ngrp >= 2){
          dimnames(param$delta.beta) <- list(sta.nam, cls.nam, str.nam)
        }
        if(ngrp == 1){
          dimnames(param$delta.beta) <- list(sta.nam, cls.nam)
        }
      }
      if(ncol.tau == 1){
        if(ngrp >= 2){
          dimnames(param$tau.beta) <- list(sta.nam, sta.nam, cls.nam, str.nam)
        }
        if(ngrp == 1){
          dimnames(param$tau.beta) <- list(sta.nam, sta.nam, cls.nam)
        }
      }
      dimnames(param$big.rho) <- list(si.nam, NULL, sta.nam, cls.nam, time.nam)
    }
    if( is.cls == F ){
      if(ngrp >= 2){
        param$delta <- apply(delta.sam[(burn.in+1):niter,,,,], c(2,3,5), mean)
        param$tau <- apply(tau.sam[(burn.in+1):niter,,,,,], c(2,3,4), mean)
      }
      if(ngrp == 1){
        param$delta <- apply(delta.sam[(burn.in+1):niter,,,,], c(2,3), mean)
        param$tau <- apply(tau.sam[(burn.in+1):niter,,,,,], c(2,3), mean)
      }
      
      if(ncol.del >= 2){
        if(ngrp >= 2){
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], c(2,3,4), mean)
        }
        if(ngrp == 1){
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], c(2,3), mean)
        }
      }      
      if(ncol.tau >= 2){
        if(ngrp >= 2){
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3,4,5), mean)
        }
        if(ngrp == 1){
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3,4), mean)
        }
      }
      
      param$delta.beta <- array(param$delta.beta, c(ncol.del, ns, nc, ngrp) )
      param$tau.beta <-
        array(param$tau.beta, c(ncol.tau, ns, ns, nc, nt-1, ngrp) )
      
      if(ncol.del == 1){
        if(ngrp >= 2){
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], c(2,3), mean)
        }
        if(ngrp == 1){
          param$delta.beta <-
            apply(delta.beta.sam[(burn.in+1):niter,,,,], 2, mean)
        }
      }     
      if(ncol.tau == 1){
        if(ngrp >= 2){
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3,4), mean)
        }
        if(ngrp == 1){
          param$tau.beta <-
            apply(tau.beta.sam[(burn.in+1):niter,,,,,,], c(2,3), mean)
        }
      }
      
      param$big.rho <-
        apply(big.rho.sam[(burn.in+1):niter,,,,,], c(2,3,4,5), mean)
      if(ngrp >= 2){
        dimnames(param$delta) <- list(sta.nam, time.nam, str.nam)
        dimnames(param$tau) <- list(sta.nam, sta.nam, str.nam)
      }
      if(ngrp == 1){
        dimnames(param$delta) <- list(sta.nam, time.nam)
        dimnames(param$tau) <- list(sta.nam, sta.nam)
      }
      
      if(ncol.del >= 2){
        dimnames(param$delta.beta) <-
          list(cov.del.nam, sta.nam, cls.nam, str.nam)
      }
      if(ncol.del == 1){
        if(ngrp >= 2){
          dimnames(param$delta.beta) <- list(sta.nam, str.nam)
        }
        if(ngrp == 1){
          names(param$delta.beta) <- sta.nam
        }
      }

      if(ncol.tau >= 2){
        dimnames(param$tau.beta) <-
          list(cov.tau.nam, sta.nam, sta.nam, cls.nam, time.nam[-1], str.nam)
      }
      if(ncol.tau == 1){
        if(ngrp >= 2){
          dimnames(param$tau.beta) <- list(sta.nam, sta.nam, str.nam)        
        }
        if(ngrp == 1){
          dimnames(param$tau.beta) <- list(sta.nam, sta.nam)
        }
      }
      
      dimnames(param$big.rho) <- list(si.nam, NULL, sta.nam, time.nam)
    }
    if( is.cls == T ){
      param$beta <- array(param$beta, c(ncol.gam,nc,ngrp))
      dimnames(param$beta)  <- list(cov.gam.nam, cls.nam, str.nam)
    }
    param$delta.beta <- array(param$delta.beta, c(ncol.del,ns,nc,ngrp))
    dimnames(param$delta.beta)  <- list(cov.del.nam, sta.nam, cls.nam, str.nam)
    param$tau.beta <- array(param$tau.beta, c(ncol.tau,ns,ns,nc,nt-1,ngrp))
    dimnames(param$tau.beta)  <-
      list(cov.tau.nam, sta.nam, sta.nam, cls.nam, time.nam[2:nt], str.nam) 
    ## save confidence intervals
    interval <- list()
    if( is.cls == T ){
      interval$rho.LCI <- rho.LCI
      interval$rho.UCI <- rho.UCI
      interval$gamma.LCI <- gamma.LCI
      interval$gamma.UCI <- gamma.UCI
      interval$beta.LCI <- beta.LCI
      interval$beta.UCI <- beta.UCI
    }
    interval$big.rho.LCI <- big.rho.LCI
    interval$big.rho.UCI <- big.rho.UCI
    interval$delta.LCI <- delta.LCI
    interval$delta.UCI <- delta.UCI
    interval$tau.LCI <- tau.LCI
    interval$tau.UCI <- tau.UCI
    interval$delta.beta.LCI <- delta.beta.LCI
    interval$delta.beta.UCI <- delta.beta.UCI
    interval$tau.beta.LCI <- tau.beta.LCI
    interval$tau.beta.UCI <- tau.beta.UCI   
    ## save sample
    sample <- list()   
    if( is.cls == T ){
      sample$gamma <- gamma.sam
      sample$beta <- beta.sam
      sample$rho <- rho.sam
    }
    sample$delta <- delta.sam
    sample$tau <- tau.sam
    sample$delta.beta <- delta.beta.sam
    sample$tau.beta <- tau.beta.sam
    sample$big.rho <- big.rho.sam
    ## save accpeptance rate
    accept <- list()
    if( is.cls == T ){
      accept$beta <- acp.beta
    }
    accept$delta.beta <- acp.delta.beta
    accept$tau.beta <- acp.tau.beta
    return(param, interval, sample, accept)
  }