inst/thesis/fig/o3ef_snapshot.R

## ipeglim/inst/thesis/fig/o3ef_snapshot.R
## 
## Three snapshots from the 3D levelset
## 
# hparam <- c(nu, tau2)
#           eta2         eta1          eta0
# eta <- c(0.5/tau2, sum(y)+nu/tau2, length(y))
#          theta^2       theta       exp(theta)
# -6 <= nu <= 6
# 0.5 <= tau2 <= 2.5

rm(list=ls())
library(ipeglim)
set.seed(16979238)

# 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")

# par(mfrow=c(2,2))

# if(FALSE){
# lmu <- c(-3, 2, 1, 0, -1, -2, 3)
lmu <- seq(-1, 1, 0.2)
# for actual compilation should use this 0.05
# e0 <- seq(0, 3, 0.05)
# e2 <- seq(1e-3, 1, 0.05)
e0 <- seq(0, 3, 0.25)
e2 <- seq(1e-3, 3, 0.25)
hcube <- expand.grid(e0=e0, e2=e2, lmu=lmu)

e1 <- apply(hcube, 1, function(x){
  x <- as.vector(x)
  rt.e1 <- tryCatch(uniroot(eta1fn, lower=-20, upper=30, tol=1e-4, eta0=x[1], eta2=x[2], logm=x[3])$root, error=function(e) return(NaN))
  return(rt.e1)
  })
hcube$e1 <- e1

save(list=ls(), file="~/Documents/PhD/ipeglim/thesis/fig/o3ef_snapshot.RData")

Try the ipeglim package in your browser

Any scripts or data that you put into this service are public.

ipeglim documentation built on May 2, 2019, 4:31 p.m.