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

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()

这个方法看似有效，但是有点重！

Answer 1

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

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

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 匹配的随机数。

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

Answer 2

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

Answer 3

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

Answer 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()

Answer 5

您对使用与密度成比例的权重进行采样的直觉作为近似值是很好的。我建议在积分密度后使用工具来反转分布函数。这是一个例子：

import numpy as np
from scipy.stats.sampling import NumericalInversePolynomial
from matplotlib import pyplot as plt
from scipy.integrate import quad


class MyDist:
    def __init__(self, a):
        self.a = a

    def support(self):
        # distribution restricted to 0, 5, can be changed as needed
        return (0, 5)

    def pdf(self, x):
        # this is not a proper pdf, the normalizing
        # constant is missing (does not integrate to one)
        # this is ok for the method presented here
        return x * (x + np.sin(5*x) + 2) * np.exp(-x**self.a)


dist = MyDist(0.5)
gen = NumericalInversePolynomial(dist)

# compute the missing normalizing constant to plot the pdf
const_pdf = quad(dist.pdf, *dist.support())[0]

r = gen.rvs(size=50000)
x = np.linspace(r.min(), r.max(), 500)

# show histogram together with the pdf
plt.plot(x, dist.pdf(x) / const_pdf)
plt.hist(r, density=True, bins=100)
plt.show()

还有更多工具可以通过提供一些有关分布的信息（例如密度/pdf）来从自定义连续或离散单变量分布中进行采样。不同方法的概述： https://docs.scipy.org/doc/scipy/reference/stats.sampling.html