inst/thesis/fig/lg_optimizers.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## 
# set.seed(18372342)
rm(list=ls())
library(ipeglim)
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/model.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/print_imprecise.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/plot_imprecise.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/summary_imprecise.R")

y <- rvg4yx(N=1e0, pars=2, ztrunc=FALSE, xreg=FALSE, kind="poisson")$y
# table(y)

fn <- function(x, ...){
  a <- x[1]
  b <- x[2]
  n <- length(y)
  ybar <- mean(y)
  val <- digamma(a+n*ybar) - log(b+n)
  return(val)
}

gfn <- function(x, ...){
  a <- x[1]
  b <- x[2]
  n <- length(y)
  ybar <- mean(y)
  val <- c( trigamma(a+n*ybar), -1/(b+n))
  return(val)
}

m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)

## my test set 
lc0 <- list(lhs=rbind(c(1,1), c(-1,-1), c(-1,1), c(1,-1), c(0,-1), c(0,1)), rhs=c(3,-7,-3,-1,-3,1)*5)
xtm0 <- iprior(eqns=lc0)
cc <- colMeans(do.call(rbind, xtm0$xtms))

## scenario 1.
# lc1 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2))), rhs=c(0,0,-10,-10))
cmfit <- iprior(obj=m2fit, eqns=lc0, verbose=FALSE)
plot(cmfit)
sop0 <- summary(update(obj=cmfit, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

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

idx <- c(sop0$inf.idx, sop0$sup.idx)
tb0 <- cbind(do.call(rbind, lapply(sop0$m1[idx], with, hparam)), c(sop0$inf, sop0$sup), idx)
tb1 <- cbind(rbind(mini$par, maxi$par), c(mini$val, maxi$val), idx)
tbs0 <- rbind(tb0, tb1)
colnames(tbs0) <- c("$\\alpha$", "$\\beta$", "$E(\\theta)$", "XID")
rownames(tbs0) <- c("$\\Lprevision{\\theta}{}$", "$\\Uprevision{\\theta}{}$", "Min.Optim", "Max.Optim")
print(xtable(tbs0, digits=3, caption="Optimizers at extreme points"), sanitize.text.function=function(x){x})