demo/xreg-shkg.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/xreg2d-cx0.R
## xreg2d-foo1.R
## 
## 2013.10.19 (Saturday)
## Template for the model including an intercept with centered covariates

library(ipeglim)
rm(list=ls())
# set.seed(1)
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/rvg4yx.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/bayesPoisXreg.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/iprior.imprecise.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/plot.summary.imprecise.R")

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

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

# random variates generation
tmp <- rvg4yx(N=1e3, cf0=0.5, pars=c(-0.5), ztrunc=TRUE, is.reg=TRUE, kind="mvnorm", mean=xm, sd=xsd, center.X=TRUE)
dat <- tmp$yX
table(dat$y)

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()
xi1 <- c(1, unlist(dat[1,-1]))
op <- update(obj=cmfit, method="LA", B=Bmat)
sop <- summary(op, HT.est=TRUE)
X11()
plot(sop)
print(sop$est)
print(proc.time()-t0)

if(FALSE){
t1 <- proc.time()
op1 <- update(obj=cmfit, method="MH", B=Bmat, xi=xi1)
sop1 <- summary(op1)
X11()
plot(sop1)
print(sop1$est)
print(proc.time()-t1)
}


lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(rep(-1,4)))
cmfit <- iprior(obj=m2fit, eqns=lc0)

# plotting extreme points on the xy-grid 
cmfit <- setGrid(obj=cmfit, len=5)

# imprecise model fitting 
xi1 <- unlist(c(1,dat[1,-1]))
op <- update(obj=cmfit, method="LA2", B=Bmat, xi=xi1)
sop <- summary(op, HT.est=TRUE)
sop$est
sop$accept.rate
plot(sop, # for m0xtms
  main="Imprecise Posterior Shrinkage", 
  sub=list(expression(-1<beta[0]~","~beta[1]<1), cex=1.25),
  xlab=expression(beta[0]),
  ylab=expression(beta[1]))