Nothing
## 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.