如何用PyMC3从无解析形式的数据集分布函数生成样本？

核心思路 #

PyMC3没有直接调用自定义密度函数生成样本的快捷方式，但可以通过自定义对数概率函数结合无梯度采样器（如Metropolis）实现。你之前的问题大概率是误用了pm.Densityhist（这是绘图工具，不是分布定义类），且默认的NUTS采样器需要梯度信息，而非解析分布无法提供，导致采样失败。

一维场景实现步骤 #

1. 准备对数概率函数

先将密度函数转换为对数概率函数（避免数值下溢），同时处理密度为0的边界情况：

import pymc3 as pm
import numpy as np
from scipy.stats import gaussian_kde

# 示例：用核密度估计模拟你从数据得到的非解析PDF
# 替换成你自己的真实密度函数即可
raw_data = np.concatenate([np.random.normal(-2, 1, 500), np.random.normal(2, 1, 500)])
custom_pdf = gaussian_kde(raw_data)

def logp(x):
    pdf_val = custom_pdf(x)
    # 避免log(0)报错，密度极小时返回负无穷
    return np.log(pdf_val) if pdf_val > 1e-10 else -np.inf

2. 构建PyMC3模型并采样

用pm.DensityDist定义自定义分布，指定Metropolis采样器（无需梯度）：

with pm.Model() as model:
    # 定义自定义分布，testval设为数据均值作为初始采样点
    x = pm.DensityDist("x", logp, testval=np.mean(raw_data))
    
    # 采样：tune期用于调整采样器参数，之后取1000个有效样本
    trace = pm.sample(
        2000,          # 总采样数
        tune=1000,     # 预热期（丢弃这部分样本）
        cores=2,
        step=pm.Metropolis(),  # 指定无梯度采样器
        return_inferencedata=False  # 兼容旧版PyMC3的trace格式
    )

# 提取1000个有效样本
samples = trace["x"][1000:]  # 跳过预热期，刚好取1000个

3. 验证样本分布

通过直方图和原PDF对比验证样本质量：

import matplotlib.pyplot as plt

x_grid = np.linspace(-5, 5, 1000)
plt.plot(x_grid, custom_pdf(x_grid), label="Original PDF")
plt.hist(samples, bins=30, density=True, alpha=0.5, label="Generated Samples")
plt.legend()
plt.show()

高维场景扩展 #

高维场景逻辑和一维一致，只需调整两点：

1. 对数概率函数：接受高维输入（比如二维输入是x = [x1, x2]），确保你的密度函数支持高维计算。
2. 初始值：testval设为高维数据的均值向量（比如二维的np.mean(raw_data, axis=0)）。

注意：高维下Metropolis采样效率会下降，可尝试pm.Slice采样器，或延长预热期（如tune=2000）、增加总采样数来保证样本质量。

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