|
# NEGATIVE BINOMIAL/GAMMA MODEL |
|
# COVERAGE COMPUTED EXAMPLE |
|
# GAMMA = 1 AND W = 1 |
|
rm(list = ls()) |
|
library(xtable) |
|
library(LearnBayes) |
|
source("https://gist.githubusercontent.com/eliardocosta/85781d93cbaae35411eed74372c7c2d4/raw/9e0703312f1d83f69f91d7685829533e3ffc90a4/logp.nbgam.R") |
|
source("https://gist.githubusercontent.com/eliardocosta/85781d93cbaae35411eed74372c7c2d4/raw/9e0703312f1d83f69f91d7685829533e3ffc90a4/rlam.nbgam.R") |
|
gam <- 1 |
|
w <- 1 |
|
phi <- 8 |
|
lam0 <- 10 |
|
eps <- 2 # prior var = eps^2 = 4 |
|
theta0 <- (lam0/eps)^2; theta0 |
|
#---------------------------------------------------------- |
|
#-- LAMBDA = 7 |
|
#--------------- |
|
# c = 0.01 |
|
# SIMULATING THE COUNTS |
|
seed <- 123456 |
|
n <- 51 # Table 2 |
|
lamR <- 7 # real concentration |
|
set.seed(seed) |
|
counts <- rnbinom(n, mu = w*lamR, size = phi); counts |
|
sum(counts) # 319 |
|
hist(counts, prob = TRUE) |
|
print(xtable(matrix(c(counts, NA), 4, 13, byrow = TRUE), |
|
digits = 0), include.rownames = FALSE, |
|
include.colnames = FALSE) |
|
|
|
x <- counts |
|
set.seed(seed) |
|
lampos1 <- rlam.nbgam(N = 1E4, x = x, phi = phi, |
|
lam0 = lam0, theta0 = theta0, |
|
w = w, burnin = 1E3, thin = 10) |
|
round(lampos1$accept, 2) # acceptance rate ~ 0.48 |
|
par(mfrow = c(2, 1)) |
|
plot.ts(lampos1$lam, xlab = "iteration", ylab = "") # traceplot |
|
acf(lampos1$lam, main = "") # ACF |
|
|
|
medpos <- mean(lampos1$lam) |
|
varpos <- var(lampos1$lam) |
|
a <- medpos - sqrt(varpos/gam); round(a, 2) # 6.10 |
|
b <- medpos + sqrt(varpos/gam); round(b, 2) # 7.05 |
|
|
|
# coverage probability ~ 0.68 |
|
count <- 0 |
|
for (k in 1:length(lampos1$lam)) { |
|
if (lampos1$lam[k] >= a && lampos1$lam[k] <= b) count <- count + 1 |
|
} |
|
count/length(lampos1$lam) # cov. probability |
|
|
|
#---------------------------------------------------------- |
|
#-- LAMBDA = 10 |
|
#--------------- |
|
# c = 0.01 |
|
# SIMULATING THE COUNTS |
|
seed <- 123456 |
|
n <- 51 # Table 2 |
|
lamR <- 10 # real concentration |
|
set.seed(seed) |
|
counts <- rnbinom(n, mu = w*lamR, size = phi); counts |
|
sum(counts) # 523 |
|
par(mfrow = c(1, 1)) |
|
hist(counts, prob = TRUE) |
|
print(xtable(matrix(c(counts, NA), 4, 13, byrow = TRUE), |
|
digits = 0), include.rownames = FALSE, |
|
include.colnames = FALSE) |
|
x <- counts |
|
set.seed(seed) |
|
lampos2 <- rlam.nbgam(N = 1E4, x = x, phi = phi, |
|
lam0 = lam0, theta0 = theta0, |
|
w = w, burnin = 1E3, thin = 10) |
|
round(lampos2$accept, 2) # acceptance rate ~ 0.58 |
|
par(mfrow = c(2, 1)) |
|
plot.ts(lampos2$lam, xlab = "iteration", ylab = "") # traceplot |
|
acf(lampos2$lam, main = "") # ACF |
|
|
|
medpos <- mean(lampos2$lam) |
|
varpos <- var(lampos2$lam) |
|
a <- medpos - sqrt(varpos/gam); round(a, 2) # 9.61 |
|
b <- medpos + sqrt(varpos/gam); round(b, 2) # 10.89 |
|
|
|
# coverage probability ~ 0.68 |
|
count <- 0 |
|
for (k in 1:length(lampos2$lam)) { |
|
if (lampos2$lam[k] >= a && lampos2$lam[k] <= b) count <- count + 1 |
|
} |
|
count/length(lampos2$lam) # cov. probability |
|
|
|
#---------------------------------------------------------- |
|
#-- LAMBDA = 13 |
|
#--------------- |
|
# c = 0.01 |
|
# SIMULATING THE COUNTS |
|
seed <- 123456 |
|
n <- 51 # Table 2 |
|
lamR <- 13 # real concentration |
|
set.seed(seed) |
|
counts <- rnbinom(n, mu = w*lamR, size = phi); counts |
|
sum(counts) # 641 |
|
par(mfrow = c(1, 1)) |
|
hist(counts, prob = TRUE) |
|
print(xtable(matrix(c(counts, NA), 4, 13, byrow = TRUE), |
|
digits = 0), include.rownames = FALSE, |
|
include.colnames = FALSE) |
|
x <- counts |
|
set.seed(seed) |
|
lampos3 <- rlam.nbgam(N = 1E4, x = x, phi = phi, |
|
lam0 = lam0, theta0 = theta0, |
|
w = w, burnin = 1E3, thin = 10) |
|
round(lampos3$accept, 2) # acceptance rate ~ 0.62 |
|
par(mfrow = c(2, 1)) |
|
plot.ts(lampos3$lam, xlab = "iteration", ylab = "") # traceplot |
|
acf(lampos3$lam, main = "") # ACF |
|
|
|
medpos <- mean(lampos3$lam) |
|
varpos <- var(lampos3$lam) |
|
a <- medpos - sqrt(varpos/gam); round(a, 2) # 11.58 |
|
b <- medpos + sqrt(varpos/gam); round(b, 2) # 13.04 |
|
|
|
# coverage probability ~ 0.68 |
|
count <- 0 |
|
for (k in 1:length(lampos3$lam)) { |
|
if (lampos3$lam[k] >= a && lampos3$lam[k] <= b) count <- count + 1 |
|
} |
|
count/length(lampos3$lam) # cov. probability |