demo/normal-tsl.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## 2013.11.11
## THIS PROGRAM IS PROCESSED INTO THE CHAPTER 4 OF THE THESIS 

## imprecise/demo/normal-translation.R
## 
## translation behaviour of a canonical parameter on its hyperparameter space.
## -- normal imprecise prior
#
# hparam <- c(nu, tau2)
#           eta2         eta1          eta0
# eta <- c(0.5/tau2, sum(y)+nu/tau2, length(y))
#          theta^2       theta       exp(theta)
# -6 <= nu <= 6
# 0.5 <= tau2 <= 2.5
# y <- rvg4yx(N=5e1, pars=1, ztrunc=TRUE, xreg=FALSE, kind="pois")$y

library(ipeglim)
rm(list=ls())

lmu <- c(-3, 2, 1, 0, -1, -2, 3)
e0 <- seq(0,2,0.25)
e2 <- seq(0,5,0.3)
hcube <- expand.grid(e0=e0, e2=e2, lmu=lmu)

e1 <- apply(hcube, 1, function(x){
  x <- as.vector(x)
  rt.e1 <- tryCatch(uniroot(eta1fn, lower=-20, upper=30, tol=1e-4, eta0=x[1], eta2=x[2], logm=x[3])$root, error=function(e) return(NaN))
  return(rt.e1)
  })
hcube$e1 <- e1

# n=0
cube0 <- subset(hcube, subset=(e0==0))
with(cube0, plot(e1, e2, type='l', xlim=c(-20,20),
  main="Translation Behaviour\n(Normal Imprecise Prior on Poisson Model)",
  xlab=expression(eta[1]),
  ylab=expression(eta[2]), col="black", lty=1
))
box0 <- rbind(c(-2,0.2), c(2,0.2), c(2,1.0), c(-2,1.0))
polygon(box0, border="black", lty=1, lwd=3)
mlab <- subset(cube0, subset=(e2=="2.7"))
with(mlab, text(x=e1, y=e2, labels=lmu)) 

# n=1
cube1 <- subset(hcube, subset=(e0==1))
with(cube1, lines(e1, e2, col="blue", lty=2))
y1 <- rpois(n=1, lambda=1)
box1 <- box0 + matrix(rep(c(y1,0),nrow(box0)), ncol=2, byrow=TRUE)
polygon(box1, border="blue", lty=2, lwd=3)

# n=2
cube2 <- subset(hcube, subset=(e0==2))
with(cube2, lines(e1,e2,col="red", lty=3))
y2 <- c(y1, rpois(n=1, lambda=1))
box2 <- box1 + matrix(rep(c(sum(y2),0), nrow(box1)), ncol=2, byrow=TRUE)
polygon(box2, border="red", lty=3, lwd=3)

legend(x=-20, y=1, 
  legend=c("n=0", "n=1", "n=2"),
  col=c("black", "blue", "red"),
  lty=c(1,2,3))
  

########################################################################
## Scenario 2 -- important
# hparam <- c(nu, tau2)
#           eta2         eta1          eta0
# eta <- c(0.5/tau2, sum(y)+nu/tau2, length(y))
#          theta^2       theta       exp(theta)
# -6 <= nu <= 6
# 0.5 <= tau2 <= 2.5

rm(list=ls())
y <- rpois(n=10, lambda=2)

lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-6, -6, 0.5, -2.5))

y1 <- y[1]
m2fit <- model(formula=y1~0, ztrunc=FALSE, dist="poisson")
cmfit <- iprior(obj=m2fit, eqns=lc0)
op <- summary(update(obj=cmfit, apriori="normal", method="LA", silent=TRUE), silent=TRUE)

mat0 <- as.data.frame(do.call(rbind, op$xtms)) # nspace
names(mat0) <- c("nu", "tau2")
xtms0 <- with(mat0, cbind(0.5/tau2, nu/tau2)) # hspace

plot(xtms0, xlim=c(0,1.5), ylim=c(-15, 15), 
  xlab=expression(eta[2]), ylab=expression(eta[1]))
polygon(xtms0[chull(xtms0),])

mat1 <- op$tslpars2
xtms1 <- op$tslpars[,1:2]
polygon(xtms1[chull(xtms1),], border="blue")

y2 <- y[1:2]
m2fit <- model(formula=y2~0, ztrunc=FALSE, dist="poisson")
cmfit <- iprior(obj=m2fit, mat=mat1)
op <- summary(update(obj=cmfit, apriori="normal", method="LA", silent=TRUE), silent=TRUE)
xtms2 <- op$tslpars[,1:2]
mat2 <- op$tslpars2
polygon(xtms2[chull(xtms2),], border="red")

y3 <- y[1:3]
m2fit <- model(formula=y3~0, ztrunc=FALSE, dist="poisson")
cmfit <- iprior(obj=m2fit, mat=mat2)
op <- summary(update(obj=cmfit, apriori="normal", method="LA", silent=TRUE), silent=TRUE)
xtms3 <- op$tslpars[,1:2]
mat3 <- op$tslpars2
polygon(xtms3[chull(xtms3),], border="darkgreen")