inst/thesis/fig/o3ef_optima.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

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

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

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


ofn <- function(hparam, ...){
	n <- length(y)
	sumy <- sum(y)
	eta <- c(hparam[2], hparam[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)
}

sfn <- function(hparam, ...){

	n <- length(y)
	sumy <- sum(y)
	eta <- c(hparam[2], hparam[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(ede1, ede2)
	return(val)
}

# default
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-1, -1, 0.1, -1))
# lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-3, -3, 0.1, -2.5))
set.seed(16979238)
y0 <- rpois(n=1, lambda=2)
y <- y0

mini <- constrOptim(theta=c(0.5,0.5), f=ofn, grad=sfn, ui=lc0$lhs, ci=lc0$rhs)
maxi <- constrOptim(theta=c(0.5,0.5), f=ofn, grad=sfn, ui=lc0$lhs, ci=lc0$rhs, control=list(fnscale=-1))

m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson")
cmfit <- iprior(obj=m2fit, eqns=lc0)
cmfit <- set.grid(obj=cmfit, len=15)
op <- update(obj=cmfit, method="LA", apriori="normal", silent=TRUE, control=list(const=30))
sop <- summary(op, silent=TRUE)

tb0min <- c(sop$m1[[sop$inf.idx]]$hparam, sop$m1[[sop$inf.idx]]$theta)
tb1min <- c(mini$par, mini$val)
tb0max <- c(sop$m1[[sop$sup.idx]]$hparam, sop$m1[[sop$sup.idx]]$theta)
tb1max <- c(maxi$par, maxi$val)
tab <- cbind(tb0min, tb1min, tb0max, tb1max)
round(tab,3)