如何从任意概率密度函数生成IID样本？Python实现正确性验证

Great question! Let's break down your implementation, talk about what works, where it falls short for arbitrary PDFs, and how to make it more robust.

First: Does your current code work for the case you tested? #

Yes! For the PDF pp(x) = x² over the interval [0, 1] (your default M=1), your approach does generate samples that roughly follow the target distribution. Here's why it works for this specific case:

• You sample N points uniformly from [0, M]
• You weight each point by its PDF value f(x), normalize those weights to sum to 1, then re-sample from the original points using those weights.
• Since uniform sampling has a constant density across [0,M], normalizing the f(x) values gives a valid discrete approximation of the continuous PDF's relative probabilities.

Key limitations that prevent it from handling any arbitrary PDF: #

Your implementation has hardcoded assumptions that break for many valid PDFs:

• Fixed sampling interval: Your function uses [0, M] as the sampling range, but arbitrary PDFs can have any domain (e.g., (-∞, ∞) for normal distributions, [2, 7] for a bounded non-zero PDF). If you don't match the sampling interval to the actual region where f(x) is non-zero, you'll either miss critical parts of the distribution (leading to biased samples) or waste computation on regions where f(x) is zero.
• Discrete approximation error: You're using N points to approximate a continuous distribution. For PDFs with sharp peaks or rapid changes, this can lead to poor accuracy unless you use far more initial candidate points than your final sample size N. For example, if f(x) has a narrow peak that only 1% of your uniform samples land on, your weighted sampling won't capture that peak properly.
• No PDF validation: A valid PDF must be non-negative everywhere on its domain and have a finite total integral. Your function doesn't check for these conditions—if someone passes a negative-valued function or a function that doesn't integrate to a finite value, the code will either produce invalid negative probabilities or nonsensical samples.
• Inefficiency for sparse distributions: If the PDF is concentrated in a tiny subset of the sampling interval, most of your initial uniform samples will have near-zero weight, leading to wasted computation and worse sampling quality.

How to improve this for arbitrary PDFs: #

Here are actionable fixes to make your function more robust:

1. Let users specify the PDF's valid domain: Modify the function to accept x_min and x_max (for bounded domains) or add logic to handle unbounded domains via practical truncation. For example:

   import numpy as np
   import matplotlib.pyplot as plt
   
   def randgen(f, N, x_min, x_max, num_candidates=10*N):
       # Use more candidates to reduce discrete approximation error
       sample_candidates = np.random.uniform(x_min, x_max, num_candidates)
       weights = f(sample_candidates)
       
       # Validate the input PDF is non-negative
       if np.any(weights < 0):
           raise ValueError("PDF must return non-negative values everywhere on its domain")
       
       # Normalize weights to form a valid probability distribution
       weights /= np.sum(weights)
       
       # Sample from candidates using the weighted probabilities
       return np.random.choice(sample_candidates, p=weights, size=N)
   

2. Use more candidate points: Generating 10x or 100x more initial samples than your final N reduces discrete approximation error and improves sampling accuracy, especially for complex PDFs.
3. Add basic validation: Checking for non-negative weights catches invalid PDF inputs early. You could also add a warning if the sum of weights is very small (hinting the sampling interval might not match the PDF's actual support).
4. Consider advanced methods for tricky PDFs: For unbounded domains or PDFs with complex shapes, methods like inverse transform sampling (if you can compute the CDF inverse), accept-reject sampling, or Metropolis-Hastings (a Markov Chain Monte Carlo method) are more reliable. Your current approach is a simple discrete approximation that works well for bounded, smooth PDFs but isn't universal.

Example test with a normal distribution: #

If you test with a truncated normal PDF, you'll see the improved function produces samples that closely match the true distribution:

def normal_pdf(x):
    return (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x**2)

# Generate samples from truncated normal (-3 to 3)
z = randgen(normal_pdf, 2000, x_min=-3, x_max=3)

# Plot histogram and overlay true PDF
plt.hist(z, bins=30, density=True, alpha=0.7)
x = np.linspace(-3, 3, 1000)
plt.plot(x, normal_pdf(x), 'r', linewidth=2)
plt.title("Sampled Normal Distribution vs True PDF")
plt.show()

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