inst/thesis/fig/o3ef_optim4.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## 2013.11.12
## with normal prior
## optimization is test

rm(list=ls())
library(ipeglim)

sfn <- function(h, nutau2=FALSE, ...){

	n <- length(y)
	sumy <- sum(y)
	if(nutau2) eta <- c(0.5/h[2], sumy+h[1]/h[2], n)
	else eta <- c(h[2], h[1]+sumy, n)

	g <- function(t) t*exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	vg <- integrate(g, lower=-Inf, upper=Inf)$value
	h <- function(t) exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	vh <- integrate(h, lower=-Inf, upper=Inf)$value

	gd2 <- function(t) -t^3*exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	vgd2 <- integrate(gd2, lower=-Inf, upper=Inf)$value
	hd2 <- function(t) -t^2*exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	vhd2 <- integrate(hd2, lower=-Inf, upper=Inf)$value
	ede2 <- (vgd2*vh-vg*vhd2)/vh^2
	
	gd1 <- function(t) t^2*exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	vgd1 <- integrate(gd1, lower=-Inf, upper=Inf)$value
	hd1 <- function(t) t*exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	vhd1 <- integrate(hd1, lower=-Inf, upper=Inf)$value
	ede1 <- (vgd1*vh-vg*vhd1)/vh^2

#	gd0 <- function(t) -t*exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t)+t)
#	vgd0 <- integrate(gd0, lower=-Inf, upper=Inf)$value
#	hd0 <- function(t) -exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t)+t)
#	vhd0 <- integrate(hd0, lower=-Inf, upper=Inf)$value
#	ede0 <-(vgd0*vh-vg*vhd0)/vh^2
	
	val <- c(ede2, ede1)
	return(val)
}


ofn <- function(h, nutau2=FALSE, ...){
	n <- length(y)
	sumy <- sum(y)
	if(nutau2) eta <- c(0.5/h[2], sumy+h[1]/h[2], n)
	else eta <- c(h[2], h[1]+sumy, n)

	f1 <- function(t) t*exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	f2 <- function(t) exp(-eta[1]*t^2+eta[2]*t-exp(log(eta[3])+t))
	v1 <- tryCatch(integrate(f1, lower=-Inf, upper=Inf)$value, error=function(e) return(NA))
	v2 <- tryCatch(integrate(f2, lower=-Inf, upper=Inf)$value, error=function(e) return(NA))
	# val <- suppressWarnings(log(v1)) - suppressWarnings(log(v2))
	val <- tryCatch(v1/v2, error=function(e) return(NA))
	return(val)
}

source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/cpef.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/xtms2eqns.R")

# natural specification by nu and tau2 
# nlc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-1.5,-1.5, 0.5,-4)) # nutau2=TRUE
nlc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-3, -3, 0.5, -2.5))  # nutau2=FALSE

nxtm <- iprior(eqns=nlc0)
nx <- as.data.frame(do.call(rbind, nxtm$xtms))
names(nx) <- c("nu", "tau2")

# identify extremes on the hyperparameter space which is mapped from natural specification
tmp <- list()
tmp$e2 <- 1/(2*nx$tau2)
tmp$e1 <- nx$nu/nx$tau2
tmp$xtms <- as.list(as.data.frame(rbind(tmp$e1, tmp$e2)))
names(tmp$xtms) <- paste("x", 1:4, sep="")

par(mfrow=c(2,2))
# find new inequalities for hyperparameter specification
x2e <- xtms2eqns(tmp, show=TRUE, type="eqns")
# its order is now hparam=(e1,e2);
hlc0 <- list(lhs=x2e$ui, rhs=x2e$ci)
hxtm <- iprior(eqns=hlc0)
plot(hxtm)
# -3 < e1 < 3, 0 < e2 < 1

y0 <- rpois(n=1, lambda=2)
y <- y0
m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson")
# cmfit <- iprior(obj=m2fit, eqns=hlc0)
cmfit <- iprior(obj=m2fit, eqns=nlc0)
cmfit1 <- set.grid(cmfit, 20)
op <- update(obj=cmfit1, apriori="normal", method="LA", silent=TRUE, nutau2=FALSE)
sop <- summary(op, silent=TRUE)
# plot(sop, bnd=tmp, xlim=c(-5,5), ylim=c(0,1))
# plot(sop, xlab=expression(nu), ylab=expression(tau2)) ## nutau2=TRUE
plot(sop, xlab=expression(eta[1]), ylab=expression(eta[2])) ## nutau2=FALSE
# sop$tslpars

mini <- constrOptim(theta=x2e$cc, f=ofn, grad=sfn, nutau2=FALSE, ui=x2e$ui, ci=x2e$ci)
maxi <- constrOptim(theta=x2e$cc, f=ofn, grad=sfn, nutau2=FALSE, ui=x2e$ui, ci=x2e$ci, control=list(fnscale=-1))
mini$par
maxi$par
sop$xtms[sop$inf.idx]
sop$xtms[sop$sup.idx]




if(FALSE){ 
## investigate the change of the shape
## nu = 1
## 1/2 <= tau2 <= 2

## eta2 = (1/2*tau2)
## eta1 = 2*eta2 + sumy

## once tau2 is given, 1 <= tau2 <= 4
nu <- seq(0.5, 1.5,0.1)
tau2 <- seq(0.5,1.5,0.1)
# e2 <- 1/(2*tau2)
# e1 <- 2*e2*nu + sumy

nsp <- expand.grid(nu=nu, tau2=tau2)
# hsp <- expand.grid(e2=e2, e1=e1)
hsp1 <- list()
hsp1$e2 <-1/(2*nsp$tau2)
hsp1$e1 <- nsp$nu/nsp$tau2
hsp1 <- as.data.frame(hsp1)
# plot(hsp1)
# plot(nsp)
# plot(hsp)
}