## load function
source("../../../pgm/r_function/ltacov_mcmc_multi.r")
## Read data
cov <- scan("../../../PreAnalysis/winLTA/cov.dat")
cov <- matrix( cov, ncol=3, byrow=T )
w <- cov[,1]==16
cov <- cov[!w,]
items <- scan("../../../PreAnalysis/winLTA/murray.dat")
items <- matrix( items, ncol=12, byrow=T )
items <- items[!w,]
data <- cbind(cov, items)
rm(cov, items)

## Tranfer data
age <- data[,1]
sex <- data[,2]  # 1="male", 2="female"
race <- data[,3] # 1="white", 2="black", 3="other"
w <- sex == 2
sex[!w] <- 0 # 0 = male
sex[w] <- 1  # 1 = female
rm(w)
w <- race == 2
race[!w] <- 0 # 0 = white & other ( white=459, other=80)
race[w] <- 1  # 1 = black (black=922)

ci <- data[,4:5]
si <- array( c(data[,6:15]), c(nrow(ci), 5, 2) )
cov.gamma <- matrix( 1, nrow(ci), 4 )
cov.gamma[,2] <- age
cov.gamma[,3] <- sex
cov.gamma[,4] <- race
cov.delta <- matrix( 1, nrow(ci), 4 )
cov.delta[,2] <- age
cov.delta[,3] <- sex
cov.delta[,4] <- race
cov.tau <- array( 1, c(nrow(ci), 3, 2) )
cov.tau[,2,1] <- age
cov.tau[,2,2] <- age
#cov.tau[,2,1] <- sex
#cov.tau[,2,2] <- sex
cov.tau[,3,1] <- race
cov.tau[,3,2] <- race
rm(data,age,sex,race,w)

## Configure model
ngrp <- 1
nc <- 2
ns <- 4
nt <- 2
nci <- 2
nrc <- c(2, 2)
nsi <- 5
nrs <- c(2, 2, 2, 2, 2)
nkase <- 1461

## Constraints ##
is.big.rho.equal.across.time <- T
is.big.rho.equal.across.cls <- T
big.rho.const <- array(1, c(nsi, max(nrs), ns, nc, nt) )
tau.const <- array(1, c(ns, ns, nc, nt-1, ngrp))
tau.const[3,1,1,1,1] <- 0
tau.const[4,1,1,1,1] <- 0
tau.const[3,2,1,1,1] <- 0
tau.const[4,2,1,1,1] <- 0
tau.const[1,3,1,1,1] <- 0
tau.const[2,3,1,1,1] <- 0
tau.const[1,4,1,1,1] <- 0
tau.const[2,4,1,1,1] <- 0

tau.const[2,1,2,1,1] <- 0
tau.const[3,1,2,1,1] <- 0
tau.const[4,1,2,1,1] <- 0
tau.const[3,2,2,1,1] <- 0
tau.const[4,2,2,1,1] <- 0
tau.const[1,3,2,1,1] <- 0
tau.const[2,3,2,1,1] <- 0
tau.const[1,4,2,1,1] <- 0
tau.const[2,4,2,1,1] <- 0

## Variable names
cls.nam <- c("wk.est", "str.est")
ci.nam <- c("ygdmath", "ygdother")
sta.nam <- c("bad.bad", "bad.good", "good.bad", "good.good")
si.nam <- c("mgrdfn", "egrdfn", "ysimport", "ysnolern", "ysnotime")
time.nam <- c("time.1", "time.2")
cov.gam.nam <- c("int", "age", "sexF", "raceBlk")
ncol.gam <- 4
cov.del.nam <- c("int", "age", "sexF", "raceBlk")
ncol.del <- 4
cov.tau.nam <- c("int", "age", "raceBlk")
ncol.tau <- 3
str.nam <- c("str")

## priors for big rho
prior.big.rho <- array(0.5, c(nsi, max(nrs), ns, nc, nt))
prior.rho <- array(0.5, c(nci, max(nrc), nc))

## prior for betas
prior.beta <- matrix(0, nc, ngrp)
prior.delta.beta <- array(0, c(ns, nc, ngrp))
prior.tau.beta <- array(0, c(ns, ns, nc, ngrp))
################################
## Data Augmentation via MCMC ##
################################
n.fixed <- array(0, c(ns, ns, nc, nt-1, ngrp))
n.fixed[4,3,1,1,1] <- 0
n.fixed[3,4,1,1,1] <- 0
w <- tau.const==0
n.fixed[w] <- 0

niter <- 7000
burn.in <- 2000
t.df <- 4
scale <- 2.4
scale.delta <- matrix( c(1.9, 1.7), nc, ngrp )
scale.tau <- array( c(2.4, 2.6, 
                      2.4, 1.5, 
                      2.1, 2.4,
                      2.4, 2.5), c(nc, ns, ngrp))
update.cov.tau.beta <- F
update.every <- 250

mcmc.17.tau <-
  ltacov.mcmc.multi(ci=ci, si=si, cov.gamma=cov.gamma, cov.delta=cov.delta,
                    cov.tau=cov.tau, nrc=nrc, nrs=nrs,
                    nc=nc, ns=ns, niter=niter, burn.in=burn.in,
                    t.df=t.df, scale=scale, scale.delta=scale.delta,
                    scale.tau=scale.tau,
                    is.big.rho.equal.across.time=is.big.rho.equal.across.time,
                    is.big.rho.equal.across.cls=is.big.rho.equal.across.cls,
                    tau.const=tau.const, big.rho.const=big.rho.const,
                    rho=em.17.tau$param$rho,
                    big.rho=em.17.tau$param$big.rho,
                    beta=em.17.tau$param$beta,
                    delta.beta=em.17.tau$param$delta.beta,
                    tau.beta=em.17.tau$param$tau.beta,
                    prior.rho=prior.rho,
                    prior.big.rho=prior.big.rho,
                    prior.beta=prior.beta,
                    prior.delta.beta=prior.delta.beta,
                    prior.tau.beta=prior.tau.beta,
                    cov.beta=em.17.tau$cov$beta, 
                    cov.delta.beta=em.17.tau$cov$delta.beta, 
                    cov.tau.beta=em.17.tau$cov$tau.beta, 
                    cls.nam=cls.nam, sta.nam=sta.nam, ci.nam=ci.nam,
                    si.nam=si.nam, time.nam=time.nam, 
                    cov.gam.nam=cov.gam.nam, cov.del.nam=cov.del.nam,
                    cov.tau.nam=cov.tau.nam,
                    str.nam=str.nam,
                    update.cov.tau.beta=update.cov.tau.beta,
                    update.every=update.every,
                    assign.sub=em.17.tau$assign.sub,
                    n.fixed=n.fixed)