inst/thesis/fig/lg_pbox.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

# using markov chain
# pbxo

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

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

size <- c(10, 20, 50, 100)

# 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.ah <- 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)
	op <- update(obj=cmfit, method="AS", apriori="lgamma", silent=TRUE)
	sop <- summary(op, silent=TRUE)
	
	# op2 <- update(obj=cmfit, method="MH", apriori="lgamma", silent=TRUE, control=list(len.chain=2e3, len.burnin=1e3))
	# when the actual compilation needs to use the 5e4 chain 
	op2 <- update(obj=cmfit, method="MH", apriori="lgamma", silent=TRUE, control=list(len.chain=5e4, len.burnin=1e3))
	sop2 <- summary(op2, silent=TRUE)

	tab <- cbind(sop$est, sop2$est)
	colnames(tab) <- c("AS", "MH")
	
	if (count == 0) {
		xlim0 <- range(sop2$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.ah[[count]] <- sop
	sop.mh[[count]] <- sop2

	MHchain <- sop2$MHchain
	MHchain <- as.data.frame(lapply(MHchain, sort))
	mdif[count] <- mean(apply(MHchain,1,function(x) max(x)-min(x)))
}

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