inst/thesis/fig/lg_Dn.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## 2013-NOV-09
## lg_Dn.R

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

set.seed(18372342)

y <- rvg4yx(N=1, pars=2, ztrunc=FALSE, xreg=FALSE, kind="poisson")$y
m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)

# scenarios
lc1 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2))), rhs=c(0,0,-10,-10))
lc2 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2)), c(-1,1)), rhs=c(1,1,-8,-8,-4))
lc3 <- list(lhs=rbind(c(1,1), c(-1,-1), c(-1,1), c(1,-1), c(0,-1), c(0,1)), rhs=c(3,-7,-3,-1,-3,1)*2)
lc4 <- list(x=5,y=5,r=5,len=25)

cmfit1 <- iprior(obj=m2fit, eqns=lc1, verbose=FALSE)
cmfit1 <- set.grid(obj=cmfit1, len=10)
cmfit2 <- iprior(eqns=lc2, verbose=FALSE)
cmfit3 <- iprior(eqns=lc3, verbose=FALSE)
cmfit4 <- iprior(circle=lc4, verboser=FALSE)

sop1 <- summary(update(obj=cmfit1, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

# par(mfrow=c(2,2), mar=c(2,3,2,1)+0.1)
z <- range(sop1$est)
z.mar <- 0.2
zlim <- c(z[1]-z.mar, z[2]+z.mar)
xlim <- ylim <- c(-1,11)
plot(sop1, zlim=zlim, xlim=xlim, ylim=ylim, main="Region 1", 
     xlab=expression(alpha), ylab=expression(beta), 
     zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90), 
     split=c(1,1,2,2), more=TRUE)

plot(sop1, bnd=cmfit2, zlim=zlim, xlim=xlim, ylim=ylim, main="Region 2", 
     xlab=expression(alpha), ylab=expression(beta), 
     zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90), 
     split=c(2,1,2,2), more=TRUE)

plot(sop1, bnd=cmfit3, zlim=zlim, xlim=xlim, ylim=ylim, main="Region 3", 
     xlab=expression(alpha), ylab=expression(beta), 
     zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90), 
     split=c(1,2,2,2), more=TRUE)

plot(sop1, bnd=cmfit4, zlim=zlim, xlim=xlim, ylim=ylim, main="Region 4", 
     xlab=expression(alpha), ylab=expression(beta), 
     zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90), 
     split=c(2,2,2,2), more=TRUE)