如何在Python中对任意复杂多元函数进行高效采样？

Alright, let's figure out how to sample from your complex multivariate function since inverse transform sampling (which only works for univariate cases) isn't going to cut it here. Here are some practical, battle-tested approaches that work for this scenario:

1. Markov Chain Monte Carlo (MCMC) Methods #

This is the go-to category for most complex multivariate sampling problems, especially when you can't compute the inverse CDF or even don't know the normalization constant of your target distribution (you just need to calculate relative probability densities).

• Metropolis-Hastings Algorithm:
  It's a foundational MCMC method that's straightforward to use (either via libraries or custom implementation). The core idea is to build a Markov chain that converges to your target distribution. The basic steps are:

  1. Start with an initial sample point $\boldsymbol{x}_0$
  2. Generate a candidate sample $\boldsymbol{x}^*$ from a proposal distribution (like a multivariate normal centered at the current sample)
  3. Calculate the acceptance rate $\alpha = \min\left(1, \frac{p(\boldsymbol{x}^)}{p(\boldsymbol{x}_t)} \cdot \frac{q(\boldsymbol{x}_t|\boldsymbol{x}^{*)}{q(\boldsymbol{x}}|\boldsymbol{x}_t)}\right)$ — here $p$ is your target distribution, and $q$ is the proposal distribution
  4. Draw a uniform random number $u \sim U(0,1)$; if $u \leq \alpha$, accept $\boldsymbol{x}^*$ as the next sample, otherwise keep $\boldsymbol{x}_t$
  5. Repeat steps 2-4 until the chain stabilizes (converges), then use subsequent samples as draws from your target distribution

• Hamiltonian Monte Carlo (HMC) / No-U-Turn Sampler (NUTS):
  If you can compute the gradient of your target distribution's log probability, HMC is way more efficient than Metropolis-Hastings. It uses Hamiltonian dynamics to generate candidate samples that are much more likely to be accepted, reducing the time needed for convergence. NUTS is an adaptive variant of HMC that automatically tunes its parameters, making it even easier to use. Most modern probabilistic programming libraries have built-in NUTS support.

2. Multivariate Rejection Sampling (and Variants) #

Rejection sampling works by using a simpler, easy-to-sample proposal distribution $g(\boldsymbol{x})$ that "covers" your target distribution $p(\boldsymbol{x})$ (i.e., there exists a constant $M$ where $p(\boldsymbol{x}) \leq M \cdot g(\boldsymbol{x})$ for all $\boldsymbol{x}$). The steps are:

1. Sample a candidate $\boldsymbol{x}^*$ from $g(\boldsymbol{x})$
2. Sample $u$ from $U(0, M \cdot g(\boldsymbol{x}^*))$
3. If $u \leq p(\boldsymbol{x}^)$, accept $\boldsymbol{x}^$; otherwise, reject it and try again

The catch here is that for high-dimensional or very complex distributions, finding a good $g(\boldsymbol{x})$ and $M$ can be tricky (it might lead to extremely low acceptance rates). Adaptive rejection sampling has multivariate extensions, but they're often more complex to implement than MCMC methods.

3. Importance Sampling (for Expectation Calculations) #

If your end goal is to compute expectations (rather than get direct samples from the distribution), importance sampling is a solid option. You sample from a proposal distribution $g(\boldsymbol{x})$, then assign each sample a weight $w(\boldsymbol{x}) = \frac{p(\boldsymbol{x})}{g(\boldsymbol{x})}$. You can then use these weighted samples to estimate expectations.

If you need approximate independent samples, you can follow up with importance resampling: draw samples from your importance sample pool, with replacement, using the weights as probabilities. Note that this works best when the weights don't vary too widely.

4. Use Off-the-Shelf Python Libraries #

You don't need to implement these algorithms from scratch! There are great Python libraries that handle the heavy lifting:

• PyMC3: A probabilistic programming library that supports NUTS, Metropolis-Hastings, and more. You just define your target distribution's log probability (even as a custom function) and let the library handle sampling.
• emcee: A lightweight library focused on MCMC, perfect for custom Metropolis-Hastings implementations if you need more control.
• PyStan: The Python interface for Stan, which has an industry-leading implementation of the NUTS algorithm, ideal for complex multivariate models.

Here's a quick example using PyMC3 for a custom multivariate distribution:

import pymc3 as pm
import numpy as np

# Define the log probability density of your custom multivariate distribution
def log_target(x):
    # Replace this with your actual complex function
    # Example: a custom bivariate distribution
    return -np.sum(x**2) - np.exp(x[0] * x[1])

with pm.Model() as model:
    # Define a multivariate variable with a uniform initial proposal range
    x = pm.Uniform("x", lower=-5, upper=5, shape=2)
    # Attach our custom log probability
    pm.Potential("log_p", log_target(x))
    # Run NUTS sampling (tune phase warms up the chain)
    trace = pm.sample(2000, tune=1000, cores=2)

# Extract the sampled points
samples = trace["x"]

Final Takeaway #

For complex multivariate functions where inverse transform sampling isn't feasible, MCMC methods (especially NUTS if you can compute gradients) are the most practical choice. Using probabilistic programming libraries will save you tons of time and avoid the pitfalls of implementing MCMC correctly from scratch.

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