如何从（任意）连续概率分布进行模拟？ [复制]答案

【问题标题】：How to simulate from an (arbitrary) continuous probability distribution? [duplicate]如何从（任意）连续概率分布进行模拟？ [复制]
【发布时间】：2015-05-19 15:27:45
【问题描述】：

我有一个这样的概率密度函数：

def p1(x):
    return ( sin(x) ** (-0.75) ) / (4.32141 * (x ** (1/5)))

我想用这个pdf 拒绝[0; 1] 上的随机值。如何做随机值？

【问题讨论】：

从任意分布采样可以通过在 [0, 1] 中均匀采样然后使用累积密度函数的倒数来完成。如果您无法解析地计算此逆 cdf，则可以使用 pdf 的数值积分（例如使用梯形规则）并将值存储在列表中；然后你可以对你的 [0, 1] 样本进行二分搜索以找到你的 pdf 样本。
@FrancisColas，这是解决我的任务的典型方法。但是，当有 scipy 这么强大的工具时，我真的可以用经典的方式解决我的问题吗？
好吧，我不知道使用 numpy/scipy 的任何具体方式。但请注意，这只是为您提供解决方案提示的评论，其他人可能会发布更令人满意的答案。
@FrancisColas，无论如何，谢谢你。
@FrancisColas 最好将其与quad 集成，然后插入反函数。

标签： python python-3.x random scipy distribution

【解决方案1】：

不使用 scipy 并给定 PDF 的数字采样，您可以使用累积分布和线性插值进行采样。下面的代码假定 x 中的间距相等。可以对其进行修改以对任意采样的 PDF 进行集成。请注意，它将 PDF 重新规范化为 x 范围内的 1。

import numpy as np

def randdist(x, pdf, nvals):
    """Produce nvals random samples from pdf(x), assuming constant spacing in x."""

    # get cumulative distribution from 0 to 1
    cumpdf = np.cumsum(pdf)
    cumpdf *= 1/cumpdf[-1]

    # input random values
    randv = np.random.uniform(size=nvals)

    # find where random values would go
    idx1 = np.searchsorted(cumpdf, randv)
    # get previous value, avoiding division by zero below
    idx0 = np.where(idx1==0, 0, idx1-1)
    idx1[idx0==0] = 1

    # do linear interpolation in x
    frac1 = (randv - cumpdf[idx0]) / (cumpdf[idx1] - cumpdf[idx0])
    randdist = x[idx0]*(1-frac1) + x[idx1]*frac1

    return randdist

【讨论】：

好答案 - 我将它用于一个非常复杂的 pdf（两个法线的比率）并且它很有魅力（而且它非常快） - 你应该在链接的 SO 帖子中发布你的答案，因为这个被标记为重复 - 我认为你的答案更容易实现 - 这是 pdf：static.stevereads.com/papers_to_read/…

【解决方案2】：

正如 Francis 所说，您最好了解您的发行版的 cdf。无论如何，scipy 提供了一种方便的方式来定义自定义分布。看起来很像

from scipy import stats
class your_distribution(stats.rv_continuous):
    def _pdf(self, x):
        return ( sin(x) ** (-0.75) ) / (4.32141 * (x ** (1/5)))

distribution = your_distribution()
distribution.rvs()

【讨论】：

(+1)，很好的答案。如果结果太慢（因为它通过通用代码路径），那么只有集成/内插 cdf 是值得的，并且它对于随机采样是相反的。
@Gioelelm，谢谢。
pdf_z 是 scipy.interpolate.interp1d 的插值对象，如果超出 x 范围，则 bounds=0 ！ from scipy import stats class TransfromPriorPDF(stats.rv_continuous): def _pdf(self, x): return pdf_z(x) dist = TransfromPriorPDF() dist.rvs() 这给了我一个运行时错误：--> 776 r = _zeros。 _brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp) 777 return results_c(full_output, r) 778 RuntimeError: 100 次迭代后无法收敛。知道哪里出了问题吗？