Nothing
## imprecise/demo/lgamma-ns.R
##
## show the numerical stabiability on estimating the canonical parameter of
## Poisson model using numerical algorithms of 'LA', 'MH', 'IS', 'AQ',
## as well as 'NA'.
##
## 2013.10.21.
library(ipeglim)
rm(list=ls())
# random variates generation
### try to change 'ztrunc' from 'TRUE' to 'FALSE'
y <- rvg4yx(N=5e1, pars=1, ztrunc=TRUE, xreg=FALSE, kind="pois")$y
table(y)
# model to be fitted
### the value of 'ztrunc' in the function 'model()' should be consistent with
### the value of 'ztrunc' in the function 'rvg4yx()'.
m2fit <- model(formula=y~0, ztrunc=TRUE, dist="poisson")
# box-constrained normal imprecise prior
lc0 <- list(lhs=rbind(c(1,0), c(-1,0), c(0,1), c(0,-1)), rhs=c(-6, -6, 0.5, -2.5))
cmfit <- impose(obj=m2fit, eqns=lc0)
# apply numerical methods
t0 <- proc.time()
ola <- summary(update(obj=cmfit, method="LA", apriori="normal"))
t1 <- proc.time()
omh <- summary(update(obj=cmfit, method="MH", apriori="normal"))
t2 <- proc.time()
ois <- summary(update(obj=cmfit, method="IS", apriori="normal"))
t3 <- proc.time()
oaq <- summary(update(obj=cmfit, method="AQ", apriori="normal"))
t4 <- proc.time()
# comparison
tt <- list(t0, t1, t2, t3, t4)
elapsed <- diff(do.call(c, lapply(tt, "[[", 3)))
tb <- cbind(ola$est, omh$est, ois$est, oaq$est)
colnames(tb) <- c("LA", "MH", "IS", "AQ")
printCoefmat(tb, tst.ind=1:4)
# quality check
ois$ess # 'IS'
omh$arate # 'MH'
ola$sratio # 'LA'
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.