demo/final.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/xreg-cc2d-x0.R
# intercept model (final test)
# 2013.10.24

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

source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/cpef2reg.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/cpef.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/model.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/iprior.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/update_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/summary_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/print_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/plot_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/plot_summary.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/pbox_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/qc.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/rvg4yx.R")
# iprior -> imprecisePrior

## case 1. log(E(y)) ~ b0 + b1*x1
## 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


SILENCE <- FALSE
ZTRUNC <- TRUE
VERBOSE <- TRUE

tmp <- rvg4yx(N=5e2, 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=ZTRUNC, dist="poisson")
cmfit <- iprior(obj=m2fit, circle=list(x0=0, y0=0, r=1, len=5))

t0 <- proc.time()

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

# LAF is not working with the model including the intercept
olaf <- update(obj=cmfit, method="LAF", B=Bmat, xi=xi1, silent=SILENCE)
solaf <- summary(olaf, verbose=VERBOSE)
t2 <- proc.time()

omh <- update(obj=cmfit, method="MH", B=Bmat, xi=xi1, silent=SILENCE)
somh <- summary(omh, verbose=VERBOSE)
t3 <- proc.time()

ois <- update(obj=cmfit, method="IS", B=Bmat, xi=xi1, silent=SILENCE)
sois <- summary(ois, verbose=VERBOSE)
t4 <- proc.time()

tb1 <- cbind(somh$est[,1], sola$est[,1], sois$est[,1], solaf$est[,1])
tb1
tb2 <- cbind(somh$est[,2], sola$est[,2], sois$est[,2], solaf$est[,2])
tb2

X11()
plot(sola, # for m0xtms
  main="Imprecise Posterior Shrinkage - Laplace", 
  sub=list(expression(-1<beta[1]~","~beta[2]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))

X11()
plot(somh, # for m0xtms
  main="Imprecise Posterior Shrinkage - Metropolis", 
  sub=list(expression(-1<beta[1]~","~beta[2]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))

X11()
plot(sois, # for m0xtms
  main="Imprecise Posterior Shrinkage - Importance", 
  sub=list(expression(-1<beta[1]~","~beta[2]<1), cex=1.25),
  xlab=expression(beta[1]),
  ylab=expression(beta[2]))