inst/thesis/fig/lg_optimizers_exam.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## 
## 
set.seed(18372342)
rm(list=ls())
library(ipeglim)

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

if(FALSE){
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))
xtm0 <- iprior(eqns=lc0)

## 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)
sop0 <- summary(update(obj=cmfit, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

mini <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc0$lhs, ci=lc0$rhs)
maxi <- constrOptim(theta=c(3,2), 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], "[[", "hparam")), c(sop0$inf, sop0$sup))
tb1 <- cbind(rbind(mini$par, maxi$par), c(mini$val, maxi$val))
tbs0 <- rbind(tb0, tb1)
tbs0
# 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})
}


#########################
## These constraints were used in the file "lg_Dn.R"
## 2013-NOV-09
y <- rvg4yx(N=1, pars=2, ztrunc=FALSE, xreg=FALSE, kind="poisson")$y

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(diag(rep(1,2)), diag(rep(-1,2)), c(1,-1)), rhs=c(1,1,-8,-8,-3))
lc4 <- list(x=5,y=5,r=5,len=20)

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

if(FALSE){
# scenario 1
cmfit1 <- iprior(obj=m2fit, eqns=lc1, verbose=FALSE)
sop1 <- summary(update(obj=cmfit1, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

mini <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc1$lhs, ci=lc1$rhs)
maxi <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc1$lhs, ci=lc1$rhs, control=list(fnscale=-1))

idx <- c(sop1$inf.idx, sop1$sup.idx)
tb0 <- cbind(do.call(rbind, lapply(sop1$m1[idx], "[[", "hparam")), c(sop1$inf, sop1$sup))
tb1 <- cbind(rbind(mini$par, maxi$par), c(mini$val, maxi$val))
tbs0 <- rbind(tb0, tb1)
round(tbs0,3)

# scenario 2
cmfit2 <- iprior(obj=m2fit, eqns=lc2, verbose=FALSE)
sop2 <- summary(update(obj=cmfit2, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

mini <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc2$lhs, ci=lc2$rhs)
maxi <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc2$lhs, ci=lc2$rhs, control=list(fnscale=-1))

idx <- c(sop2$inf.idx, sop2$sup.idx)
tb0 <- cbind(do.call(rbind, lapply(sop2$m1[idx], "[[", "hparam")), c(sop2$inf, sop2$sup))
tb1 <- cbind(rbind(mini$par, maxi$par), c(mini$val, maxi$val))
tbs0 <- rbind(tb0, tb1)
round(tbs0,3)



# scenario 3
cmfit3 <- iprior(obj=m2fit, eqns=lc3, verbose=FALSE)
sop3 <- summary(update(obj=cmfit3, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)

mini <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc3$lhs, ci=lc3$rhs)
maxi <- constrOptim(theta=c(3,2), f=fn, grad=gfn, ui=lc3$lhs, ci=lc3$rhs, control=list(fnscale=-1))

idx <- c(sop3$inf.idx, sop3$sup.idx)
tb0 <- cbind(do.call(rbind, lapply(sop3$m1[idx], "[[", "hparam")), c(sop3$inf, sop3$sup))
tb1 <- cbind(rbind(mini$par, maxi$par), c(mini$val, maxi$val))
tbs0 <- rbind(tb0, tb1)
round(tbs0,3)
}

# scenario 4
lc4 <- list(x=5,y=5,r=5,len=101) # odd number of lengths is fine
cmfit4 <- iprior(obj=m2fit, circle=lc4, verbose=FALSE)
sop4 <- summary(update(obj=cmfit4, method="AS", apriori="lgamma", silent=TRUE), silent=TRUE)


xtms2eqns <- function(obj=cmfit, show=FALSE){
	
	xy1 <- xy <- cmfit4$xtms
	xy0 <- do.call(rbind, xy)
	xy1[[length(xy)+1]] <- xy[[1]]
	xy1 <- do.call(rbind, xy1)

	cc <- colMeans(xy0)

	ui <- ci <- list()
	if(show) plot(0, xlim=c(0,10), ylim=c(0,10))

	for(i in 1:length(xy)){
		x1 <- xy1[i,1]
		y1 <- xy1[i,2]
		x2 <- xy1[i+1,1]
		y2 <- xy1[i+1,2]
		b <- (y2-y1)/(x2-x1)
		a <- -b*x1 + y1
		
		if ( (a+b*cc[1] <= cc[2]) ) {
			ui[[i]] <- c(1,-b)
			ci[[i]] <- a
			if (show) abline(a=a,b=b, col="red")
		} else {
			ui[[i]] <- -c(1,-b)
			ci[[i]] <- -a
			if (show) abline(a=a,b=b, col="blue")
		}
		if (show) polygon(xy0[chull(xy0),], col=alpha("azure2", 0.2), border=alpha("darkblue",0.5))
		if (show) points(x=cc[1], y=cc[2])
		if (show) points(xy0)
	}
	rvals <- list(cc=cc, ui=do.call(rbind, ui), ci=do.call(c, ci))
}

ilc <- xtms2eqns(obj=cmfit4, show=TRUE)
mini <- constrOptim(theta=ilc$cc, f=fn, grad=gfn, ui=ilc$ui, ci=ilc$ci)
maxi <- constrOptim(theta=ilc$cc, f=fn, grad=gfn, ui=ilc$ui, ci=ilc$ci, control=list(fnscale=-1))

idx <- c(sop4$inf.idx, sop4$sup.idx)
tb0 <- cbind(do.call(rbind, lapply(sop4$m1[idx], "[[", "hparam")), c(sop4$inf, sop4$sup))
tb1 <- cbind(rbind(mini$par, maxi$par), c(mini$val, maxi$val))
tbs0 <- rbind(tb0, tb1)
print(round(tbs0,3))