根据分布函数生成随机数

Question

如果我们想要均匀地生成区间 [a,b] 内的随机数，我们可以很容易地生成它：

A=rand()*(b-a) + a

其中

rand()

是一个可以生成 0 到 1 之间均匀随机数的函数。所以

A

是

[a,b]

中的随机数。

为了根据

y=x-x^2

这样的分布函数生成随机数，我遇到了一个问题。

我想使用here提到的方法。但我对使用 python 函数不感兴趣

inverse_cdf(np.random.uniform())

。

我可以通过 0 和 X 的积分来计算函数“y”的 CDF，我称之为“f”。但是当我将 rand() 函数（0 到 1 之间的数字）代入 f 的反函数时，我得到了一个复数！ 这意味着：

A=f^(-1) (rand())

返回一个复数。

基于分布函数生成随机数的正确方法是吗？

我使用这个网站来计算

f=x^2/2 - x^3/3

的逆，下面的代码是计算的一部分，显示

tmp1

始终为负

for i=1:10
rnd1=rand;
tmp1  = 2*sqrt(6)*sqrt(6*rnd1^2-rnd1)-12*rnd1+1
cTmp1 = tmp1^(1/3)
end

Answer 1

问题是：“基于分布函数生成随机数的正确方法是吗？”因此，基本上您有连续概率密度函数 (PDF) 或带有符号 f(x) 的离散概率质量函数 (PMF)，并且您正在寻找一种找到随机变量 x 的方法。

至少有两种方法可以做到这一点。

1. 使用逆变换分布。
2. 使用拒绝方法

使用逆变换： 如果我们知道概率分布函数，那么对于某些累积分布函数（CDF），我们可以找到随机变量的闭集。假设你的概率函数是f(x)，CDF是F(x)，那么假设你能得到反函数，你就能得到随机变量

x=inverse F(U)

其中 U 是随机均匀分布

使用拒绝方法： 如果 CDF 没有封闭形式逆，那么您始终可以使用拒绝方法。拒绝法的思想是生成一个随机的二维点（一对随机数）：（r1，r2），然后该点要么在PDF的曲线下方，要么在PDF的曲线上方。如果该点位于曲线下方，则我们将该值作为随机数，否则我们对另一个点进行采样。 假设 PDF f(x) 以最大值 M

为界

r1 is generated within interval [a, b] of horizontal axis
r2 is generated within interval [0, M] of vertical axis
If r2 <f (r1) then accept r1 and output r1
   Else reject r1 and generate new random point

如果你能找到CDF的逆，逆变换方法就优于拒绝方法，因为你可以得到封闭形式。例如，比率为 1 的指数分布的 CDF 为 F(x) = 1-exp(-x)。逆变换将是

x= inverse F(U) = -ln(1-U) = -ln(U)

因为 1-U2=U2

Answer 2

在过去的几年里，SciPy 中添加了一些不错的新工具来解决 Python 中的此类问题。 只需提供一些有关分布的信息（例如密度/pdf），您就可以轻松地从自定义连续或离散单变量分布生成样本。

不同方法的概述： https://docs.scipy.org/doc/scipy/reference/stats.sampling.html

教程： https://docs.scipy.org/doc/scipy/tutorial/stats/sampling.html

这是自定义密度的示例：

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