I want to run a simulation study and I need to simulate data from the Conway–Maxwell–Poisson distribution.
However, it seems like the probability mass function is not available in closed-form, so, I do not know how to simulate from this distribution.
Is there a method that I can use to simulate from this distribution?
Is there any R package that implements this simulation?
The CMP distribution generalizes the Poisson distribution. It is a discrete distribution on the natural numbers $0,1,2,\ldots$ with probabilities
$$\Pr(x;\lambda,\nu)\ \propto\ \frac{\lambda^x}{\left(x!\right)^\nu}$$
where $\lambda \gt 0$ and $\nu \ge 0$ are its parameters. When $\nu = 1,$ this is the Poisson distribution with parameter $\lambda.$
Here are some examples, plotted as probability mass functions (bar heights equal the probabilities):
The constant of proportionality for the probabilities is not easily computed, especially when $\nu$ is not integral. The challenging cases to simulate are for $\nu\ll 1,$ as I will show with an example. Here is a plot of the counts in a simulation of a million iid values from the CMP distribution with parameters $\lambda = 50, \nu = 1/4:$
The red curve plots the expected counts (it's proportional to the underlying probabilities).
The range of the simulated values is from less than 6.23 million to over 6.27 million, comprising over 30,000 distinct values. Within this range the individual probabilities range by five orders of magnitude from almost $0.000\,08$ to less than $0.000\,000\,000\,1.$ Nevertheless, this sample was produced in a fraction of a second using the general function s posted at https://stats.stackexchange.com/a/606853/919. That algorithm scales well: generating a larger sample of $10^{30}$ values required only one more second. The computation time depends primarily on the expected range of the sample, not on the sample size.
The implementation challenge is twofold: first, to overcome the lack of a closed form for the normalizing constant; second, to handle potential numerical overflow. The raw values given by the formula are largest at the mode of $6\,250\,000,$ where they get as large as
$$\Pr(6250000;50,1/4)\ \propto\ \frac{50^{6250000}}{(6250000!)^{1/4}} \approx 10^{678584.2}.$$
That will overflow your calculator ;-).
The solution (coded in R as the dCMPZ and dCMP functions) given below overcomes these difficulties by searching from the mode (which occurs near $\lceil \lambda ^ {1/\nu}\rceil$) outwards into each tail until the individual probabilities are tiny compared to the probability at the mode. These are taken to be the entire support of the distribution, essentially truncating it at lower and upper limits beyond which your simulation is unlikely to produce any values. To overcome the numerical overflow problems, these probabilities are first divided by the modal probability (the largest possible one) and then summed up -- and all this is done on a logarithmic scale.
Finally, if you need an iid sequence of values for the simulation, rather than a tabulation of the frequencies of its values, it's straightforward to convert this tabulated output into the individual values with their multiplicities and randomly permute them. Just don't try this with a sample of size $10^{30}$!
#
# Un-normalized CMP probability at `x`.
# Requires lambda > 0, nu >= 0.
#
dCMP <- function(x, lambda = 1, nu = 1, log.p = FALSE) {
q <- x * log(lambda) - nu * lfactorial(x)
if(!isTRUE(log.p)) q <- exp(q)
q
}
#
# Return essentially all the probabilities, ignoring probabilities less than `tol`
# compared to the mode. `n.max` assures the function won't hang up when the
# search is too long: it limits the range of the search.
#
dCMPZ <- function(lambda = 1, nu = 1, tol = 1e-16, n.max = 1e6) {
m <- ceiling(lambda ^ (1/nu)) # A mode + 1
p.max <- dCMP(m, lambda, nu, log.p = TRUE)
if (is.infinite(p.max)) stop("CMP parameters cause logarithmic overflow.")
#
# Search into the upper tail.
# It proceeds in blocks of values of length `stride` to capitalize on
# vectorized computation. These blocks are (asymptotically in lambda)
# one SD long.
#
v <- lambda^(1/nu) / nu * (1 + (nu^2 - 1)/(24 * nu^2) * lambda^(-2/nu))
if (is.na(v) | is.infinite(v) | v <= 0) v <- 1
stride <- min(ceiling(n.max/2), m + 1, ceiling(sqrt(v)))
m.start <- m - 1
x <- m.start + seq_len(stride)
p <- p1 <- dCMP(x, lambda, nu, log.p = TRUE)
p.min <- p.max + log(tol)
while(isTRUE(p1[stride] > p.min) && isTRUE(length(x) < n.max)) {
m.start <- m.start + stride
x1 <- m.start + seq_len(stride)
p1 <- dCMP(x1, lambda, nu, log.p = TRUE)
x <- c(x, x1)
p <- c(p, p1)
}
if (length(x) >= n.max)
warning("The problem is too large; a partial distribution is returned.")
#
# Compute the lower tail.
# This assumes all the CMP distributions have positive skewness, implying
# it suffices to search no further out into the lower tail than into the
# upper tail.
#
x0 <- seq(max(0, m - x[length(x)] + x[1]), x[1] - 1)
x <- c(x0, x)
p <- c(dCMP(x0, lambda, nu, log.p = TRUE), p) - p.max
#
# Order by increasing probability for maximum precision in the summation.
#
i <- order(p, decreasing = FALSE)
x <- x[i]
p <- exp(p[i]) # Convert from logs to relative probabilities (max is 1).
p <- p / sum(p) # The normalized probability distribution.
#
# Discard probabilities reduced to zero by underflow in the normalization.
#
j <- which(p > 0)
x <- x[j]
p <- p[j]
#
# Return the values, their probabilities, and a record of the caller's
# arguments to the function.
#
return(list(x = x, probs = p, tol = tol, lambda = lambda, nu = nu))
}
#
# Examples of dCMPZ.
#
theta <- list(c(1e-2, 1e-4),
c(1, 2),
c(5, 1/2),
c(1e6, 3))
pars <- par(mfrow = c(2,2))
for (q in theta)
with(dCMPZ(q[1], q[2], tol = 1e-6),
plot(x, probs, type = "h", col = hsv(0.01, 1, 0.7),
xlab = "Value", ylab = "Probability",
main = bquote(paste(lambda == .(q[1]), ", ", nu == .(q[2])))))
par(pars)
#==============================================================================#
#
# Sample from a discrete distribution on the values 1, 2, 3, ..., length(p).
# Returns the *unordered* sample as a vector tabulation of the counts.
#
s <- function(m, p) {
# Initialization.
P <- rev(cumsum(rev(p)))
k <- rep(0, length(p))
# The algorithm.
for (i in seq_along(p)) {
k[i] <- rbinom(1, m, p[i]/P[i])
m <- m - k[i]
}
k
}
#==============================================================================#
#
# Simulating from a CMP distribution.
#
set.seed(17) # For reproducibility
n <- 1e6 # Sample size
lambda <- 50; nu <- 1/4
system.time({
obj <- dCMPZ(lambda, nu) # Generate the full distribution
k <- s(n, obj$probs) # Sample from it and summarize the sample
})
i <- k > 0 # Remove zero counts
plot(obj$x[i], k[i], col = gray(.9), xlab = "Value", ylab = "Count",
main = bquote(paste("CMP Simulation for ", lambda == .(lambda),
" and ", nu == .(nu))))
j <- order(obj$x) # Need to place values in order for drawing a curve!
lines(obj$x[j], n * obj$probs[j], lwd = 2, col = hsv(0.01, 1, .9))
#------------------------------------------------------------------------------#
#
# A quick test that the sample algorithm `s` is correct.
# The distribution is broken approximately into equal quantiles, then a sample
# is summarized by those pre-specified bins and compared to the distribution
# with a chi-squared test.
#
n.bins <- min(10, length(obj$x)) # Specify how many bins to use
P <- cumsum(obj$probs)
bins <- ceiling(P * n.bins) # Assign a bin to each possible value
probs <- sapply(by(obj$probs, bins, `+`), sum)
n <- n.bins * 5e3 # Sample size -- should be at least five times the number of bins
samp <- s(n, obj$probs)
counts <- sapply(by(samp, bins, `+`), sum) # Count by bin
chisq.test(counts, p = probs)
#==============================================================================#
#
# How to generate an ordered *iid* sample.
#
obj <- dCMPZ(5, 0.7) # Precompute the distribution
k <- s(5000, obj$probs) # Sample from it
x <- sample(rep(obj$x, k)) # Replicate each value and randomly permute the vector
# -- Display a histogram of the sample to check.
hist(x, freq = FALSE, breaks = seq(0, length(k)+1) - 1/2)
j <- order(obj$x)
lines(obj$x[j], obj$probs[j], lwd = 2, col = hsv(0.01, 1, .9))
