The function dnbmix() constructs a (vectorized) probability mass
function from a matrix of (simulated) means and corresponding size
parameters, as a function of the time point (row of means) within
the simulation period. The distribution at each time point is obtained
as a mixture of negative binomial (or Poisson) distributions.
dnbmix(means, size = NULL)Arguments
- means
a
n.aheadxn.simmatrix of means.- size
the dispersion parameter of the
dnbinom()distribution orNULL(Poisson forecasts). Can also be time-varying (of lengthn.ahead).
Value
a function(x, tp = 1, log = FALSE), which takes a vector of
counts x and calculates the (log-)probabilities of observing
each of these numbers at the tp'th time point of the simulation
period (indexing the rows of means).
See also
logs_nbmix() where this function is used.
Examples
## a GLARMA example
library("glarma")
y <- as.vector(CHILI)
## fit a simple NegBin-GLARMA model
X <- t(sapply(2*pi*seq_along(y)/52.1775,
function (x) c(sin = sin(x), cos = cos(x))))
X <- cbind(intercept = 1, X)
fit <- glarma(y = y[1:883], X = X[1:883,], type = "NegBin", phiLags = 1)
## simulate the last four weeks (only 500 runs, for speed)
set.seed(1)
means <- replicate(500, {
forecast(fit, n.ahead = 4, newdata = X[884:887,], newoffset = rep(0,4))$mu
})
## derive the weekly forecast distributions
dfun <- dnbmix(means, coef(fit, type = "NB"))
dfun(4000, tp = 1)
#> [1] 0.0001458309
dfun(4000, tp = 4)
#> [1] 7.550036e-05
curve(dfun(x, tp = 4), 0, 30000, type = "h",
main = "4-weeks-ahead forecast",
xlab = "No. infected", ylab = "Probability")
# \dontshow{
## verify distribution at the first time point (i.e., one-step-ahead NegBin)
stopifnot(all.equal(
dfun(0:100, tp = 1),
dnbinom(0:100,
mu = forecast(fit, n.ahead=1, newdata=X[884,,drop=FALSE])$mu,
size = coef(fit, type = "NB"))
))
## check that we have a probability distribution at the second time point
.xgrid <- seq(0, 200000, by = 500)
stopifnot(abs(1 -
integrate(approxfun(.xgrid, dfun(.xgrid, tp = 2)), 0, 200000)$value
) < 0.01)
# }