inst/thesis/fig/lg_pdconflict.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(2,3,2,1)+0.1)
set.seed(16979238)

## Plot 1. Top left 
## complete ignorance (or vacuous) prior 
xlim <- ylim <- 10
n <- 10
a <- seq(1e-3, xlim, 0.1)
b <- seq(0.1, ylim, 0.1)
th0 <- outer(X=a, Y=b, FUN=function(X,Y) digamma(X) - log(Y))
contour(x=a, y=b, z=th0, levels=seq(-3, 3, 0.5), xlim=c(0, xlim), ylim=c(0, ylim), method="edge", labcex=1, lty="dashed", xlab=expression(alpha), ylab=expression(beta))

th1 <- xtms <- list()
lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(0,0,-1,-1))
R0 <- iprior(eqns=lc0)
xtms[[1]] <- R0 <- as.matrix(do.call(rbind, R0$xtms))
polygon(R0[chull(R0), ], col=alpha("azure2",0.75), border="darkblue")
axis(side=1, col="darkred")
abline(v=0, h=0, col="darkred", lwd=1.5)

## uniform (alpha=1, rate=0)
points(x=1, y=0, cex=1.5, pch=16, col="red")
# text(1,0, label="Uniform")
## Jeffery (alpha=1/2, rate=0)
points(x=1/2, y=0, cex=1.5, pch=17, col="blue")
# text(1,1/2, label="Jeffery")
legend("bottomright", legend=c("Uniform", "Jeffery"), col=c("red", "blue"), pch=c(16, 17))


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, sop$delta)
}
cnt <- do.call(rbind, lapply(xtms, colMeans))

## Plot 2. Top right
## taking 1 sample
contour(x=a, y=b, z=th0, levels=seq(-3, 3, 0.5), xlim=c(0, xlim), ylim=c(0, ylim), method="edge", labcex=1, lty="dashed", xlab=expression(alpha), ylab=expression(beta))
lapply(xtms, function(x) polygon(x[chull(x),], col=alpha("azure2", 0.5)))
text(cnt, labels=paste("y",0:8,sep=""))


th2 <- xtms1 <- list()
lc1 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(3,0,-4,-1))
R1 <- iprior(eqns=lc1)
xtms1[[1]] <- R1 <- as.matrix(do.call(rbind, R1$xtms))
polygon(R0[chull(R0), ], col=alpha("azure2",0.75), border="darkblue")
axis(side=1, col="darkred")
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, sop1$delta)
}
cnt1 <- do.call(rbind, lapply(xtms1, colMeans))
lapply(xtms1, function(x) polygon(x[chull(x),], border="darkred"))
text(cnt1, labels=paste("y",0:8,sep=""))

# points(cnt)
it <- rowSums(cnt)
for(i in 1:nrow(cnt)) {
	arrows(x0=cnt[i,1], y0=cnt[i,2], x1=cnt1[i,1], y1=cnt1[i,2], length=0.1)
}

th1 <- do.call(rbind, th1)
th2 <- do.call(rbind, th2)