如何从Python中给定CDF的分布中采样

Question

我想从概率分布中抽取样本 CDF

1 - e^(-x^2)

。

python/scipy/等中有没有方法？使您能够从仅给定 CDF 的概率分布中进行采样？

Answer 1

要创建给定 CDF 的自定义随机变量类，您可以子类

scipy.rv_continuous

并覆盖

rv_continuous._cdf

。然后，这将自动生成相应的 PDF 和有关您的发行版的其他统计信息，例如

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

class MyRandomVariableClass(stats.rv_continuous):
    def __init__(self, xtol=1e-14, seed=None):
        super().__init__(a=0, xtol=xtol, seed=seed)

    def _cdf(self, x):
        return 1-np.exp(-x**2)


if __name__ == "__main__":
    my_rv = MyRandomVariableClass()

    # sample distribution
    samples = my_rv.rvs(size = 1000)

    # plot histogram of samples
    fig, ax1 = plt.subplots()
    ax1.hist(list(samples), bins=50)

    # plot PDF and CDF of distribution
    pts = np.linspace(0, 5)
    ax2 = ax1.twinx()
    ax2.set_ylim(0,1.1)
    ax2.plot(pts, my_rv.pdf(pts), color='red')
    ax2.plot(pts, my_rv.cdf(pts), color='orange')

    fig.tight_layout()
    plt.show()

Answer 2

逆变换采样

要添加到Heike的解决方案中，您可以使用逆变换采样通过CDF进行采样：

import math, random
import matplotlib.pyplot as plt

def inverse_cdf(y):
    # Computed analytically
    return math.sqrt(math.log(-1/(y - 1)))

def sample_distribution():
    uniform_random_sample = random.random()
    return inverse_cdf(uniform_random_sample)

x = [sample_distribution() for i in range(10000)]
plt.hist(x, bins=50)
plt.show()

SciPy 是如何做到的

我也很好奇这在 SciPy 中是如何工作的。实际上看起来它所做的事情与上面的非常相似。基于 SciPy 文档：

默认方法 _rvs 依赖于应用于均匀随机变量的 cdf 的逆 _ppf。为了有效地生成随机变量，要么需要覆盖默认的 _ppf（例如，如果逆 cdf 可以以显式形式表示），要么需要在自定义 _rvs 方法中实现采样方法。

并且基于 SciPy 源代码，如果未明确指定，

_ppf

（即 CDF 的逆）实际上看起来是数字近似值。非常酷！

Answer 3

我的答案建立在@Heike 和@Jeff N 的答案的基础上（这都是正确的，不要误会我的意思）。但是，我必须混合使用定义 both

rv_continuous._ppf

、and

rv_continuous._pdf

来处理更大的样本值，而且这个解决方案也相当快。

就我而言，我必须模拟超过 50,000 个值的随机变量。当仅定义

rv_continuous._pdf

时，输出需要很长时间，并给出 ValueError，遇到 NaN 值，并停止运行器。

当仅定义

rv_continuous._ppf

时，函数的折线图绘制时出现问题，因为它超出了 Python 中的递归限制（尽管可以设置更高的递归限制，但它并没有为我解决该问题）。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

class dist_gen(stats.rv_continuous):
  def __init__(self, seed=233423):
        super().__init__(seed=seed)
  
  # Both _pdf and _ppf overriden in scipy.stats.rv_continuous
  def _pdf(self, x):
    return 0.5*np.exp(-1*(abs(x)))
  
  def _ppf(self, x):
    return -1*np.sign(x-0.5)*np.log(1-(2*abs(x-0.5)))

if __name__ == "__main__":
  # Number of variables to be simulated
  N = 100000
  dist = dist_gen()

  samples = dist.rvs(size=N)

  # Plot histogram
  fig, ax = plt.subplots()
  ax.hist(samples, bins=200, density=True)

  # Plotting the line function
  domain = np.linspace(-7.5, 7.5)
  ax.set_ylim(0, 0.75)
  ax.plot(domain, dist.pdf(domain), color='purple')

  # Showing the final plot
  plt.show()