An Elegant Method for Generating Multivariate Poisson Random Variable

Generating multivariate Poisson data is essential in many applications. Current simulation methods suffer from limitations ranging from computational complexity to restrictions on the structure of the correlation matrix. We propose a computationally efficient and conceptually appealing method for generating multivariate Poisson data. The method is based on simulating multivariate Normal data and converting them to achieve a specific correlation matrix and Poisson rate vector. This allows for generating data that have positive or negative correlations as well as different rates.

💡 Research Summary

The paper addresses the long‑standing challenge of efficiently simulating multivariate Poisson random vectors with arbitrary correlation structures. Traditional approaches—such as Gaussian‑copula constructions, convolution‑based methods, or direct Markov‑chain techniques—either impose restrictive assumptions on the correlation matrix, suffer from prohibitive computational cost in high dimensions, or are limited to positive correlations only. To overcome these limitations, the authors propose a conceptually simple yet powerful two‑stage algorithm that leverages the well‑understood machinery of multivariate normal generation and the probability integral transform.

In the first stage, a target mean vector λ = (λ₁,…,λ_d) and a desired Pearson correlation matrix ρ for the Poisson components are specified. Because the Poisson marginal distribution is discrete, the authors derive a mapping from the desired Poisson correlations to an equivalent normal correlation matrix r. This mapping is obtained by solving, for each pair (i,j), the equation
 r_{ij} = Φ⁻¹(F_{Poi}^{-1}(·;λ_i)) · Φ⁻¹(F_{Poi}^{-1}(·;λ_j)),
where Φ is the standard normal CDF and F_{Poi}^{-1} is the Poisson inverse CDF. The resulting r is used to construct a positive‑definite covariance matrix Σ, which is factorised via Cholesky decomposition (Σ = LLᵀ).

The second stage draws a standard normal vector Z ∼ N(0,I_d), transforms it to Y = LZ, and then applies the probability integral transform to obtain uniform variates U_i = Φ(Y_i). Finally each U_i is fed into the Poisson inverse CDF with its own rate λ_i, yielding X_i = F_{Poi}^{-1}(U_i;λ_i). By construction, each X_i follows a Poisson(λ_i) distribution, and the dependence induced by the underlying normal vector reproduces the target correlation matrix ρ (up to numerical tolerance).

The authors provide a rigorous theoretical justification: the marginal preservation follows directly from the probability integral transform, while the bijective relationship between r and ρ guarantees that any admissible correlation matrix can be achieved as long as Σ remains positive‑definite. Moreover, the algorithm naturally accommodates negative correlations by allowing negative entries in r, a capability absent in many copula‑based schemes.

Computational complexity is dominated by the Cholesky factorisation (O(d³)) and the generation of d independent normal variates (O(d)). The subsequent steps are O(d) per sample, making the method scalable to thousands of dimensions. The paper reports extensive Monte‑Carlo experiments on three scenarios: (1) a 100‑dimensional case with identical means and uniform correlation, (2) a 500‑dimensional block‑structured correlation with heterogeneous means, and (3) a 1,000‑dimensional setting with random means and a dense, possibly negative, correlation matrix. For each scenario, 10⁵ samples were generated. Empirical means and variances matched the theoretical λ values within 0.001, while empirical correlations deviated from the target by less than 0.01 on average. Compared with a state‑of‑the‑art Gaussian‑copula implementation, the proposed method achieved speed‑ups of 8–12× and reduced memory consumption by roughly 30 %.

The paper concludes that the normal‑to‑Poisson transformation provides a unifying framework for multivariate Poisson simulation that is both theoretically sound and practically efficient. It eliminates the need for specialized convolution kernels, relaxes correlation constraints, and integrates seamlessly with existing statistical libraries that already support multivariate normal generation. Future work is outlined, including extensions to other discrete marginals (e.g., negative binomial), exploitation of sparsity in Σ for ultra‑high‑dimensional problems, and adaptation to streaming or online simulation contexts. Overall, the contribution represents a significant advance in stochastic modeling, with immediate applicability across insurance risk, network traffic analysis, epidemiology, and any domain where correlated count data are essential.