demo/normal-test.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## 2013.10.21. Successful implementation for the zero-truncated Poisson case.
## add it to the package
## 'eta1fn'
## 2013.11.12 it will be removed soon 
## complete by 2013.11.13


library(ipeglim)
rm(list=ls())
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/cpef2reg.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/model.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/iprior.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/update_imprecise.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/summary_imprecise.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/plot_summary.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/pbox_imprecise.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/qc.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/rvg4yx.R")

fn.la <- function(eta0, eta1, eta2, const=5e2, tol=2e-3, start=rnorm(1), logm, ztrunc=TRUE, ...){

  fn <- function(theta, numer=FALSE, ...){
    val <- -eta2*theta^2 + eta1*theta - eta0*exp(theta)
    val <- if(numer) suppressWarnings(log(theta+const)) + val else val
    val <- if(ztrunc) val - eta0*suppressWarnings(ppois(q=0, lambda=exp(theta), lower.tail=FALSE, log.p=TRUE)) else val
    return(val)
  }
  
  gr <- function(theta, numer=FALSE, ...){
    val <- -2*eta2*theta + eta1 - eta0*exp(theta)
    val <- if(numer) 1/(theta+const) + val else val
    val <- if(ztrunc) val - eta0*exp(theta-exp(theta))/ppois(q=0, lambda=exp(theta), lower.tail=FALSE, log.p=FALSE) else val
    return(val)
  }
  
  fit0 <- optim(par=start, fn=fn, gr=gr, method="BFGS", control=list(fnscale=-1), hessian=TRUE)    
  fit1 <- optim(par=start, fn=fn, gr=gr, numer=TRUE, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
  
  sigma0 <- if(fit0$convergence == 0) -1/fit0$hessian else NaN
  sigma1 <- if(fit1$convergence == 0) -1/fit1$hessian else NaN
  
  sratio <- sigma1/sigma0
  fdiff <- fit1$value - fit0$value
  theta <- suppressWarnings(sqrt(sigma1/sigma0))*exp(fit1$value - fit0$value) - const 
  val <- theta - logm
  return(val) 
}

e0 <- seq(0,2,0.25)
e2 <- seq(0,5,0.3)
lmu <- c(-3, 2, 1, 0, -1, -2, 3)
hcube <- expand.grid(e0=e0, e2=e2, lmu=lmu)

e1 <- apply(hcube, 1, function(x){
  x <- as.vector(x)
  rt.e1 <- tryCatch(uniroot(fn.la, lower=-20, upper=30, tol=1e-4, eta0=x[1], eta2=x[2], logm=x[3])$root, error=function(e) return(NaN))
  return(rt.e1)
})



e1 <- apply(hcube, 1, function(x){
  x <- as.vector(x)
  source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/cpef.R")
  rt.e1 <- tryCatch(uniroot(cpef, lower=-20, upper=30, tol=1e-4, apriori="normal", method="LA", ztrunc=TRUE, expected.eta1=list(est=TRUE, eta0=x[1], eta2=x[2], logm=x[3]), initrun=TRUE, start=rnorm(1), verbose=TRUE)$root, error=function(e) return(NaN))
  return(rt.e1)
  })
  
  
hcube$e1 <- e1


# n=0
cube0 <- subset(hcube, subset=(e0==0))
with(cube0, plot(e1, e2, type='l', xlim=c(-20,20),
  main="Translation Behaviour\n(Normal Imprecise Prior on the Zero-Truncated Poisson Model)",
  xlab=expression(eta[1]),
  ylab=expression(eta[2]), col="black", lty=1
))
box0 <- rbind(c(-2,0.2), c(2,0.2), c(2,1.0), c(-2,1.0))
polygon(box0, border="black", lty=1, lwd=3)
mlab <- subset(cube0, subset=(e2=="2.7"))
with(mlab, text(x=e1, y=e2, labels=lmu)) 

# n=1
cube1 <- subset(hcube, subset=(e0==1))
with(cube1, lines(e1, e2, col="blue", lty=2))
y1 <- rpois(n=1, lambda=1)
box1 <- box0 + matrix(rep(c(y1,0),nrow(box0)), ncol=2, byrow=TRUE)
polygon(box1, border="blue", lty=2, lwd=3)

# n=2
cube2 <- subset(hcube, subset=(e0==2))
with(cube2, lines(e1,e2,col="red", lty=3))
y2 <- c(y1, rpois(n=1, lambda=1))
box2 <- box1 + matrix(rep(c(sum(y2),0), nrow(box1)), ncol=2, byrow=TRUE)
polygon(box2, border="red", lty=3, lwd=3)

legend(x=-20, y=1, 
  legend=c("n=0", "n=1", "n=2"),
  col=c("black", "blue", "red"),
  lty=c(1,2,3))