Nothing
lvl <- seq(from=0.25, to=0.75, by=0.25)
mu <- -log(-lvl+1)
theta <- log(mu)
## R> mu
## [1] 0.2876821 0.6931472 1.3862944
##
## Testing
if(FALSE){
ppois(q=0, lambda=mu, lower.tail=FALSE)
## Interpretation of 0.25
## 75% of cases are missed (i.e., this system captures 25% of cases)
}
## Use this vector of 'mu' as the null values that gives 0.25%, 0.50%, 0.75%
## of capturing performance
mu0 <- c(0.288, 0.694, 1.386)
##
## Models without covariats
## imprecise gamma prior
## (a,b) = { 0.01 <= a <= 10, 0.01 <= b <= 10}
##
## 'alpha' (shape parameter of gamma-poisson parametrization)
## 'alpha' goes to 'infty' (i.e., 'dispersion' goes to '0');
## thus, 'f(y)' goes to Poisson from NegBinom.
##
## meta-parameters
##
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()
for(idx.s in alpha){
for(idx.m in mu0){
for(idx.n in N){
for(i in seq_len(niters)){
id <- id + 1
# data and its summary stats
tmp <- simulateYX(N=idx.n, Xreg=FALSE, ztrunc=TRUE,
param=idx.m, shape=idx.s)
y <- tmp$y #$
# classical ZTP model
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
fit.ztnb <- try(summary(ztpr(formula=y~1, ztrunc=TRUE, dist="nbinom"),
HT.est=TRUE))
if(class(fit.ztnb) == "try-error") fit.ztnb <- 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)
rvals <- list(d0=tmp, ztp=fit.ztp, ztnb=fit.ztnb, fit=fit)
rm(tmp, y, fit.ztp, fit.ztnb, fit)
op[[id]] <- rvals
print(paste("id=", id, "; alpha=", idx.s,": mu0=", idx.m, ": N=", idx.n, ": niter=",
i, sep=""))
}
}
}
}
warnings()
## Potential errors catched by 'try()'
##
## Error in optim(par = init, fn = fnllik, gr = grllik, method = method, :
## L-BFGS-B needs finite values of 'fn'
## Error in solve.default(as.matrix(fit$hessian)) :
## Lapack routine dgesv: system is exactly singular: U[1,1] = 0
save(op,
file="~/Documents/PhD/ipeglim/thesis/chapters/ch04/sim-01-lgamma.RData")
op1 <- do.call(rbind, op)
# lapply(op, with, ztp$N)
if(FALSE){
# z <- range(sfit$est)
# zlim <- c(z[1]-0.05, z[2]+0.05)
# plot(sfit, zlim=zlim)
# cmf1 <- setGrid(cmf, len=10)
# fit1 <- update(obj=cmf1, method="LA", priorType="lgamma")
# sfit1 <- summary(fit1, HT.est=TRUE)
# z1 <- range(sfit1$est)
# zlim1 <- c(z1[1]-0.05, z1[2]+0.05)
# plot(sfit1, zlim=zlim1, xlim=c(-1,11), ylim=c(-1,11), clevels=10)
################################################################################
## Conventional model
################################################################################
## Final simulation test
# log(mu) = 1 + x1 + x2
#
# What if the skewness about X? (in multivariate)
# the multivariate skew-normal (MSN)
rm(list=ls())
library(ipeglim)
# source("~/Documents/PhD/ipeglim/pkg/ipeglim/R/simulateYX.R")
# btb <- expand.grid(N=c(100), b0=seq(-1.0, 1.0, len=3), b1=1, b2=-1, x1sd=seq(0.5, 1.5, len=3), x2sd=seq(0.5, 1.5, len=3), r=seq(0.0, 0.8, len=3))
# show a robustnesss over correlations
btb <- expand.grid(N=c(100), b0=seq(-1.0, 1.0, len=3), b1=1, b2=-1, x1sd=c(1), x2sd=c(1), r=seq(0.0, 0.8, len=3))
op <- list()
niters <- 1
i <- 1
for(j in 1:nrow(btb)){
while(i <= niters){
Xmean <- c(1, 0, 0)
Xsd <- c(0, btb$x1sd[j], btb$x2sd[j])
XR <- matrix(c(1,0,0, 0,1,btb$r[j], 0,btb$r[j],1), ncol=length(Xmean))
N0 <- btb$N[j]
bcoef <- as.vector(c(btb$b0[j], btb$b1[j], btb$b2[j]))
Xvcv <- diag(Xsd) %*% XR %*% t(diag(Xsd))
mu <- exp((Xvcv+diag(Xmean))%*%bcoef)[1]
miss <- dpois(0, lambda=mu)
tmp <- simulateYX(N=N0, Xreg=TRUE, param=c(btb$b0[j], btb$b1[j], btb$b2[j]), Xstr=list(type="mvnorm", mean=Xmean, sd=Xsd, R=XR), ztrunc=TRUE)
dat <- tmp$yX
yt <- c(tmp$y)
fit <- ztpr(formula=y~x1+x2, data=dat, dist="poisson", ztrunc=TRUE, method="BFGS")
fit <- summary(fit)
sds <- sqrt(diag(fit$vcv))
cover <- all(c(fit$cil < N0, N0 < fit$ciu))
op[[niters*(j-1)+i]] <- c(N0=N0, b0=btb$b0[j], b1=btb$b1[j], b2=btb$b2[j], x1sd=btb$x1sd[j], x2sd=btb$x2sd[j], r=btb$r[j], n=nrow(dat), nr=nrow(dat)/N0, scm1=mean(yt), scm2=var(yt), sg1=skewness(yt), scv=cv(yt), N=fit$N, Nr=(fit$N)/N0, cil=fit$cil, cilr=fit$cil/N0, ciu=fit$ciu, ciur=fit$ciu/N0, cover=cover, cf0=fit$cfs[1], cf1=fit$cfs[2], cf2=fit$cfs[3], cf0se=sds[1], cf1se=sds[2], cf2se=sds[3])
i <- i + 1
}
i <- 1
print(j)
}
out <- as.data.frame(do.call(rbind, op))
out
ss <- subset(out, subset=(b0==0 & x1sd==1 & x2sd==1 & r==0))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.