inst/thesis/fig/lg_linear_1.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

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

rm(list=ls())
options(warn=2)
library(ipeglim)

set.seed(18372342)
n <- 10
y <- rvg4yx(N=n, pars=2, ztrunc=FALSE, xreg=FALSE, kind="poisson")$y

y1 <- y[1]
y2 <- y[1:2]
y3 <- y[1:5]
y4 <- y[1:10]


lc1 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2))), rhs=c(0,0,-10,-10))

## with y1
m2fit1 <- model(formula=y1~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
cmfit1 <- iprior(obj=m2fit1, eqns=lc1, verbose=FALSE)
cmfit1 <- set.grid(obj=cmfit1, len=10)
sop1 <- summary(update(obj=cmfit1, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)


## with y2 
m2fit2 <- model(formula=y2~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
cmfit2 <- iprior(obj=m2fit2, eqns=lc1, verbose=FALSE)
cmfit2 <- set.grid(obj=cmfit2, len=10)
sop2 <- summary(update(obj=cmfit2, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

## with y3
m2fit3 <- model(formula=y3~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
cmfit3 <- iprior(obj=m2fit3, eqns=lc1, verbose=FALSE)
cmfit3 <- set.grid(obj=cmfit3, len=10)
sop3 <- summary(update(obj=cmfit3, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

## with y4
m2fit4 <- model(formula=y4~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
cmfit4 <- iprior(obj=m2fit4, eqns=lc1, verbose=FALSE)
cmfit4 <- set.grid(obj=cmfit4, len=10)
sop4 <- summary(update(obj=cmfit4, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)


par(mfrow=c(2,2), mar=c(4,4,0,0))
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=expression(n==1), 
     xlab=expression(alpha), ylab=expression(beta), 
     zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90), 
     split=c(1,1,2,2), more=TRUE)

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

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

plot(sop4, zlim=zlim, xlim=xlim, ylim=ylim, main=expression(n==10), 
     xlab=expression(alpha), ylab=expression(beta), 
     zlab=list(expression(E~"("~theta~"|"~y~")"), rot=90), 
     split=c(2,2,2,2), more=TRUE)