inst/thesis/fig/lg_tsleta11.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## ipeglim/demo/lgamma-translation.R
##
## Demonstration of translation behaviour (log-gamma imprecise prior)


rm(list=ls())

set.seed(16979238)
# background
lim <- 10
a <- seq(1e-3, lim, by=0.1)
b <- seq(0.1, lim, by=0.1)
z <- outer(X=a, Y=b, FUN=function(X,Y) { 
	v <- digamma(X) - log(Y); 
	v
	# exp(v) # for mu 
	} )
contour(x=a, y=b, z=z, 
				# levels=c(0.5, 1, 2, 4), # for mu
			  levels=seq(-3, 3, 0.5), 
        xlim=c(0,lim), ylim=c(0,lim), 
        method="edge", 
        labcex=1
        )

y0 <- rvg4yx(N=lim, pars=1, xreg=FALSE, ztrunc=FALSE, kind="poisson")$y
	
sz <- seq(1, lim, 1)
xtms <- list()
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(0, -1, 0, -1))

xtms0 <- iprior(eqns=lc0)
xtms[[1]] <- do.call(rbind, xtms0$xtms)
ym <- list(length(sz))

for(i in seq_len(length(sz))){
  y <- y0[seq_len(sz[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]
	ym[[i]] <- c(sop$inf, sop$sup)
}

lapply(xtms, function(x) polygon(x[chull(x),]))
# lapply(xtms, function(x) abline(h=x[1,]))
centered <- do.call(rbind, lapply(xtms, colMeans))
text(x=centered, labels=paste("n=", c(0,sz), sep=""))

## with lc1
sz1 <- sz
lc1 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(3, -4, 0, -1))

xtms1 <- list()
xtms00 <- iprior(eqns=lc1)
xtms1[[1]] <- do.call(rbind, xtms00$xtms)
ym1 <- list(length(sz))

for(i in seq_len(length(sz1))){
  y1 <- y0[seq_len(sz1[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]
	ym1[[i]] <- c(sop1$inf, sop1$sup)
}
# print(xtms1)
lapply(xtms1, function(x) polygon(x[chull(x),], lwd=2, border="red"))
centered1 <- do.call(rbind, lapply(xtms1, colMeans))

for(i in 1:nrow(centered1)) {
	arrows(x0=centered[i,1], y0=centered[i,2], x1=centered1[i,1], y1=centered1[i,2], length=0.1)
}

t0 <- do.call(rbind, ym)
t1 <- do.call(rbind, ym1)
dev.off()