inst/thesis/fig/lg_tsleta2.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

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

library(ipeglim)
rm(list=ls())
options(warn=2)
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/iprior.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/s4xtms.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/model.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/print_imprecise.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/plot_imprecise.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/summary_imprecise.R")


# background
a <- b <- seq(1,2,0.1)
plot(a,b,xlim=c(0,10), ylim=c(0,10), 
  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/0.75, 1/0.5, 1/0.25)
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, 0.75, 0.5, 0.25)*9.5
y <- c(0.25, 0.5, 0.75, 1, 1, 1, 1)*9.5
text(x=x, y=y, labels=round(rev(m0),2))

# random variates from Poisson
set.seed(16979238)
lam <- 1
n <- 1e1
y0 <- rvg4yx(N=1e2, pars=lam, xreg=FALSE, ztrunc=FALSE, kind="poisson")$y
	
sz <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
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 <- numeric(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]
}

lapply(xtms, function(x) polygon(x[chull(x),]))
centered <- do.call(rbind, lapply(xtms, colMeans))
text(x=centered, labels=paste("n=", c(0,sz), sep=""))
# title(main="Log-Gamma Imprecise Prior")
text(x=9,y=9, labels=expression(mu==1), pos=4, cex=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))