Python中如何根据给定连续PDF生成单个服从分布的随机数

这个单变量连续分布没有对应numpy、scipy的内置采样接口，直接用逆变换采样实现即可，步骤如下：

1. 前置推导：归一化常数与CDF #

首先化简PDF的被积项：
$$\frac{1+(1-z)^2}{z} = z - 2 + \frac{2}{z}$$
z的定义域为$[z_0, 1]$，根据归一化条件$\int_{z_0}^1 p(z)dz=1$，可直接算出归一化常数N：
不定积分结果为$\int (z-2+2/z)dz = 0.5z^2 -2z + 2\ln z + C$，代入上下限可得：
$$N = \frac{1}{-0.5z_0^2 + 2z_0 - 2\ln z_0 - 1.5}$$
对应累积分布函数（CDF，即随机变量取值小于等于z的概率）为：
$$F(z) = N\cdot\left[(0.5z^2 -2z + 2\ln z) - (0.5z_0^2 -2z_0 + 2\ln z_0)\right]$$
由于p(z)在定义域内恒正，F(z)在$[z_0,1]$上严格单调递增，反函数存在。

2. 采样逻辑 #

逆变换采样的核心规则：若u服从[0,1]区间的均匀分布，则$F^{-1}(u)$服从目标分布。
由于F(z)的反函数没有解析闭式，直接用有界区间数值求根方法求解$F(z)=u$即可，效率和精度都能满足需求。

3. 可直接复用的代码实现 #

import numpy as np
from scipy.optimize import root_scalar

def build_z_sampler(z0):
    # 预计算固定常量，避免采样时重复计算
    integral_const = 0.5 * z0**2 - 2 * z0 + 2 * np.log(z0)
    norm_N = 1 / ((0.5 - 2) - integral_const)  # 代入z=1时的积分值

    def _cdf(z):
        return norm_N * (0.5 * z**2 - 2*z + 2*np.log(z) - integral_const)

    def sample(count=1):
        uniform_rands = np.random.uniform(0, 1, size=count)
        results = []
        for u in uniform_rands:
            # 用Brent法在[z0,1]区间求根，收敛速度快、精度稳定
            root_res = root_scalar(lambda z: _cdf(z) - u, bracket=[z0, 1], method='brentq')
            results.append(root_res.root)
        return results[0] if count == 1 else np.array(results)
    
    return sample

# 用法示例
if __name__ == "__main__":
    z0 = 0.2  # 替换为实际的z0取值
    z_sampler = build_z_sampler(z0)
    # 生成单个符合分布的随机数
    single_num = z_sampler()
    print(f"单个随机数结果：{single_num:.4f}")
    # 批量生成10000个随机数可用于分布校验
    # batch_nums = z_sampler(10000)

注意事项 #

• 上述实现生成单个随机数的耗时在微秒级，常规使用场景下完全不需要做性能优化。
• 若z0取值极小（小于1e-10），$2\ln z_0$项可能出现浮点数下溢，常规取值范围下无精度问题。
• 不需要强行匹配numpy/scipy的内置分布，自定义截断分布用逆变换采样的实现成本远低于适配通用分布接口。

内容的提问来源于stack exchange，提问作者Aleksandar Knežević