demo/lg-tsl.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())
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/print_imprecise.R")
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/plot_imprecise.R")


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

# random variates from Poisson
TRUNCATION <- FALSE
lam <- 1
n <- 1e2
y0 <- 
	if(!TRUNCATION){
		rpois(n=n, lambda=lam) 
	} else {
		rvg4yx(N=1e2, pars=lam, Xreg=FALSE, ztrunc=TRUNCATION, kind="pois")$y
	}
	
sz <- c(6, 8, 10, 12, 14)
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(0, -4, 0, -4))
xtms <- list()
i <- 0

for(s in sz){
  i <- i + 1
  y <- y0[seq_len(s)]
  if(i==1){
    fit <- summary(update(iprior(model(formula=y~0, ztrunc=TRUNCATION, dist="poisson"), eqns=lc0), method="LA", apriori="lgamma"))
    xtms[[i]] <- do.call(rbind, fit$xtms)
  }
  else fit <- summary(update(iprior(model(formula=y~0, ztrunc=TRUNCATION, dist="poisson"), eqns=lc0), method="LA", apriori="lgamma"))
  xtms[[i+1]] <- lc <- do.call(rbind, lapply(fit$m1, with, tslpars))
  polygon(x=xtms[[i]][chull(xtms[[i]]),]) # cvxh() is deprecated (2013.10.15)
}

centered <- do.call(rbind, lapply(xtms, colMeans))
text(x=centered, labels=paste("n=",c(0,sz),sep=""))
title(main="Translation Behaviour\n(Log-Gamma Imprecise Prior on Poisson model)")
text(x=45,y=45, labels=expression(mu==1))