inst/thesis/fig/o3ef_optim3.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

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

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

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

	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], 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=TRUE, ...){
	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], 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)
}
# v <- sfn(eta=c(0.5,3,5))
# print(v)
# o <- ofn(eta=c(0.5, 10,2))
# print(o)
## what if 'constrOptim()' is used?
# }

# set.seed(16979238)
# lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-3, -3, 0.5, -2.5))
# lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-0.2, -0.5, 0.4, -1))  # Good
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-0.5, -0.5, 0.4, -1)) 


y0 <- rpois(n=10, lambda=2)
y <- y0
m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson")
cmfit <- iprior(obj=m2fit, eqns=lc0)
cmfit <- set.grid(cmfit, 10)
op <- update(obj=cmfit, apriori="normal", method="LA", silent=TRUE, nutau2=TRUE)
sop <- summary(op, silent=TRUE)

# rv1 <- optim(par=c(1,1), fn=ofn, gr=sfn, nutau2=TRUE, method="L-BFGS-B", lower=c(-3,0.5), upper=c(3,2.5)) 
# rv2 <- optim(par=c(1,1), fn=ofn, gr=sfn, nutau2=TRUE, method="L-BFGS-B", lower=c(-3,0.5), upper=c(3,2.5), control=list(fnscale=-1)) 

# rv1 <- optim(par=c(1,1), fn=ofn, gr=sfn, nutau2=FALSE, method="L-BFGS-B", lower=c(-3,0.5), upper=c(3,2.5)) 
# rv2 <- optim(par=c(1,1), fn=ofn, gr=sfn, nutau2=FALSE, method="L-BFGS-B", lower=c(-3,0.5), upper=c(3,2.5), control=list(fnscale=-1)) 


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

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


# cmfit <- iprior(eqns=lc0)
# ntbox <- as.data.frame(do.call(rbind, cmfit$xtms)) ## region on the intuitive space
# names(ntbox) <- c("nu", "tau2")
# xtms0 <- with(ntbox, cbind(0.5/tau2, nu/tau2)) ## region on the hyperparameter space
# xtms0 <- xtms0[chull(xtms0),]
# > xtms0
#     [,1] [,2]
#[1,]  1.0 -6.0
#[2,]  0.2 -1.2
#[3,]  0.2  1.2
#[4,]  1.0  6.0
# 
# 
if(FALSE){
sumy <- sum(y)
n <- length(y)
h1 <- hparam[1]
h2 <- hparam[2]

eta <- c(0.5/h2, sumy+h1/h2, n)


source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/xtms2eqns.R")
xtms1 <- sop$
mm <- list()
mm$xtms <- as.list(as.data.frame(t(xtms0)))
x2e <- xtms2eqns(obj=mm, show=TRUE)
x2e

mini <- constrOptim(theta=x2e$cc, f=fn, grad=gfn, ui=x2e$lhs, ci=x2e$rhs)
maxi <- constrOptim(theta=x2e$cc, f=fn, grad=gfn, ui=x2e$lhs, ci=x2e$rhs, control=list(fnscale=-1))
}


if(FALSE){
sfn <- function(eta, ...){

	fn11d2 <- function(t) -t^3*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v11d2 <- integrate(fn11d2, lower=-Inf, upper=Inf)$value
	
	fn10d2 <- function(t) t*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v10d2 <- integrate(fn10d2, lower=-Inf, upper=Inf)$value
	
	fn21d2 <- function(t) -t^2*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v21d2 <- integrate(fn10d2, lower=-Inf, upper=Inf)$value
	
	fn20d2 <- function(t) exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v20d2 <- integrate(fn10d2, lower=-Inf, upper=Inf)$value
	
	ede2 <- v11d2/v10d2 - v21d2/v20d2



	fn11d1 <- function(t) t^2*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v11d1 <- integrate(fn11d1, lower=-Inf, upper=Inf)$value
	
	fn10d1 <- function(t) t*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v10d1 <- integrate(fn10d1, lower=-Inf, upper=Inf)$value
	
	fn21d1 <- function(t) t*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v21d1 <- integrate(fn10d1, lower=-Inf, upper=Inf)$value
	
	fn20d1 <- function(t) exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v20d1 <- integrate(fn10d1, lower=-Inf, upper=Inf)$value
	
	ede1 <- v11d1/v10d1 - v21d1/v20d1



	fn11d0 <- function(t) -t*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t)+t)
	v11d0 <- integrate(fn11d0, lower=-Inf, upper=Inf)$value
	
	fn10d0 <- function(t) t*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v10d0 <- integrate(fn10d0, lower=-Inf, upper=Inf)$value
	
	fn21d0 <- function(t) -exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t)+t)
	v21d0 <- integrate(fn10d0, lower=-Inf, upper=Inf)$value
	
	fn20d0 <- function(t) exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v20d0 <- integrate(fn10d0, lower=-Inf, upper=Inf)$value
	
	ede0 <- v11d0/v10d0 - v21d0/v20d0
	
	val <- c(ede2, ede1, ede0)
	return(val)
}


ofn <- function(eta, ...){
	f1 <- function(t) t*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	f2 <- function(t) exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t))
	v1 <- integrate(f1, lower=-Inf, upper=Inf)$value
	v2 <- integrate(f2, lower=-Inf, upper=Inf)$value
	val <- suppressWarnings(log(v1)) - suppressWarnings(log(v2))
	return(val)
}


v <- sfn(eta=c(0.5,3,5))
print(v)

o <- ofn(eta=c(0.5, 10,2))
print(o)

rv <- optim(par=c(1,0,1.5), fn=ofn, gr=sfn, method="L-BFGS-B", lower=c(0.5,-3,0), upper=c(2,3,5)) 
print(rv)
}