关于使用Metropolis-Hastings算法对Conway-Maxwell-Poisson（CMP）分布进行采样的技术咨询

Hey Pedro, totally get where you're stuck here—CMP’s a tricky distribution without a closed-form CDF, so ditching the inversion method makes total sense. Let’s walk through exactly how to set up your Metropolis-Hastings (MH) sampler with the Poisson proposal you’ve already picked, plus some practical tips to keep things running smoothly.

核心思路：MH算法在离散CMP分布上的适配 #

First, let’s recap the key bits you already have right: your target CMP distribution has a kernel $ \pi(x) \propto \lambda^x / (x!)^\nu $, and using a Poisson proposal is a smart call because it matches the "shape" of the CMP kernel closely (especially since CMP reduces to Poisson when $\nu=1$). For discrete distributions, the MH acceptance probability formula is straightforward—we just need to compare the target distribution and proposal distribution ratios correctly.

步骤1：定义提议分布与接受概率 #

Let’s formalize this:

• Target kernel: $\pi(x) \propto \lambda^x / (x!)^\nu$
• Proposal distribution: We’ll use a Poisson distribution centered at the current state $x$, so $Q(y|x) = \text{Poisson}(y; x) = e^{-x} x^y / y!$ (this is a symmetric random-walk style proposal, which is stable for discrete counts)

The acceptance probability $\alpha(x,y)$ when proposing $y$ from current state $x$ is:
$$
\alpha(x,y) = \min\left(1, \frac{\pi(y) Q(x|y)}{\pi(x) Q(y|x)}\right)
$$

We can simplify this ratio to avoid messy calculations and numerical overflow:

1. First, compute the target ratio: $\frac{\pi(y)}{\pi(x)} = \lambda^{y-x} \cdot \left( \frac{x!}{y!} \right)^\nu$
2. Then the proposal ratio (since we're using Poisson($x$) for $Q(y|x)$ and Poisson($y$) for $Q(x|y)$): $\frac{Q(x|y)}{Q(y|x)} = e^{x-y} \cdot \frac{y^x y!}{x^y x!}$

Multiply these together and simplify to get a more computable form:
$$
\frac{\pi(y) Q(x|y)}{\pi(x) Q(y|x)} = \left( \frac{\lambda}{e} \right)^{y-x} \cdot \left( \frac{x!}{y!} \right)^{\nu-1} \cdot \frac{y^x}{xy}
$$

步骤2：用对数计算避免数值问题 #

Directly computing factorials and large powers will cause numerical underflow/overflow fast, so always work in log space for these calculations. Taking the natural log of the ratio gives:
$$
\log(r) = (y-x)\log\left(\frac{\lambda}{e}\right) + (\nu-1)\left(\log(x!) - \log(y!)\right) + x\log(y) - y\log(x)
$$

Then exponentiate this value to get the ratio $r$, and set $\alpha = \min(1, r)$. If $\log(r) > 0$, $\alpha$ is just 1 (we always accept the proposal).

步骤3：完整的MH采样流程 #

Here’s a step-by-step breakdown of the sampler:

• Initialize: Pick an initial state $x_0$ (a good starting point is the mode of the CMP, which is roughly $\lfloor \lambda^{1/\nu} \rfloor$—this speeds up convergence)
• Burn-in & Iteration:
  1. For each step from 1 to $N + \text{burn_in}$:
     • Sample a candidate $y$ from $\text{Poisson}(x_{\text{current}})$
     • Compute $\log(r)$ using the formula above (handle $x=0$ or $y=0$ by adding a tiny epsilon like $10^{-10}$ to avoid $\log(0)$ errors)
     • Calculate $\alpha = \min(1, \exp(\log(r)))$
     • Generate a uniform random number $u \sim U(0,1)$:
       • If $u \leq \alpha$, set $x_{\text{current}} = y$
       • Else, keep $x_{\text{current}}$ as is
  2. Discard burn-in samples: The first $M$ samples (e.g., 10% of total, or use convergence diagnostics to pick $M$) are discarded because the sampler needs time to stabilize to the target distribution. The remaining samples are your CMP draws.

实用技巧与注意事项 #

• Log factorial calculation: Use the gamma log function (e.g., scipy.special.gammaln in Python) since $\log(x!) = \text{gammaln}(x+1)$—this is fast and numerically stable.
• Convergence checks: Use trace plots (plot the sequence of $x_{\text{current}}$ over time) to ensure the sampler has stabilized. You can also use the Gelman-Rubin diagnostic if running multiple chains.
• Adjusting the proposal: If you’re working with extreme $\nu$ values (far from 1), you could tweak the proposal’s mean (e.g., use $\mu = x \cdot c$ for some constant $c$) to balance acceptance rate and mixing, but the Poisson($x$) proposal works well for most cases.

示例伪代码（Python风格） #

import numpy as np
from scipy.special import gammaln

def sample_cmp_mh(lambda_, nu, n_samples, burn_in=1000):
    # 初始化：用CMP的众数作为起点，加速收敛
    x_current = int(np.floor(lambda_ ** (1/nu)))
    samples = []
    
    for step in range(n_samples + burn_in):
        # 从Poisson(x_current)采样候选y
        y = np.random.poisson(x_current)
        
        # 计算对数形式的比值，避免数值问题
        # 处理x=0或y=0的情况，防止log(0)
        log_x = np.log(x_current + 1e-10)
        log_y = np.log(y + 1e-10)
        log_fact_x = gammaln(x_current + 1)
        log_fact_y = gammaln(y + 1)
        
        log_ratio = (y - x_current) * np.log(lambda_ / np.e)
        log_ratio += (nu - 1) * (log_fact_x - log_fact_y)
        log_ratio += x_current * log_y - y * log_x
        
        # 计算接受概率
        alpha = min(1.0, np.exp(log_ratio))
        
        # 接受/拒绝
        if np.random.uniform() < alpha:
            x_current = y
        
        # 收集样本（跳过burn-in）
        if step >= burn_in:
            samples.append(x_current)
    
    return np.array(samples)

备注：内容来源于stack exchange，提问作者Pedro Lemes