## ----lib-----------------------------------------------------------------
library(CNPBayes)

## ----McmcParams----------------------------------------------------------
mp <- McmcParams(iter=1000, burnin=100, thin=1)

## ----graph-marginal, echo=FALSE------------------------------------------
library(grid)
bayesianGrob <- function(x, r=unit(0.25, "inches")){
  tg <- textGrob(x)
  cg <- circleGrob(r=r)
  boxedText <- gTree(children=gList(tg, cg))
}

grobY <- bayesianGrob(expression(y[i]))
grobThetas <- bayesianGrob(expression(theta[h]))
grobSigma2 <- bayesianGrob(expression(sigma[h]^2))
grobH <- bayesianGrob(expression(z[i]))
grobNu0 <- bayesianGrob(expression(nu[0]))
grobSigma02 <- bayesianGrob(expression(sigma[0]^2))
grobPi <- bayesianGrob(expression(pi[h]))
grobMu <- bayesianGrob(expression(mu[h]))
grobTau2 <- bayesianGrob(expression(tau[h]^2))

h <- unit(0.25, "inches")
e <- unit(0.05, "inches")
d <- unit(0.025, "inches")

ar <- arrow(ends="last", length=unit(0.075, "inches"), type="closed")
grid.newpage()
y.x <- 0.5; y.y <- 0.1
h.x <- 0.1; h.y <- 0.3
theta.x <- 0.5; theta.y <- 0.3
sigma2.x <- 0.7; sigma2.y <- 0.3
pi.x <- 0.1; pi.y <- 0.5
mu.x <- 0.45; mu.y <- 0.5
tau2.x <- 0.55; tau2.y <- 0.5
nu0.x <- 0.65; nu0.y <- 0.5
s20.x <- 0.75; s20.y <- 0.5
grid.draw(editGrob(grobY, vp=viewport(y.x, y.y), gp=gpar(fill="gray")))
grid.draw(editGrob(grobY, vp=viewport(y.x, y.y), gp=gpar(fill="transparent")))
grid.draw(editGrob(grobH, vp=viewport(h.x, h.y)))

grid.draw(editGrob(grobThetas, vp=viewport(theta.x, theta.y)))
grid.draw(editGrob(grobSigma2, vp=viewport(sigma2.x, sigma2.y)))

## theta -> y
grid.move.to(unit(theta.x, "npc"), unit(theta.y, "npc") - h)
grid.line.to(unit(y.x, "npc"), unit(y.y, "npc") + h, arrow=ar,
             gp=gpar(fill="black"))
## sigma2 -> y
grid.move.to(unit(sigma2.x, "npc") - e, unit(sigma2.y, "npc") -h)
grid.line.to(unit(y.x, "npc") + h, unit(y.y, "npc") + h, arrow=ar,
             gp=gpar(fill="black"))

## h -> theta
grid.move.to(unit(h.x, "npc") + h, unit(h.y, "npc") - h)
grid.line.to(unit(y.x, "npc") - h, unit(y.y, "npc") + h,  arrow=ar,
             gp=gpar(fill="black"))

##pi
grid.draw(editGrob(grobPi, vp=viewport(pi.x, pi.y)))
## pi -> h
grid.move.to(x=unit(pi.x, "npc"), y=unit(pi.y, "npc") - h)
grid.line.to(x=unit(h.x, "npc"),
             y=unit(h.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))


## mu_h
grid.draw(editGrob(grobMu, vp=viewport(mu.x, mu.y)))
## mu_h -> theta
grid.move.to(x=unit(mu.x, "npc")+e, y=unit(mu.y, "npc") - h)
grid.line.to(x=unit(theta.x, "npc")-e, y=unit(theta.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))


## sigma2_h
grid.draw(editGrob(grobTau2, vp=viewport(tau2.x, tau2.y)))
## sigma2_h -> theta_h
grid.move.to(x=unit(tau2.x, "npc")-e, y=unit(tau2.y, "npc") - h)
grid.line.to(x=unit(theta.x, "npc")+e, y=unit(theta.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))

grid.draw(editGrob(grobNu0, vp=viewport(nu0.x, nu0.y)))
## nu_0 -> sigma2_0
grid.move.to(x=unit(nu0.x, "npc")+e, y=unit(nu0.y, "npc") - h)
grid.line.to(x=unit(sigma2.x, "npc")-e, y=unit(sigma2.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))

grid.draw(editGrob(grobSigma02, vp=viewport(s20.x, s20.y)))
## nu_0 -> sigma2_0
grid.move.to(x=unit(s20.x, "npc")-e, y=unit(s20.y, "npc") - h)
grid.line.to(x=unit(sigma2.x, "npc")+e, y=unit(sigma2.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))

## ----Hyperparameters-----------------------------------------------------
hypp <- Hyperparameters(type="marginal", k=3)

## ----simulateData--------------------------------------------------------
sim.data <- simulateData(N=2500, p=rep(1/3, 3),
                         theta=c(-1, 0, 1),
                         sds=rep(0.1, 3))

## ----model---------------------------------------------------------------
model <- MarginalModel(data=y(sim.data), k=3,
                       hypp=hypp,
                       mcmc.params=mp)

model <- posteriorSimulation(model)

## ----plot----------------------------------------------------------------
plot(DensityModel(model), y(sim.data), 
     main="Marginal Model posterior")

## ----batch---------------------------------------------------------------
## Create McmcParams for batch model
mp <- McmcParams(iter=2000, burnin=1000, thin=1)

## Create Hyperparameters for batch model
hypp <- Hyperparameters(type="batch", k=3)

## simulate batch data
k <- 3
nbatch <- 3
means <- matrix(c(-1.2, -1.0, -0.8,
                -0.2, 0, 0.2,
                0.8, 1, 1.2), nbatch, k, byrow=FALSE)
sds <- matrix(0.1, nbatch, k)
N <- 1500
sim.data <- simulateBatchData(N=N,
                              batch=rep(letters[1:3], length.out=N),
                              theta=means,
                              sds=sds,
                              p=c(1/5, 1/3, 1-1/3-1/5))

# create BatchModel and run posteriorSimulation
model <- BatchModel(data=y(sim.data), k=3,
                    batch=batch(sim.data),
                    hypp=hypp,
                    mcmc.params=mp)

model <- posteriorSimulation(model)
plot(DensityModel(model), y(sim.data),
     breaks=100)

## ----marginalLikelihood--------------------------------------------------
mp <- McmcParams(iter=2e3, burnin=1e3)
model <- MarginalModel(data=y(sim.data), k=2, mcmc.params=mp)

mlist <- list(posteriorSimulation(model),
              posteriorSimulation(model, k=3))

x <- marginalLikelihood(mlist)
x

## ----bic-----------------------------------------------------------------
## simulate k=2 model
sim.data <- simulateData(N=2500, p=rep(1/3, 3),
                         theta=c(-1, 0, 1),
                         sds=rep(0.1, 3))
hypp1 <- Hyperparameters(k=2)
m1 <- MarginalModel(data=y(sim.data), k=2,
                    hypp=hypp1,
                    mcmc.params=mp)
m1 <- posteriorSimulation(m1)

## simulate k=3 model
hypp2 <- Hyperparameters(k=3)
m2 <- MarginalModel(data=y(sim.data), k=3,
                    hypp=hypp2,
                    mcmc.params=mp)
m2 <- posteriorSimulation(m2)

bic(m1)
bic(m2)

## ----merge---------------------------------------------------------------
mp <- McmcParams(iter=5000, burnin=1000, thin=1)
model <- MarginalModel(data=y(sim.data), k=4,
                       hypp=Hyperparameters(k=4),
                       mcmc.params=mp)
model <- posteriorSimulation(model)
dm <- DensityModel(model)
modes(dm)
dm_merged <- DensityModel(model, merge=TRUE)
k(dm_merged)

par(mfrow=c(1,2), las=1)
plot(dm, y(sim.data))
plot(dm_merged, y(sim.data))

## ----merge-under---------------------------------------------------------
model <- MarginalModel(data=y(sim.data), k=3,
                       hypp=Hyperparameters(k=3),
                       mcmc.params=mp)
model <- posteriorSimulation(model)
dm <- DensityModel(model)
dm_merged <- DensityModel(model, merge=TRUE)
k(dm)
k(dm_merged)

## ----collapseBatch-------------------------------------------------------
k <- 3
nbatch <- 3
means <- matrix(c(-1.2, -1.0, -1.0,
                -0.2, 0, 0,
                0.8, 1, 1), nbatch, k, byrow=FALSE)
sds <- matrix(0.1, nbatch, k)
N <- 1500
sim.data <- simulateBatchData(N=N,
                              batch=rep(letters[1:3], length.out=N),
                              theta=means,
                              sds=sds,
                              p=c(1/5, 1/3, 1-1/3-1/5))

plot(sim.data)

model <- BatchModel(data=y(sim.data), k=3,
                    batch=batch(sim.data))

model <- BatchModel(data=y(sim.data), k=3,
                    batch=collapseBatch(model))

model <- posteriorSimulation(model)
plot(model)