demo/xrc2.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/xrc2.R
##
## This program shows all features of the package 'ipeglim' 
## in the context of (zero-truncated) Poisson regression model.
##
## * multivariate normal imprecise prior
## * no intercept sampling model
## * X is centered
## 
## * log(E(y)) ~ b1*x1 + b2*x2
# param <- c(-1, -1) # very good parameters
# param <- c(-1, -1.5) # very good parameters
# param <- c(-1.5, 0.5)


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

source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/cpef2reg.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/pbox_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/qc.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/extract_MHchain.R")

# structure of X
xm <- c(3,-2)
xsd <- c(1,1)
xr <- diag(2)
xr[2] <- xr[3] <- r <- 0
xsigma <- t(diag(xsd))%*%xr%*%diag(xsd)

# 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

SILENT <- TRUE
ZTRUNC <- TRUE

# random variates generation
tmp <- rvg4yx(N=3e2, pars=bcfs4, xreg=TRUE, ztrunc=ZTRUNC, kind="mvnorm", mean=xm, sigma=xsigma, center.X=TRUE)
dat <- tmp$yX
table(dat$y)
xi1 <- unlist(dat[1,-1])

m2fit <- model(formula=y~x1+x2-1, 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=SILENT)
sla <- summary(ola, HT.est=TRUE, silent=SILENT)
t1 <- proc.time()


if(FALSE){
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)
}