如何在MATLAB中手动生成Conway-Maxwell-Poisson（COM-Poisson）分布的随机数？

Since you're building your entire simulation pipeline in MATLAB and want to avoid switching to R packages, let's implement a custom solution for COM-Poisson random number generation using inverse transform sampling—this is a reliable, straightforward approach for discrete distributions like COM-Poisson.

Quick Recap of the COM-Poisson PMF #

First, let's restate the PMF you referenced from Guikema & Goffelt (2008) for clarity:

P(Y=y) = (1/S(mu, nu)) * (mu^y / y!)^nu, where S(mu, nu) = ∑ₙ=₀^∞ (muⁿ / n!)^nu

The main hurdles here are safely calculating the infinite series normalization constant S(mu, nu) and building a cumulative distribution function (CDF) for sampling. Let's break this down step by step.

Step-by-Step Implementation #

1. Reusable MATLAB Function

Here's a custom function that handles all the heavy lifting:

function r = compoissrnd(mu, nu, n_samples)
    % COMPOISSRND Generate COM-Poisson random numbers
    % Inputs:
    %   mu        - Mean parameter (positive scalar)
    %   nu        - Dispersion parameter (positive scalar)
    %   n_samples - Number of random samples to generate
    % Output:
    %   r         - Vector of COM-Poisson random samples
    
    % Step 1: Calculate normalization constant S(mu, nu) via truncated series
    tol = 1e-12; % Stop adding terms once they're negligible
    S = 0;
    term = 1; % First term (n=0): (mu^0/0!)^nu = 1^nu = 1
    n = 0;
    
    while term > tol
        S = S + term;
        n = n + 1;
        % Use log-transforms to avoid numerical overflow
        log_term = nu * (n*log(mu) - sum(log(1:n)));
        term = exp(log_term);
    end
    
    % Step 2: Build the CDF from PMF values
    cdf_vals = zeros(1, n);
    cdf_vals(1) = 1/S; % PMF for y=0
    for y = 1:n-1
        log_pmf = nu * (y*log(mu) - sum(log(1:y)));
        pmf_val = exp(log_pmf)/S;
        cdf_vals(y+1) = cdf_vals(y) + pmf_val;
    end
    
    % Step 3: Inverse transform sampling
    u = rand(1, n_samples);
    r = zeros(1, n_samples);
    
    for i = 1:n_samples
        % Find the smallest index where CDF >= the uniform random value
        r(i) = find(cdf_vals >= u(i), 1) - 1; % Adjust to 0-based y values
    end
end

2. Key Details to Note

• Numerical Stability: Using log-transforms for PMF calculations prevents overflow, which is crucial when working with large values of mu or y.
• Truncation Tolerance: The tol = 1e-12 ensures we only stop adding terms to S(mu, nu) when they're mathematically negligible, keeping the approximation accurate for simulation.
• Validation Check: When nu = 1, the COM-Poisson distribution reduces to the standard Poisson distribution. Test this by comparing samples from compoissrnd(mu, 1, 1000) with MATLAB's built-in poissrnd(mu, 1, 1000)—they should match closely.

Example Usage #

% Generate 1000 COM-Poisson samples with mu=3, nu=0.5
samples = compoissrnd(3, 0.5, 1000);

% Visualize the sample distribution
histogram(samples, 'Normalization', 'probability');
xlabel('Y Value');
ylabel('Empirical Probability');
title('COM-Poisson Distribution (mu=3, nu=0.5)');

This will produce a histogram that aligns with the theoretical PMF of your chosen COM-Poisson parameters.

内容的提问来源于stack exchange，提问作者UlduzM