inst/thesis/fig/normal-test.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## with normal prior
## optimization is test

rm(list=ls())

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

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

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

	gd0 <- function(t) -t*exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t)+t)
	vgd0 <- integrate(gd0, lower=-Inf, upper=Inf)$value
	hd0 <- function(t) -exp(-eta[3]*t^2+eta[2]*t-eta[1]*exp(t)+t)
	vhd0 <- integrate(hd0, lower=-Inf, upper=Inf)$value
	ede0 <-(vgd0*vh-vg*vhd0)/vh^2
	
	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))
	val <- suppressWarnings(v1/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,3), fn=ofn, gr=sfn, method="L-BFGS-B", lower=c(0.5,-3,2), upper=c(1,3,5)) 
rv <- optim(par=c(1,0,3), fn=ofn, gr=sfn, method="L-BFGS-B", lower=c(0.5,-3,2), upper=c(1,3,5), control=list(fnscale=-1)) 
print(rv)
## what if 'constrOptim()' is used?
# }


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)
}