inst/thesis/fig/lg_tsl.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/inst/thesis/fig/lg_tsl.R
##
## Chel Hee Lee 
## 2013.NOV.09
##
## Translation of hyperparameters on the hyperparameter space
## concerning the natural parameter
## Four figures

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

## Plot 1. Top left 
## complete ignorance (or vacuous) prior 
n <- xlim <- ylim <- 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")

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[1:2], function(x) polygon(x[chull(x),]))
text(cnt[1:2,], labels=paste("y",c(0,1),sep=""))

## Plot 3. Bottom left
## taking the first four samples
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[1:5], function(x) polygon(x[chull(x),]))
text(cnt[1:5,], labels=paste("y",seq(0,4),sep=""))


## Plot 4. Bottom right 
## taking the first nine samples
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[1:10], function(x) polygon(x[chull(x),]))
text(cnt[1:10,], labels=paste("y",seq(0,9),sep=""))