如何利用蒙特卡洛模拟生成给定解析形式的概率密度函数（PDF）

首先，你给出的概率密度函数（PDF）是Gamma分布的标准形式，其中：

• 形状参数为 d
• 尺度参数为 g

生成样本的常用方法 #

方法1：利用编程语言内置工具

大多数科学计算库都内置了Gamma分布的随机数生成函数，这是最便捷的实现方式：

• Python：使用 numpy.random.gamma()，参数对应 shape=d、scale=g
• R：使用 rgamma(n, shape=d, scale=g)

方法2：手动实现（无内置函数时）

如果需要手动编写逻辑，可根据参数类型选择方案：

• 当d为整数时：Gamma分布等价于d个独立指数分布的和，指数分布的尺度参数为 g。只需生成d个服从 Exp(g) 的样本并求和，即可得到一个Gamma分布样本。
• 当d为非整数时：推荐使用Marsaglia-Tsang算法（高效通用），或接受-拒绝采样（需选择合适的提议分布，如正态分布，计算接受概率筛选样本）。

代码示例（Python） #

1. 内置函数生成样本并验证PDF

import numpy as np
import math
import matplotlib.pyplot as plt

# 定义参数
d = 2.5  # 形状参数
g = 1.2  # 尺度参数

# 生成蒙特卡洛样本
sample_size = 10000
samples = np.random.gamma(shape=d, scale=g, size=sample_size)

# 计算解析PDF
def gamma_pdf(x, d, g):
    return (x**(d-1) * math.exp(-x/g)) / (math.gamma(d) * (g**d))

# 可视化对比样本直方图与解析PDF
x = np.linspace(0, np.max(samples), 100)
pdf_values = [gamma_pdf(xi, d, g) for xi in x]

plt.hist(samples, bins=50, density=True, alpha=0.6, label='蒙特卡洛样本直方图')
plt.plot(x, pdf_values, 'r-', label='解析PDF')
plt.xlabel('x')
plt.ylabel('概率密度')
plt.legend()
plt.show()

2. 手动实现（整数d场景）

import numpy as np

d = 3  # 整数形状参数
g = 1.5
sample_size = 10000

# 生成d个独立指数分布样本并求和
exponential_samples = np.random.exponential(scale=g, size=(sample_size, d))
gamma_samples = np.sum(exponential_samples, axis=1)

核心说明 #

蒙特卡洛模拟的目标是生成服从目标分布的随机样本，而非直接计算PDF——你之前通过代入x计算PDF的操作，是用于验证样本是否符合分布的手段。对于Gamma分布，优先使用内置函数或成熟算法，避免重复造轮子。

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