demo/lg-ptt.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## imprecise/demo/lgamma.R
##
## illustration how to perform the imprecise inferential framework by
## the use of log-gamma imprecise prior on the canonical parameter of 
## Poisson model.
##
## 2013.10.21.

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

## random variates generation from the Poisson with mean 'lambda'
y <- rvg4yx(N=5e1, pars=1, ztrunc=TRUE, xreg=FALSE, kind="pois")$y
table(y)

## model to be fitted
m2fit <- model(formula=y~0, ztrunc=TRUE, dist="pois")
attributes(m2fit)

## model fit with the imprecise prior  
## 0.05 <= alpha <= 10
## 0.05 <= beta <= 10
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(0.05, -10, 0.05, -10))
cmfit <- impose(obj=m2fit, eqns=lc0)
attributes(cmfit)

## visualization - imprecise prior
## hyperparameter space (M0) 
plot(cmfit, xlim=c(-2,12), ylim=c(-2,12), 
  main="Hyperparameter Space\n constrained by a set of inequalities",
  xlab=expression(alpha), ylab=expression(beta))

## output from model fit
op <- update(obj=cmfit, method="LA", apriori="lgamma")
attributes(op)

## summarized output 
sop <- summary(op, HT.est=TRUE)
attributes(sop)

## quality check
sop$sratio

## print method
sop

## show the probability box
pbox(x=sop, smooth=TRUE, 
  main="Probability Box of Canonical Parameter")

pbox(x=sop, which.xtms=c(sop$infs.idx, sop$sups.idx), smooth=TRUE, 
  main="Probability Box of Canonical Parameter")

## visualization (minimal) - imprecise posterior
z <- range(sop$est)
mar <- 0.5
plot(sop, xlim=c(-2,12), ylim=c(-2,12), zlim=c(z[1]-mar, z[2]+mar),
  xlab=list(expression(alpha), cex=1.5, font=1), 
  ylab=list(expression(beta), cex=1.5, font=1), 
  zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90, cex=1.5, font=1),
  main="Poisson with Imprecise Log-Gamma Prior")

## visualization (on grid) - imprecise prior
cmfit1 <- setGrid(cmfit, len=10)
plot(cmfit1, xlim=c(-2,12), ylim=c(-2,12), 
  main="Hyperparameter specification",
  xlab=expression(alpha), ylab=expression(beta))


## try to change 'method=LA' to 'LA2', 'IS', or 'MH'.
op1 <- update(obj=cmfit1, method="LA", apriori="lgamma")
sop1 <- summary(op1, HT.est=TRUE)
sop1

## visualization (on grid) - imprecise posterior
z <- range(sop$est)
mar <- 0.5
plot(sop1, xlim=c(-2,12), ylim=c(-2,12), zlim=c(z[1]-mar, z[2]+mar),
  xlab=list(expression(alpha), cex=1.5, font=1), 
  ylab=list(expression(beta), cex=1.5, font=1), 
  zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90, cex=1.5, font=1),
  main="Surface of Canonical Parameter\nLog-Gamma Imprecise Prior")

# probability box
pbox(sop1, which.xtms=c(sop1$infs.idx, sop1$sups.idx, 31, 59), smooth=TRUE,
  main="Probability Box of Canonical Parameter\non the use of Log-Gamma Imprecise Prior")

## Other imprecise prior specification and its imprecise posterior visualization
bounds <- impose(circle=list(x0=5,y0=5, r=5, len=20))
plot(sop1, bnd=bounds, xlim=c(-2,12), ylim=c(-2,12), zlim=c(z[1]-mar, z[2]+mar),
  xlab=list(expression(alpha), cex=1.5, font=1), 
  ylab=list(expression(beta), cex=1.5, font=1), 
  zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90, cex=1.5, font=1),
  main="Surface of Canonical Parameter\non constratined log-gamma imprecise prior hyperparameter space")