inst/thesis/fig/lg_tsleta.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())

# background
a <- b <- seq(1,2,0.1)
bound <- 1e2
plot(a,b,xlim=c(0,bound), ylim=c(0,bound), 
 xlab=expression(alpha~"(shape)"), 
 ylab=expression(beta~"(rate)"),
 type="n")

# lines for expectations on mu
m0 <- c(0.25, 0.5, 0.75, 1, 1.5, 2.5, 4, 8)
# for(i in seq_len(length(m0))) abline(a=0, b=1/m0[i], lty="dashed")
# abline(a=0, b=1)
# abline(v=0)
# abline(h=0)

x <- c(1, 1, 1, 1, 1, 0.75, 0.5, 0.25)*bound
y <- c(1/8, 1/4, 1/2.5, 1/1.5, 1, 1, 1, 1)*bound

# text(x=x, y=y, labels=rev(m0))

# random variates from Poisson
# set.seed(16979238)
lam <- 1
y0 <- rvg4yx(N=bound, pars=lam, xreg=FALSE, ztrunc=FALSE, kind="poisson")$y
	
# sz <- c(5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 150, 200, 250, 300, 350, 400)
# sz <- seq(5, 900, 50)
sz <- seq(1, bound, 1)
xtms <- list()
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(0, -1, 0, -5))

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),], lwd=3))
# lapply(xtms, function(x) abline(h=x[1,]))
# centered <- do.call(rbind, lapply(xtms, colMeans))

## with lc1
sz1 <- sz
lc1 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(10, -15, 0, -5))

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=3, border="red"))
# centered1 <- do.call(rbind, lapply(xtms1, colMeans))


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






# text(x=centered, labels=paste("n=", c(0,sz), sep=""))
# title(main="Log-Gamma Imprecise Prior")
# text(x=45,y=45, labels=expression(mu==1))

# points(centered)
# it <- rowSums(centered)
# abline(a=intercept, b=rep(-1, length(intercept)))
# abline(a=it[1], b=-1)
# abline(a=it[2], b=-1)
# pp <- cbind(it/2, it/2)
# points(pp)

# edist <- sqrt(rowSums((centered-pp)^2))
# k <- diff(xtms[[1]][1,])*sqrt(2)
# print(mean(edist<k))