inst/thesis/fig/lg_linear_2.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## 
## 

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

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

set.seed(18372342)
y <- rvg4yx(N=1, pars=2, ztrunc=FALSE, xreg=FALSE, kind="poisson")$y
table(y)
m2fit <- model(formula=y~0, ztrunc=FALSE, dist="poisson", verbose=FALSE)

# scenarios
lc1 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2))), rhs=c(0,0,-10,-10))
lc2 <- list(lhs=rbind(diag(rep(1,2)), diag(rep(-1,2)), c(-1,1)), rhs=c(1,1,-8,-8,-4))
lc3 <- 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)*2)
lc4 <- list(x=5,y=5,r=5,len=25)

## R0 - 1st item
cmfit1 <- iprior(obj=m2fit, eqns=lc1, verbose=FALSE)
# cmfit1 <- set.grid(obj=cmfit1, len=10) # xtms2eqns is not ok
cmfit2 <- iprior(obj=m2fit, eqns=lc2, verbose=FALSE)
cmfit3 <- iprior(obj=m2fit, eqns=lc3, verbose=FALSE)
cmfit4 <- iprior(obj=m2fit, circle=lc4, verboser=FALSE)

sop1 <- summary(update(obj=cmfit1, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)
min1 <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc1$lhs, ci=lc1$rhs)
max1 <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc1$lhs, ci=lc1$rhs, control=list(fnscale=-1))

tb01min <- c(sop1$m1[[sop1$inf.idx]]$hparam, sop1$m1[[sop1$inf.idx]]$theta)
tb11min <- c(min1$par, min1$val)
tb01max <- c(sop1$m1[[sop1$sup.idx]]$hparam, sop1$m1[[sop1$sup.idx]]$theta)
tb11max <- c(max1$par, max1$val)
tab1 <- cbind(tb01min, tb11min, tb01max, tb11max)
round(tab1,3)


## R0 - 2st item
sop2 <- summary(update(obj=cmfit2, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)
min2 <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc2$lhs, ci=lc2$rhs)
max2 <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc2$lhs, ci=lc2$rhs, control=list(fnscale=-1))

tb02min <- c(sop2$m1[[sop2$inf.idx]]$hparam, sop2$m1[[sop2$inf.idx]]$theta)
tb12min <- c(min2$par, min2$val)
tb02max <- c(sop2$m1[[sop2$sup.idx]]$hparam, sop2$m1[[sop2$sup.idx]]$theta)
tb12max <- c(max2$par, max2$val)
tab2 <- cbind(tb02min, tb12min, tb02max, tb12max)
round(tab2,3)


## R0 - 3rd item
sop3 <- summary(update(obj=cmfit3, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)
min3 <- constrOptim(theta=c(5,5), f=fn, grad=gfn, ui=lc3$lhs, ci=lc3$rhs)
max3 <- constrOptim(theta=c(5,5), f=fn, grad=gfn, ui=lc3$lhs, ci=lc3$rhs, control=list(fnscale=-1))

tb03min <- c(sop3$m1[[sop3$inf.idx]]$hparam, sop3$m1[[sop3$inf.idx]]$theta)
tb13min <- c(min3$par, min3$val)
tb03max <- c(sop3$m1[[sop3$sup.idx]]$hparam, sop3$m1[[sop3$sup.idx]]$theta)
tb13max <- c(max3$par, max3$val)
tab3 <- cbind(tb03min, tb13min, tb03max, tb13max)
round(tab3,3)

## R0 - 4th item
source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/xtms2eqns.R")
# par(mfrow=c(2,2))
x2e <- xtms2eqns(obj=cmfit4, show=TRUE, type="circle")
# plot(iprior(eqns=list(lhs=x2e$ui, rhs=x2e$ci)))
sop4 <- summary(update(obj=cmfit4, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

min4 <- constrOptim(theta=x2e$cc, f=fn, grad=gfn, ui=x2e$ui, ci=x2e$ci)
max4 <- constrOptim(theta=x2e$cc, f=fn, grad=gfn, ui=x2e$ui, ci=x2e$ci, control=list(fnscale=-1))

tb04min <- c(sop4$m1[[sop4$inf.idx]]$hparam, sop4$m1[[sop4$inf.idx]]$theta)
tb14min <- c(min4$par, min4$val)
tb04max <- c(sop4$m1[[sop4$sup.idx]]$hparam, sop4$m1[[sop4$sup.idx]]$theta)
tb14max <- c(max4$par, max4$val)
tab4 <- cbind(tb04min, tb14min, tb04max, tb14max)
round(tab4,3)