demo/xrc2i.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/xreg-circle2d.R
##
## Poisson regression model with a multivariate normal imprecise prior
## * X is centered;
## * log(E(y)) ~ b0 + b1*x1
## * log(E(y)) ~ b1*x1 + b2*x2
## * log(E(y)) ~ b0 + b1*x1 + b2*x2
## * log(E(y)) ~ b1*x1 + b2*x2 + b3*x3
## tested by methods = 'LA', 'MH', 'IS'.
## 

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

if(FALSE){
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/update_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/plot_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/qc.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/rvg4yx.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/extract_MHchain.R")
}

# seed trial 1-20
# set.seed(4)  
# - 'LA' fails due to various reason (error handled 2013.10.27)
# - 
# set.seed(5)  
# - 'LA' gives bad points (when n < 5e2 in general)
# 
# set.seed(6)

## 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

# random variates generation
tmp <- rvg4yx(N=1e3, cf0=bcfs4[1], pars=bcfs4[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]))

source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/model.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/print_imprecise.R")
m2fit <- model(formula=y~x1, data=dat, ztrunc=TRUE, dist="poisson")
names(m2fit)
print(m2fit)

# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/iprior.R")

cmfit <- iprior(obj=m2fit, circle=list(x=0, y=0, len=10), verbose=TRUE)
names(cmfit)
print(cmfit)


op <- update(obj=cmfit, method="LA", B=Bmat, xi=xi1, silent=FALSE)
names(op)
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/model.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/print_imprecise.R")
print(op)

if(FALSE){
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/summary_imprecise.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/print_imprecise.R")
sop <- summary(op, HT.est=TRUE, silent=FALSE)
names(sop)

sop

plot(sop)

#source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/pbox_imprecise.R")
pbox(x=sop, which=list(beta=1, xtms=c(5,2)), pretty=TRUE, verbose=TRUE)
pbox(x=sop, which=list(beta=2), pretty=TRUE, verbose=TRUE)

omh <- update(obj=cmfit, method="MH", B=Bmat, xi=xi1, silent=SILENT)
smh <- summary(omh, HT.est=TRUE, silent=SILENT)
t2 <- proc.time()

ois <- update(obj=cmfit, method="IS", B=Bmat, xi=xi1, silent=SILENT)
sis <- summary(ois, HT.est=TRUE, silent=SILENT)
t3 <- proc.time()

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

tb.x1 <- cbind(sla$est$x1, smh$est$x1, sis$est$x1)
tb.x2 <- cbind(sla$est$x2, smh$est$x2, sis$est$x2)
print(tb.x1)
print(tb.x2)

X11()
plot(sla, # 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(smh, # 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(sis, # 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]))

print(ola)
print(omh)
print(ois)

print(sla)
print(smh)
print(sis)
}