### This file reproduces example 9.3 from Bayes Nets in Educational
### Assessment by Almond, Mislevy, Steinberg, Williamson, and Yan


Theta <- varDescription("theta","theta",role="SM", levels=c("1","0"))
X2 <- varDescription("X2","X2",role="Observable", levels=c("1","0"))
X3 <- varDescription("X3","X3",role="Observable", levels=c("2","1","0"))

## Cant quite figure how to get all the priors into one table.
##priorTable <- read.csv("priorTable.csv")

Theta.ptable <- matrix(c(3,3),nrow=1,dimnames=list(NULL,c("1","0")))
Theta.dtable <- Theta.ptable/rowSums(Theta.ptable)
## Workaround for "feature" of R that coerces names of data frames to
## be proper identifiers.
Theta.ptf <- data.frame(Theta.ptable)
Theta.dtf <- data.frame(Theta.dtable)
names(Theta.ptf) <- c("1","0")
names(Theta.dtf) <- c("1","0")


Theta.dist <- distribution(cond=character(0), cons=Theta$nodeName,
                           type="HyperDirichlet Distribution",
                           table=Theta.dtf, parameterTable = Theta.ptf)

lc2.sm <- gmModel("ecd://bninea/exMC/lc2/", type="SM", var=list(Theta),
               dist=list(Theta.dist))

X2.ptable <- matrix(c(2,4,
                      4,2),
                    nrow=2,byrow=TRUE,
                    dimnames=list(c("0","1"),c("1","0")))
X2.dtable <- X2.ptable/rowSums(X2.ptable)
X2.ptf <- data.frame(theta=c("0","1"), X2.ptable)
X2.dtf <- data.frame(theta=c("0","1"), X2.dtable)
## Workaround for "feature" of R that coerces names of data frames to
## be proper identifiers.
names(X2.ptf) <- c("theta","1","0")
names(X2.dtf) <- c("theta","1","0")

X2.dist <- distribution(cond=Theta$nodeName, cons=X2$nodeName,
                           type="HyperDirichlet Distribution",
                           table=X2.dtf, parameterTable = X2.ptf)

mc.em <- gmModel("ecd://bninea/exMC/binary/","ecd://bninea/exMC/MC/",
               type="EM",
               var=list(Theta,X2),
               dist=list(X2.dist))


X3.ptable <- matrix(c(1,2,3,
                      3,2,1),
                    nrow=2,byrow=TRUE,
                    dimnames=list(c("0","1"),c("2","1","0")))
X3.dtable <- X3.ptable/rowSums(X3.ptable)
X3.ptf <- data.frame(theta=c("0","1"), X3.ptable)
X3.dtf <- data.frame(theta=c("0","1"), X3.dtable)
## Workaround for "feature" of R that coerces names of data frames to
## be proper identifiers.
names(X3.ptf) <- c("theta","2","1","0")
names(X3.dtf) <- c("theta","2","1","0")

X3.dist <- distribution(cond=Theta$nodeName, cons=X3$nodeName,
                           type="HyperDirichlet Distribution",
                           table=X3.ptf, parameterTable=X3.dtf)

pc3.em <- gmModel("ecd://bninea/exMC/level3/","ecd://bninea/exMC/PC3/",
               type="EM",
               var=list(Theta,X3),
               dist=list(X3.dist))

mcLatent <-
  amd(list("lc2_SM"=lc2.sm),list("mc_EM"=mc.em,"pc3_EM"=pc3.em),list(),
      missingCodes=defaultMissingCodes)

mcLatent$stat <- addStatsForVars(mcLatent,c("margin"="Bayes Net Margin",
                                              "MAP"="Bayes Net Mode"))

writeAMD(mcLatent)


dataRep <- read.csv("dataTable.csv")
sum(dataRep$Count)
dataUnrolled <- do.call(rbind,
                        apply(dataRep,1,function (row)
                              matrix(rep(c(row[1],NA,NA,NA,row[2:4]),row[5]),
                                         ncol=7,byrow=TRUE)))
colnames(dataUnrolled) <- c("TrueTheta","MAPTheta",
                            "marginTheta0","marginTheta1",
                            "X1","X2","X3")
write.table(dataUnrolled,"dataUnrolled.csv",row.names=TRUE,col.names=FALSE,
            sep=",")


newDataRep <- read.csv("newDataTable.csv")
sum(newDataRep$Count)
newDataUnrolled <- do.call(rbind,
                        apply(newDataRep,1,function (row)
                              matrix(rep(c(NA,NA,NA,row[1:3]),row[4]),
                                         ncol=6,byrow=TRUE)))
colnames(newDataUnrolled) <- c("MAPTheta",
                            "marginTheta0","marginTheta1",
                            "X1","X3","X4")
write.table(newDataUnrolled,"newDataUnrolled.csv",row.names=TRUE,
            col.names=FALSE, sep=",")