inst/thesis/fig/o3ef.R

## ipeglim/inst/thesis/fig/o3ef.R
## 
## o3ef = order 3 exponential family of distribution (i.e., normal prior)

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

lmu <- c(-3, 2, 1, 0, -1, -2, 3)
e0 <- seq(0,2,0.25)
e2 <- seq(0,5,0.3)
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
z <- range(e1, na.rm=TRUE)

tobj <- wireframe(e1~e0*e2, data=hcube, 
	group=lmu, 
  screen=list(y=60, x=15, y=20), 
#  zoom=0.85, 
  pretty=TRUE, 
  xlim=c(0,2),
  ylim=c(0,5),
  zlim=c(z[1],z[2]),
  main=expression("Level Sets of "~E~"("~theta~"|"~eta[0]~","~eta[1]~","~eta[2]~")"),
  sub=expression("log"~"("~mu~")"~"="~-3 ~"to"~3~"by"~1~"from left to right"),
  xlab=expression(eta[0]), 
  ylab=expression(eta[2]), 
  zlab=expression(eta[1]), 
  par.settings = list(axis.line = list(col = "transparent")),
  par.box=list(col='black', lwd=0.25), 
	scales=list(arrows=FALSE, col="black"),  # show axese with ticks (z.ticks=2), 
	type="l"
	)
print(tobj)


if(FALSE){
# xyplot(e1~e2|e0, data=hcube[hcube$e0 %in% c(0,0.25),])
}

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.