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

Question

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

我的直觉是，一种方法是对曲线进行采样，然后使用这些矩形的面积来提供一个 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()

这种方法似乎可行，但有些繁重！

- Ben

---

你是说你只有一个带有那条曲线的图像文件吗？还是你实际上有代表曲线上点坐标的数字？ - BrenBarn

可能是两者之一。它可能是一个图像文件，但更有可能是手绘的曲线。可以是SVG格式或在触摸屏上进行的某种墨水绘制。 - Ben

SVG是一种图像文件。如果它在屏幕上绘制，那么你的程序如何访问它？我想知道你的程序将使用什么数据格式，而不是这个东西是如何创建的。 - BrenBarn

目前还只是假设性的。我正在 CAD 程序中进行原型设计，但最终可能会到任何地方。我假设你指的是位图，在 SVG 曲线中可以访问坐标。（最终！） - Ben

从数学角度来看，我会对概率密度函数进行积分，以得到累积分布函数。然后，如果你将其反转，就可以得到一个函数，可以将[0,1]范围内的随机数插入其中，并有效地从原始分布中获取一个值。如何实际操作取决于数据的格式。 - BrenBarn

Answer 1

一种方法是使用来自scipy.stats的rv_continuous。开始的简单方法是使用rv_continuous用一组样条函数来近似其中一个概率密度函数。事实上，您可以使用此工具通过定义pdf或cdf之一来生成伪随机变量。

- Bill Bell

Answer 2

对于你的情况，似乎基于直方图的方法肯定是最容易的，因为你有一个用户绘制的线。

但是，如果你只是想从该分布中生成随机数，可以直接在下面的函数中使用标准化的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来更改行为，以便不替换已抽取的项目。

在此处查看numpy文档

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

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

- Greg Jennings

Answer 3

另一种方法是对CDF的倒数进行抽样。然后，您可以使用均匀分布的随机生成器在倒数CDF的x轴上生成p值，以生成PDF的随机抽样。

请参阅此文章：http://matlabtricks.com/post-44/generate-random-numbers-with-a-given-distribution

- Jan

---

你提供的链接已经失效了。为了让你的建议更加具体明确，我在我的回答中添加了一个自包含代码。 - François Landes

Answer 4

我在使用rv_continuous时遇到了问题，所以我自己编写了一个小程序来从任何具有紧支撑的连续分布中进行抽样，例如从两个指数和或任何已知离散概率密度函数中进行抽样（如问题所述）。 这本质上是@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()

- François Landes