famtnb.R
library(aster)
alpha <- 2.222
k <- 2
ifam <- fam.truncated.negative.binomial(alpha, k)
aster:::setfam(list(ifam))
aster:::getfam()
aster:::clearfam()
##### usual size theta #####
p <- seq(0.9, 0.1, -0.1)
theta <- log(1 - p)
qq <- exp(theta)
pp <- (- expm1(theta))
all.equal(p, pp)
all.equal(pp, 1 - qq)
beeta <- pnbinom(k + 1, size = alpha, prob = pp, lower.tail = FALSE) /
dnbinom(k + 1, size = alpha, prob = pp)
zeroth <- double(length(theta))
first <- double(length(theta))
second <- double(length(theta))
for (i in seq(along = theta)) {
zeroth[i] <- famfun(ifam, 0, theta[i])
first[i] <- famfun(ifam, 1, theta[i])
second[i] <- famfun(ifam, 2, theta[i])
}
all.equal(zeroth, alpha * (- log(pp)) +
pnbinom(k, size = alpha, prob = pp, lower.tail = FALSE, log.p = TRUE))
all.equal(first, alpha * qq / pp + (k + 1) / (1 + beeta) / pp)
all.equal(second, alpha * qq / pp^2 - (k + 1) / (1 + beeta) / pp^2 *
(- qq + (k + 1 + alpha) * qq / (1 + beeta) +
(alpha - pp * (k + 1 + alpha)) * beeta / (1 + beeta)))
##### usual size theta (continued) #####
##### check by numerical derivative #####
epsilon <- 1e-8
zeroth.minus <- zeroth
zeroth.plus <- zeroth
for (i in seq(along = theta)) {
zeroth.minus[i] <- famfun(ifam, 0, theta[i] - epsilon)
zeroth.plus[i] <- famfun(ifam, 0, theta[i] + epsilon)
}
all.equal(first, (zeroth.plus - zeroth.minus) / (2 * epsilon),
tolerance = sqrt(epsilon))
first.minus <- first
first.plus <- first
for (i in seq(along = theta)) {
first.minus[i] <- famfun(ifam, 1, theta[i] - epsilon)
first.plus[i] <- famfun(ifam, 1, theta[i] + epsilon)
}
all.equal(second, (first.plus - first.minus) / (2 * epsilon),
tolerance = sqrt(epsilon))
##### very large negative theta #####
rm(p)
theta <- seq(-100, -10, 10)
qq <- exp(theta)
pp <- (- expm1(theta))
beeta.up <- pnbinom(k + 1, size = alpha, prob = pp, lower.tail = FALSE)
beeta.dn <- dnbinom(k + 1, size = alpha, prob = pp)
beeta <- beeta.up / beeta.dn
beeta[beeta.up == 0] <- 0
zeroth <- double(length(theta))
first <- double(length(theta))
second <- double(length(theta))
for (i in seq(along = theta)) {
zeroth[i] <- famfun(ifam, 0, theta[i])
first[i] <- famfun(ifam, 1, theta[i])
second[i] <- famfun(ifam, 2, theta[i])
}
all.equal(zeroth, alpha * (- log(pp)) +
pnbinom(k, size = alpha, prob = pp, lower.tail = FALSE, log.p = TRUE))
all.equal(first, alpha * qq / pp + (k + 1) / (1 + beeta) / pp)
all.equal(second, alpha * qq / pp^2 - (k + 1) / (1 + beeta) / pp^2 *
(- qq + (k + 1 + alpha) * qq / (1 + beeta) +
(alpha - pp * (k + 1 + alpha)) * beeta / (1 + beeta)))
##### very small negative theta #####
theta <- (- 10^(- c(1:9, seq(10, 100, 10))))
qq <- exp(theta)
pp <- (- expm1(theta))
beeta <- pnbinom(k + 1, size = alpha, prob = pp, lower.tail = FALSE) /
dnbinom(k + 1, size = alpha, prob = pp)
zeroth <- double(length(theta))
first <- double(length(theta))
second <- double(length(theta))
for (i in seq(along = theta)) {
zeroth[i] <- famfun(ifam, 0, theta[i])
first[i] <- famfun(ifam, 1, theta[i])
second[i] <- famfun(ifam, 2, theta[i])
}
all.equal(zeroth, alpha * (- log(pp)) +
pnbinom(k, size = alpha, prob = pp, lower.tail = FALSE, log.p = TRUE))
all.equal(first, alpha * qq / pp + (k + 1) / (1 + beeta) / pp)
all.equal(second, alpha * qq / pp^2 - (k + 1) / (1 + beeta) / pp^2 *
(- qq + (k + 1 + alpha) * qq / (1 + beeta) +
(alpha - pp * (k + 1 + alpha)) * beeta / (1 + beeta)))
##### random #####
nind <- 50
theta <- rep(- 1.75, nind)
pred <- 0
fam <- 1
root <- rep(1:3, length.out = nind)
theta <- cbind(theta)
root <- cbind(root)
qq <- exp(theta)
pp <- (- expm1(theta))
mu <- alpha * qq / pp
set.seed(42)
rout <- raster(theta, pred, fam, root, famlist = list(ifam))
set.seed(42)
rout.too <- rktnb(nind, alpha, k, mu, root)
all.equal(as.numeric(rout), rout.too)