demo/norm-ptt.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/poisNormal.R
## 
## Normal imprecise prior on the use of models without covariates
##
## 2013.10.20

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

# zero-truncated random variates from Poisson with mean 1
y <- rvg4yx(N=5e0, pars=1, ztrunc=TRUE, xreg=FALSE, kind="pois")$y
table(y)

## imprecise prior specification
## eqn1. -6 <= nu <= 6
## eqn2. 0.5 <= tau^2 <= 2.5
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-6, -6, 0.5, -2.5))
m2fit <- model(formula=y~0, ztrunc=TRUE, dist="poisson")
cmfit <- impose(obj=m2fit, eqns=lc0)
plot(cmfit, xlim=c(-7,7), ylim=c(-0.5,3.5), 
  main="Hyperparameter Space\nconstrained by a set of inequalities",
  xlab=expression(alpha), ylab=expression(beta))

## model fitting and its summary
## try to change 'method=LA' to 'LA2', 'MH', or 'IS'.
op <- update(obj=cmfit, method="LA1", apriori="normal")
sop <- summary(op, HT.est=TRUE)


## visualization - imprecise posterior
z <- range(sop$est)
mar <- 0.3
plot(sop, xlim=c(-7,7), ylim=c(0,3), zlim=c(z[1]-mar, z[2]+mar),
  xlab=list(expression(nu), cex=1.5, font=1), 
  ylab=list(expression(tau^2), cex=1.5, font=1), 
  zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90, cex=1.5, font=1),
  main="Poisson with Imprecise Normal Prior")

# probability box
pbox(sop, which.xtms=c(sop$infs.idx, sop$sups.idx), smooth=TRUE)


## visualization (on grid) - imprecise prior
cmfit1 <- setGrid(obj=cmfit, len=20)
plot(cmfit1, xlim=c(-7,7), ylim=c(-0.5,3.5), 
  main="Imprecise Prior Specification",
  xlab=expression(alpha), ylab=expression(beta))


## visualization (on grid) - imprecise posterior
## try to change 'method=LA' to 'LA2', 'MH', or 'IS'.
op1 <- update(obj=cmfit1, method="LA1", apriori="normal")
sop1 <- summary(op1)
z <- range(sop1$est)
mar <- 0.3
plot(sop1, xlim=c(-7,7), ylim=c(0,3), zlim=c(z[1]-mar, z[2]+mar),
  xlab=list(expression(nu), cex=1.5, font=1), 
  ylab=list(expression(tau^2), cex=1.5, font=1), 
  zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90, cex=1.5, font=1),
  main="Poisson with Imprecise Normal Prior")

# probability box
pbox(sop1, which.xtms=c(sop1$infs.idx, sop1$sups.idx), smooth=TRUE)

## Other imprecise prior specification and visualization of imprecise posterior
## eqn1. tau^2 >= 0
## eqn2. -nu -tau^2 >= -5
## eqn3. -nu -tau^2 >= -5
lc1 <- list(lhs=rbind(c(0,1), c(-1,-2), c(1,-2)), rhs=c(0, -5, -5))
bounds1 <- impose(eqns=lc1)
plot(sop1, bnd=bounds1, xlim=c(-7,7), ylim=c(0,3), zlim=c(z[1]-mar, z[2]+mar),
  xlab=list(expression(nu), cex=1.5, font=1), 
  ylab=list(expression(tau^2), cex=1.5, font=1), 
  zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90, cex=1.5, font=1),
  main="Poisson with Imprecise Normal Prior")