demo/xreg-a-single.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## this is the test file
## 2013.10.23

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

## 
xm <- c(1,0)
xsd <- c(0,1)
xr <- diag(2)
xr[2] <- xr[3] <- r <- 0.0
xsigma <- t(diag(xsd))%*%xr%*%diag(xsd)
tmp <- rvg4yx(N=1e3, pars=c(-0.5, 1.5), ztrunc=TRUE, xreg=TRUE, kind="mvnorm", mean=xm, sigma=xsigma, center.X=FALSE)
dat <- tmp$yX
print(table(dat$y))
y <- dat$y

fit0s <- summary(fit0)
fit0s

# ipeglim with a single prior
m2fit <- model(formula=y~x1+x2-1, data=dat, ztrunc=TRUE, dist="poisson")
cmfit <- impose(obj=m2fit, x=fit0$cfs)
op <- update(obj=cmfit, method="LA", B=diag(2)*1e-2)
sop <- summary(op)
sop

# ipeglim with a set of priors (LA or LA2)
m2fit <- model(formula=y~x1+x2-1, data=dat, ztrunc=TRUE, dist="poisson")
lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(rep(-1,4)))
cmfit <- impose(obj=m2fit, eqns=lc0)
op <- update(obj=cmfit, method="LA2", B=diag(2)*1e-2, const=1e5)
sop <- summary(op)
sop

# ipeglim with a set of priors (MH)
m2fit <- model(formula=y~x1+x2-1, data=dat, ztrunc=TRUE, dist="poisson")
lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(rep(-1,4)))
cmfit <- impose(obj=m2fit, eqns=lc0)
op <- update(obj=cmfit, method="LA2", B=diag(2)*1e-2, const=1e5)
sop <- summary(op)
sop