inst/thesis/fig/lg_Dn_numeric.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## lg_Dn.R

rm(list=ls())

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

tbni <- list()
idx <- 0
size <- c(5, 10, 20, 50, 100, 300, 500, 1000)
for(i in size){ 
idx <- idx + 1
y <- y0[1:i]

m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)

## scenario 1.
lc1 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2))), rhs=c(0,0,-10,-10))
cmfit <- iprior(obj=m2fit, eqns=lc1, verbose=FALSE)
cmfit <- set.grid(obj=cmfit, len=10)
sop1 <- summary(update(obj=cmfit, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)


## scenario 2. 
lc2 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2)), c(-1,1)), rhs=c(1,1,-8,-8,-4))
cmfit2 <- iprior(obj=m2fit, eqns=lc2)
sop2 <- summary(update(obj=cmfit2, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

## scenario 3.
lc3 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2)), c(1,-1)), rhs=c(1,1,-8,-8,-3))
cmfit3 <- iprior(obj=m2fit, eqns=lc3)
sop3 <- summary(update(obj=cmfit3, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)


## scenario 4. 
lc4 <- list(x=5,y=5,r=5,len=30)
cmfit4 <- iprior(obj=m2fit, circle=lc4)
sop4 <- summary(update(obj=cmfit4, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

sop <- list(sop1, sop2, sop3, sop4)
delta <- do.call(c, lapply(sop, "[[", "delta"))
inf <- do.call(c, lapply(sop, "[[", "inf"))
sup <- do.call(c, lapply(sop, "[[", "sup"))

tab <- as.data.frame(rbind(inf, sup, delta))
names(tab) <- c("S1", "S2", "S3", "S4")
tbni[[idx]] <- tab
}


names(tbni) <- paste("$n=", size, "$", sep="")
sc1 <- do.call(cbind, lapply(tbni, with, S1))
rownames(sc1) <- c("$\\Lprevision{\\theta}{}$", "$\\Uprevision{\\theta}{}$", "$\\Delta_{n}$")
print(xtable(sc1, caption="Scenario 1", digits=3), sanitize.text.function=function(x){x})

sc2 <- do.call(cbind, lapply(tbni, with, S2))
rownames(sc2) <- c("$\\Lprevision{\\theta}{}$", "$\\Uprevision{\\theta}{}$", "$\\Delta_{n}$")
print(xtable(sc2, caption="Scenario 2", digits=3), sanitize.text.function=function(x){x})

sc3 <- do.call(cbind, lapply(tbni, with, S3))
rownames(sc3) <- c("$\\Lprevision{\\theta}{}$", "$\\Uprevision{\\theta}{}$", "$\\Delta_{n}$")
print(xtable(sc3, caption="Scenario 3", digits=3), sanitize.text.function=function(x){x})

sc4 <- do.call(cbind, lapply(tbni, with, S4))
rownames(sc4) <- c("$\\Lprevision{\\theta}{}$", "$\\Uprevision{\\theta}{}$", "$\\Delta_{n}$")
print(xtable(sc4, caption="Scenario 4", digits=3), sanitize.text.function=function(x){x})