inst/tests/sim-02-1.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

## sim-02-1.R
## 
## Change Delta_n in terms of the sample size 
## in order to show a shrinkage

rm(list=ls())
library(ipeglim)
options(warn=0) # all warning message as errors 2
## options(warn=-1) >> ignore all warnings
## options(warn=0) >> default
## options(warn=1) >> print warnings immediately (make stack empty) 
## options(warn=2) >> make warning behave as error
## 
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/ztpr.R")
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/summary.ztpr.R")
# source("~/Documents/R-ko/countreg/pkg/R/zerotrunc.R")

mu0 <- c(0.288, 0.694, 1.386) 	# performance 
N <- c(5e1,1e2,3e2,1e3)  				# pop. size
alpha <- c(2, 5, 10, 100) 			# shape (gamma) for over dispsersion
# dispersion = c(0.50, 0.20, 0.10, 0.01)
# nb -> pois
lc.lgamma <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), 
	rhs=c(1e-3, -10, 1e-3, -10)) 	# constrained M0
niters <- 5e2
id <- 0
op <- list()

# data and its summary stats
tmp <- simulateYX(N=N[1], Xreg=FALSE, ztrunc=TRUE, param=mu0[2])
y <- tmp$y 	#$

# classical ZTP model
## the use of classical ZTP is important to identify centeral point of M0
fit.ztp <- try(summary(ztpr(formula=y~1, ztrunc=TRUE, dist="poisson"), 
  LM.test=TRUE, HT.est=TRUE))
if(class(fit.ztp) == "try-error") fit.ztp <- FALSE

# imprecise ZTP model
fit <- model(formula=y~0, ztrunc=TRUE, dist="poisson")
fit <- impose(obj=fit, eqns=lc.lgamma)
fit <- update(obj=fit, method="LA", priorType="lgamma")
fit <- summary(fit, HT.est=TRUE)

## senario 1. 
## The true value is right to M0
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), 
	rhs=c(4, -6, 4, -6))
fit <- model(formula=y~0, ztrunc=TRUE, dist="poisson")
fit <- impose(obj=fit, eqns=lc.lgamma)
fit <- update(obj=fit, method="LA", priorType="lgamma")
fit <- summary(fit, HT.est=TRUE)

## What are the mean values when Gamma(alpha, beta) 
## (1,1): extremes
## (9,9): extremes
## (5,5): centeral point (TRUE)
## (3,5): left on x,  (7,5): right on x
## (5,3): below on y, (5,7): top on y
##
## Y~Gamma(a,b) -> E(Y)=a/b, V(Y)=a/b^2
##
## 1-e^{-m} = ppois(0,m, lower.tail=FALSE)
## e^{-m} = dpois(0,m)
##
## Y~ZTP(m)
## E(Y)=m/ppois(0,m,FALSE)
## V(Y)=m*(1-e^{-m}-m*e^{-m})/(1-e^{-m})^2
##     =m*(ppois(0,m,FALSE)-m*dpois(0,m))/ppois(0,m,FALSE)^2
##