inst/thesis/fig/lg_pdconflict2.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/inst/thesis/fig/lg_pdconflict.R
##
## prior-data conflict
## 2013.NOV.08

rm(list=ls())
options(warn=2)
library(ipeglim)
par(mfrow=c(1,2), mar=c(4,4,0,0))
set.seed(16979238)

lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(0,0,-1,-1))
lc1 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(3, 0,-4,-1))
lc2 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(0, 5,-1,-6))

## natural imprecise prior
th1 <- xtms <- list()
R0 <- iprior(eqns=lc0)
xtms[[1]] <- R0 <- as.matrix(do.call(rbind, R0$xtms))

n <- 2e2
y0 <- rvg4yx(N=n, pars=(mu <- 1), xreg=FALSE, ztrunc=FALSE, kind="poisson")$y

for(i in 1:n){
	y <- y0[1:i]
	m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
	cmfit <- iprior(obj=m2fit, eqns=lc0, silent=TRUE)
	op <- update(obj=cmfit, method="AS", apriori="lgamma", silent=TRUE)
	sop <- summary(op, silent=TRUE)
	xtms[[i+1]] <- sop$tslpars[,-1]
	th1[[i]] <- c(sop$inf, sop$sup)
}
cnt <- do.call(rbind, lapply(xtms, colMeans))

## Plot 2. change on alpha 
th2 <- xtms1 <- list()
R1 <- iprior(eqns=lc1)
xtms1[[1]] <- R1 <- as.matrix(do.call(rbind, R1$xtms))

for(i in 1:n){
	y1 <- y0[1:i]
	m2fit1 <- model(formula=y1~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
	cmfit1 <- iprior(obj=m2fit1, eqns=lc1, silent=TRUE)
	op1 <- update(obj=cmfit1, method="AS", apriori="lgamma", silent=TRUE)
	sop1 <- summary(op1, silent=TRUE)
	xtms1[[i+1]] <- sop1$tslpars[,-1]
	th2[[i]] <- c(sop1$inf, sop1$sup)
}
cnt1 <- do.call(rbind, lapply(xtms1, colMeans))

## Plot 2. change on beta
th3 <- xtms2 <- list()
R1 <- iprior(eqns=lc2)
xtms2[[1]] <- R1 <- as.matrix(do.call(rbind, R1$xtms))

for(i in 1:n){
	y2 <- y0[1:i]
	m2fit2 <- model(formula=y2~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)
	cmfit2 <- iprior(obj=m2fit2, eqns=lc2, silent=TRUE)
	op2 <- update(obj=cmfit2, method="AS", apriori="lgamma", silent=TRUE)
	sop2 <- summary(op2, silent=TRUE)
	xtms2[[i+1]] <- sop2$tslpars[,-1]
	th3[[i]] <- c(sop2$inf, sop2$sup)
}
cnt2 <- do.call(rbind, lapply(xtms2, colMeans))

## plotting
th1 <- do.call(rbind, th1)
th2 <- do.call(rbind, th2) # change on alpha
th3 <- do.call(rbind, th3) # change on beta

plot(th1[,1], type="n", ylim=c(-0.5, 0.5), xlab="The size of a sample", ylab=expression(theta))
lines(th1[,1], type="l", col="black", lty="dashed", lwd=1.5) # lower (Natural)
lines(th1[,2], type="l", col="black", lwd=1.5)
lines(th2[,1], type="l", col="darkblue", lty="dashed", lwd=1.5) # lower (alpha)
lines(th2[,2], type="l", col="darkblue", lwd=1.5)
abline(h=0)
# abline(h=c(-0.05, 0.05), lty="dashed")
legend("topright", col=c("darkblue", "darkblue", "black", "black"), lty=c("solid", "dashed", "solid", "dashed"), legend=c("MaxPE from P1", "MinPE from P1", "MaxPE from P0", "MinPE from P0"), lwd=2)


plot(th1[,1], type="n", ylim=c(-0.5, 0.5), xlab="The size of a sample", ylab=expression(theta))
lines(th1[,1], type="l", col="black", lty="dashed", lwd=1.5) # lower (Natural)
lines(th1[,2], type="l", col="black", lwd=1.5)
lines(th3[,1], type="l", col="darkred", lty="dashed", lwd=1.5) # lower (beta)
lines(th3[,2], type="l", col="darkred", lwd=1.5)
abline(h=0)
# abline(h=c(-0.05, 0.05), lty="dashed")

legend("topright", col=c("darkred", "darkred", "black", "black"), lty=c("solid", "dashed", "solid", "dashed"), legend=c("MaxPE from P1", "MinPE from P1", "MaxPE from P0", "MinPE from P0"), lwd=2)