inst/thesis/fig/o3ef_pbox.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/inst/thesis/fig/o3ef_pbox.R
## 2013.11.13

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

# default
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-3, -3, 0.1, -2.5))
# lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-3, -3, 0.1, -2.5))
# set.seed(16979238)
set.seed(18372342)
y0 <- rvg4yx(N=200, pars=2, ztrunc=FALSE, xreg=FALSE, kind="poisson")$y

# size <- c(15, 20, 50, 100)
size <- c(20, 40, 60)

# par(mfrow=c(2,2), mar=c(2,3,2,1)+0.1)
lc0 <- 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)

count <- 0
mdif <- numeric(length(size))
sop.mh <- list()
sop.aq <- list()
sop.la <- list()
tab <- list()
for (i in size) { 
	y <- y0[1:i]
	
	m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
	cmfit <- iprior(obj=m2fit, eqns=lc0, verbose=FALSE)
	
	saq <- summary(update(obj=cmfit, method="AQ", apriori="normal", silent=TRUE), silent=TRUE)
	sla <- summary(update(obj=cmfit, method="LA", apriori="normal", silent=TRUE, control=list(const=3e2)), silent=TRUE)
	smh <- summary(update(obj=cmfit, method="MH", apriori="normal", silent=TRUE, control=list(len.chain=3e4, len.burnin=2e3)), silent=TRUE)

	if (count == 0) {
		xlim0 <- range(smh$MHchain)
	}
	# print(pbox(x=sop2, pretty=TRUE, control=list(xtms=c(1:6), xlim=xlim0), main=paste("n=", i)), split=c((count %/% 2 + 1), (count %% 2 + 1), 2, 2), more=TRUE)
#	pbox(x=sop2, pretty=TRUE, control=list(xtms=c(1:6), xlim=xlim0), main=paste("n=", i))
#	abline(v=log(2), lty="dashed")
#	text(x=log(2), y=0.1, labels=expression(theta==log~"("~2~")"))
	count <- count + 1
	sop.aq[[count]] <- saq
	sop.mh[[count]] <- smh
	sop.la[[count]] <- sla

#	MHchain <- smh$MHchain
#	MHchain <- as.data.frame(lapply(MHchain, sort))
#	mdif[count] <- mean(apply(MHchain,1,function(x) max(x)-min(x)))
	tmp.tab <- cbind(saq$est, sla$est, smh$est)
	tmp.tab <- rbind(tmp.tab, c(saq$delta, sla$delta, smh$delta))
	colnames(tmp.tab) <- c("AQ", "LA", "MH")
	tab[[count]] <- tmp.tab
}

save(list=ls(), file="~/Documents/PhD/ipeglim/thesis/fig/o3ef_pbox.RData")