demo/lg-ns.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## imprecise/demo/lgamma-ns.R
##
## show the numerical stabiability on estimating the canonical parameter of
## Poisson model using numerical algorithms of 'LA', 'MH', 'IS', 'AQ',
## as well as 'NA'.
##
## 2013.10.21.

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

# random variates generation
### try to change 'ztrunc' from 'TRUE' to 'FALSE'
y <- rvg4yx(N=5e1, pars=1, ztrunc=TRUE, xreg=FALSE, kind="pois")$y
table(y)

# model to be fitted 
### the value of 'ztrunc' in the function 'model()' should be consistent with 
### the value of 'ztrunc' in the function 'rvg4yx()'. 
m2fit <- model(formula=y~0, ztrunc=TRUE, dist="pois")

# box-constrained log-gamma imprecise prior
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)

# apply numerical methods 
t0 <- proc.time()

oaq <- summary(update(obj=cmfit, method="AQ", apriori="lgamma"))
t1 <- proc.time()

ola <- summary(update(obj=cmfit, method="LA", apriori="lgamma"))
t2 <- proc.time()

omh <- summary(update(obj=cmfit, method="MH", apriori="lgamma"))
t3 <- proc.time()

ois <- summary(update(obj=cmfit, method="IS", apriori="lgamma"))
t4 <- proc.time()

tt <- list(t0, t1, t2, t3, t4)
elapsed <- diff(do.call(c, lapply(tt, "[[", 3)))
tb <- cbind(oaq$est, ola$est, omh$est, ois$est, elapsed)
colnames(tb) <- c("AQ", "LA", "MH", "IS", "Time")
printCoefmat(tb, tst.ind=1:5)

## quality check 
ois$ess     # 'IS'
ola$arate   # 'MH'
ois$sratio  # 'LA'