demo/xreg-ns.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/xreg-ns.R
## 
## * Test a numberical stability between different numerical methods of 'LA',
## 'MH', and 'IS'.
## * Measure a computing time on each numerical method
##
## 2013.10.22

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

## default structure of X
xm <- c(0)
xsd <- c(1)

## default structure of B
bm <- c(0,0)
bsd <- c(1,1)
br <- diag(2)
br[1,2] <- br[2,1] <- r1 <- 0
Bmat <- t(diag(bsd))%*%br%*%diag(bsd)
Bmat <- Bmat*1e-2

# parameter specification
bcfs1 <- c(0.5, 0.5)      # 1 quadrant
bcfs2 <- c(-0.5, 0.5)     # 2 quadrant
bcfs3 <- c(-0.5, -0.5)    # 3 quadrant
bcfs4 <- c(0.5, -0.5)     # 4 quadrant

# random variates generation
tmp <- rvg4yx(N=1e2, cf0=bcfs1[1], pars=bcfs1[2], ztrunc=TRUE, xreg=TRUE, kind="mvnorm", mean=xm, sd=xsd, center.X=TRUE)
dat <- tmp$yX
table(dat$y)
xi1 <- c(1,unlist(dat[1,-1]))

m2fit <- model(formula=y~x1, data=dat, ztrunc=TRUE, dist="poisson")
cmfit <- iprior(obj=m2fit, circle=list(x0=0, y0=0, r=1, len=20))

t0 <- proc.time()

ola <- summary(update(obj=cmfit, method="LA", B=Bmat, xi=xi1))
t1 <- proc.time()

omh <- summary(update(obj=cmfit, method="MH", B=Bmat, xi=xi1))
t2 <- proc.time()

ois <- summary(update(obj=cmfit, method="IS", B=Bmat, xi=xi1))
t3 <- proc.time()

tt <- list(t0, t1, t2, t3)
elapsed <- diff(do.call(c, lapply(tt, "[[", 3)))
names(elapsed) <- c("LA", "MH", "IS")
elapsed

tb.x1 <- cbind(ola$est$x1, omh$est$x1, ois$est$x1)
tb.x2 <- cbind(ola$est$x2, omh$est$x2, ois$est$x2)
tb.x1
tb.x2


X11()
plot(ola,
  main="Laplace Approximation", 
  sub=list(expression(-1<beta[0]~","~beta[1]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))

X11()
plot(omh,
  main="Metropolis-Hastings Algorithm", 
  sub=list(expression(-1<beta[0]~","~beta[1]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))

X11()
plot(ois,
  main="Importance Sampling",
  sub=list(expression(-1<beta[0]~","~beta[1]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))