3
$\begingroup$

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.

  1. Is there a method that I can use to simulate from this distribution?

  2. Is there any R package that implements this simulation?

$\endgroup$
11
  • 2
    $\begingroup$ Sure, even naively trying Metropolis-Hastings would probably work. The rest of your question looks like a software request, which is not on topic per se. You can calculate an approximation of $Z(\lambda, \nu)$ rather than using an exact value. Your computer only has finite numerical precision, even for a given choice of 'arbitrary' floating point precision. $\endgroup$
    – Galen
    Mar 7, 2023 at 16:35
  • 1
    $\begingroup$ Welcome to CV. Do you need to simulate data from a fixed distribution repeatedly, or will the parameters change frequently? What are the expected ranges of the parameters? How much data do you need to simulate? Such issues are critical to the specification, design, and recommendation of almost any software. $\endgroup$
    – whuber
    Mar 7, 2023 at 16:39
  • 1
    $\begingroup$ Metropolis-Hastings is an algorithm. I wouldn't say any algorithm is automatic per se, but rather a software implementation of the algorithm automates performing the algorithm. $\endgroup$
    – Galen
    Mar 7, 2023 at 17:10
  • 1
    $\begingroup$ For the simulation, since you don't intend to explorate "extreme or pathological" aspects of the distribution, there won't be too many distinct values in any simulation -- shall we suppose less than a million of them? Use the algorithm (with code) I recently posted at stats.stackexchange.com/a/606853/919 (illustrated with a Poisson distribution, but described for any discrete distribution). The calculation of the individual probabilities is quick and easy. $\endgroup$
    – whuber
    Mar 7, 2023 at 19:41
  • 1
    $\begingroup$ @Xi'an Thanks. I think an inversion method is what I can use/implement. Case closed :). $\endgroup$
    – CMP
    Mar 7, 2023 at 19:43

1 Answer 1

5
$\begingroup$

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):

enter image description here

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:$

enter image description here

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))
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.