从任意概率密度函数生成随机数答案

【问题标题】：Generating random numbers from arbitrary probability density function从任意概率密度函数生成随机数
【发布时间】：2017-05-30 04:31:38
【问题描述】：

我希望能够生成具有来自绘制曲线的概率密度函数的随机数。下面这两个具有相同的曲线下面积，但应生成具有不同特征的随机数列表。

我的直觉是这样做的一种方法是对曲线进行采样，然后使用这些矩形的区域来提供np.random.choice 以选择一个范围以在该矩形范围的范围内进行普通随机.

感觉这不是一种非常有效的方法。有没有更“正确”的方法呢？

我真的很擅长做这件事：

import matplotlib.pyplot as plt
import numpy as np

areas = [4.397498, 4.417111, 4.538467, 4.735034, 4.990129, 5.292455, 5.633938,
         6.008574, 6.41175, 5.888393, 2.861898, 2.347887, 2.459234, 2.494357,
         2.502986, 2.511614, 2.520243, 2.528872, 2.537501, 2.546129, 7.223747,
         7.223747, 2.448148, 1.978746, 1.750221, 1.659351, 1.669999]
divisons = [0.0, 0.037037, 0.074074, 0.111111, 0.148148, 0.185185, 0.222222,
            0.259259, 0.296296, 0.333333, 0.37037, 0.407407, 0.444444, 0.481481,
            0.518519, 0.555556, 0.592593, 0.62963, 0.666667, 0.703704, 0.740741,
            0.777778, 0.814815, 0.851852, 0.888889, 0.925926, 0.962963, 1.0]
weights = [a/sum(areas) for a in areas]
indexes = np.random.choice(range(len(areas)), 50000, p=weights)
samples = []
for i in indexes:
    samples.append(np.random.uniform(divisons[i], divisons[i+1]))

binwidth = 0.02
binSize = np.arange(min(samples), max(samples) + binwidth, binwidth)
plt.hist(samples, bins=binSize)
plt.xlim(xmax=1)
plt.show()

方法貌似可行，就是有点重！

【问题讨论】：

您是说您只有一个带有该曲线的图像文件吗？或者你真的有代表曲线上点坐标的数字吗？
也可以。它可能是图像文件，但更可能是绘制的曲线。 svg 或触摸屏上的某种墨水。
SVG 是一个图像文件。如果它是在屏幕上绘制的，那么您的程序如何访问它？我问的是您的程序将使用的数据 format 是什么，而不是事物的创建方式。'
目前只是假设性的。我在 CAD 程序中对其进行原型设计，但它可能会出现在任何地方。我假设您的意思是位图，可以访问 SVG 曲线中的坐标。（最终！）
从数学上讲，我要做的是整合 PDF 以获得累积分布函数。如果然后将其反转，您将获得一个函数，您可以在其中插入一个随机数 [0, 1] 并有效地从原始分布中获取一个值。您实际如何做到这一点取决于您的数据格式。

标签： python random statistics

【解决方案1】：

对于您的情况，似乎基于直方图的方法肯定是最简单的，因为您有一条用户已绘制的线。

但由于您只是试图从该分布中生成随机数，您可以在下面的函数中直接使用归一化的 y 值（将所有像素的 y 位置相加并除以总数）作为概率分布只取用户绘制的像素数大小的数组。

from numpy.random import choice
pde = choice(list_of_candidates, number_of_items_to_pick, p=probability_distribution)

probability_distribution（归一化像素 y 值）是与 list_of_candidates（关联的 x 值）顺序相同的序列。您还可以使用关键字 replace=False 来更改行为，以便绘制的项目不会被替换。

see numpy docs here

这应该快得多，因为您实际上并没有生成整个 pde，只是绘制与 pde 匹配的随机数。

编辑：您的更新看起来是一种可靠的方法。如果您确实想生成 pde，您可以考虑研究 numba (http://numba.pydata.org) 以矢量化您的 for 循环。

【讨论】：

【解决方案2】：

一种方法是使用 scipy.stats 中的 rv_continuous。最直接的方法是用一组带有 rv_continuous 的样条线来近似其中一个 pdf。事实上，你可以通过用这个东西定义一个 pdf 或一个 cdf 来生成伪随机偏差。

【讨论】：

【解决方案3】：

另一种方法是对 CDF 的逆进行采样。然后使用统一随机生成器在逆 CDF 的 x 轴上生成 p 值，以生成 PDF 的随机绘制。见这篇文章：http://matlabtricks.com/post-44/generate-random-numbers-with-a-given-distribution

【讨论】：

您的链接现在已损坏。我在答案中添加了一个明确的自包含代码，以使您的提议更加具体

【解决方案4】：

我在使用 rv_continuous 时遇到了问题，所以我做了自己的小例程，从任何具有紧凑支持的连续分布中采样，例如来自两个指数的总和，或来自任何已知的离散 pdf（如问题中所问）。这本质上是@Jan 的解决方案（一个非常经典的解决方案）。

我的代码是完全独立的。要使其适应任何其他分布，您只需更改 unnormalized_pdf 中的公式，并确保正确设置支持的边界（在我的情况下，从 0 到 10/lambda_max 就足够了。

import numpy as np
import matplotlib.pyplot as plt

plt.ion()

## The function may be any function, so long as it is with FINITE Support
def unnormalized_pdf(T, lambda1, intercept1, lambda2, intercept2):
    return np.exp(-lambda1 * T - intercept1) + np.exp(-lambda2 * T - intercept2)


lambda1, intercept1, lambda2, intercept2 = (
    0.0012941708402716523,
    8.435217547457713,
    0.0063804460354380385,
    6.712937938322769,
)

## defining the support of the pdf by hand
x0 = 0
xmax = max(1 / lambda1, 1 / lambda2) * 10

## the more bins, the higher the precision
Nbins = 1000000
xs = np.linspace(x0, xmax, Nbins)
dx = xs[1] - xs[0]
## other way to specify it:
# dx = min(1/lambda1, 1/lambda2)/100
# xs = np.arange(x0, xmax, dx)

## compute the (approximate) pdf and cdf of the thing to sample:
pdf = unnormalized_pdf(xs, lambda1, intercept1, lambda2, intercept2)
normalized_pdf = pdf / pdf.sum()
cdf = np.cumsum(normalized_pdf)

## sampling from the distro
Nsamples = 100000
r = np.random.random(Nsamples)
indices_in_cdf = np.searchsorted(cdf, r)
values_drawn = xs[indices_in_cdf]
histo, bins = np.histogram(values_drawn, 1000, density=True)
plt.semilogy(bins[:-1], histo, label="drawn from distro", color="blue")
plt.semilogy(xs, normalized_pdf / dx, label="exact pdf from which we sample", color="k", lw=3)
plt.legend()
plt.show()

【讨论】：