###################################################
### chunk number 5: periodicPosteriorBeta
###################################################
rm(list = ls())

library(NMMAPSdata)
source(file.path("data", "utils.R"))
library(tlnise)
library(ellipse)
library(MASS)

exclude <- c("anch", "hono")

## Load the results
load(file.path("results", "seasonal.periodic-ortho.lag.012.pm10.rda"))
r <- results[[2]][setdiff(names(results[[2]]), exclude)] 

## Get region indicators for each city
data(cities)
cityRegion <- with(cities, region[match(names(r), city)])
regionNames <- sort(unique(cityRegion))
regionInd <- 1 * outer(cityRegion, regionNames, "==")
colnames(regionInd) <- regionNames

## Pool estimates
bc <- extractBetaCov(r, pollutant = "pm10")
initTLNise()
g <- tlnise(bc$beta, bc$cov, regionInd, prnt = FALSE, labelY = names(r),
            intercept = FALSE, seed = 10)

## Get sample from posterior distribution
set.seed(20)
gs <- tlnise:::sampleGamma(g, bc$cov, bc$beta, regionInd)
set.seed(30)
ds <- drawSample0(gs, n = 1500)
ads <- array(t(ds), c(7, 3, nrow(ds)))

## Overall estimate
initTLNise()
g0 <- tlnise(bc$beta, bc$cov, rep(1, nrow(bc$beta)), seed = 100, prnt = FALSE)
set.seed(200)
gs0 <- tlnise:::sampleGamma(g0, bc$cov, bc$beta)
set.seed(300)
ds0 <- drawSample0(gs0, n = 1500)

## Compute posterior probs.
postprob1 <- matrix(nrow = 8, ncol = 2)
postprob2 <- matrix(nrow = 8, ncol = 2)
n <- dim(ads)[3]

for(i in 1:7) {
    postprob1[i, ] <- c(sum(ads[i, 2, ] > 0) / n, sum(ads[i, 2, ] < 0) / n)
    postprob2[i, ] <- c(sum(ads[i, 3, ] > 0) / n, sum(ads[i, 3, ] < 0) / n)
}
postprob1[8, ] <- c(sum(ds0[, 2] > 0)/nrow(ds0), sum(ds0[, 2] < 0)/nrow(ds0))
postprob2[8, ] <- c(sum(ds0[, 3] > 0)/nrow(ds0), sum(ds0[, 3] < 0)/nrow(ds0))

## Draw ellipses
x1 <- NULL
y1 <- NULL
x2 <- NULL
y2 <- NULL

for(i in 1:7) {
    V <- cov(t(ads[i, -1, ]))
    m <- colMeans(t(ads[i, -1, ]))
    e <- ellipse(V, centre = m, level = c(.75, .95))
    e1 <- t(matrix(e[,1], nrow = 2))
    e2 <- t(matrix(e[,2], nrow = 2))
    x1 <- c(x1, e1[,1])
    y1 <- c(y1, e2[,1])
    x2 <- c(x2, e1[,2])
    y2 <- c(y2, e2[,2])    
}
V0 <- cov(ds0[, -1])
m0 <- colMeans(ds0[, -1])
e <- ellipse(V0, centre = m0, level = c(.75, .95))
e1 <- t(matrix(e[,1], nrow = 2))
e2 <- t(matrix(e[,2], nrow = 2))
x1 <- c(x1, e1[,1])
y1 <- c(y1, e2[,1])
x2 <- c(x2, e1[,2])
y2 <- c(y2, e2[,2])

b <- vector("list", length = 8)

## Only plot a few points
set.seed(1025301)
for(i in 1:7) {
    b[[i]] <- t(ads[i, 2:3, sample(1500, 100)])
}
set.seed(20)
b[[8]] <- ds0[sample(1500, 15), 2:3]

regionFullNames <- c("Industrial Midwest", "Northeast", "Northwest",
                     "Southern California", "Southeast", "Southwest",
                     "Upper Midwest", "All Regions")
f <- gl(8, 50, ordered = TRUE)
levels(f) <- regionFullNames

library(lattice)
## trellis.device(dev = x11, height = 5, width = 8, color = FALSE)

yrng <- c(-.1, .1)
xrng <- c(-.1, .09)
p <- xyplot(y2 ~ x2 | f, subscripts = TRUE, type = "l", lty = "43", lwd = 1.5,
            panel = function(x, y, subscripts, ...) {
                i <- (subscripts[1] - 1) / 50 + 1
                lpoints(b[[i]], pch = 1, cex = 0.5) ## col = gray(0.5))
                panel.xyplot(x, y, ...)
                panel.abline(h = 0, v = 0, lty = 3, lwd = 1)
                llines(x1[subscripts], y1[subscripts], lty = 1, col = 1, lwd = 1.5)
                ltext(-0.05, -0.075,
                      labels = eval(substitute(expression(paste(italic(p), "(", beta[2] > 0, ")") == k), list(k = formatC(round(postprob2[i, 1], 2), digits = 2, flag = "#")))), cex = 0.7)
                ltext(-0.05, -0.09,
                      labels = eval(substitute(expression(paste(italic(p), "(", beta[1] > 0, ")") == k), list(k = formatC(round(postprob1[i, 1], 2), digits = 2, flag = "#")))), cex = 0.7)
            }, as.table = TRUE, ylim = yrng, aspect = 1, xlim = xrng,
            xlab = expression(beta[1]), ylab = expression(beta[2]),
            scales = list(cex = 0.8, alternating = 1, tck = c(0.8, 0)),
            par.strip.text = list(cex = 0.8),
            strip = strip.custom(bg = 0)
            )
print(p)