Nothing
## 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]))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.