以 Notebook 形式打开

随机变量转换指南

背景

在 SciPy 1.15 之前，SciPy 的所有连续概率分布（例如 scipy.stats.norm）都是 scipy.stats.rv_continuous 子类的实例。

from scipy import stats
dist = stats.norm
type(dist)

scipy.stats._continuous_distns.norm_gen

isinstance(dist, stats.rv_continuous)

True

有两种明显的方法来操作这些对象。

根据更常见的方式，“参数”（例如 x）和“分布参数”（例如 loc、scale）都作为对象方法的输入提供。

x, loc, scale = 1., 0., 1.
dist.pdf(x, loc, scale)

0.24197072451914337

不太常见的方法是调用分布对象的 __call__ 方法，它返回 rv_continuous_frozen 的实例，而与原始类无关。

frozen = stats.norm()
type(frozen)

scipy.stats._distn_infrastructure.rv_continuous_frozen

这个新对象的方法只接受参数，而不接受分布参数。

frozen.pdf(x)

0.24197072451914337

从某种意义上说，像 norm 这样的 rv_continuous 实例代表“分布族”，需要参数来标识特定的概率分布，而 rv_continuous_frozen 的实例类似于“随机变量”——一个遵循特定概率分布的数学对象。

这两种方法都是有效的，并且在某些情况下具有优势。例如，对于这种简单的调用或者当方法被视为分布参数的函数而不是通常的参数（例如，似然）时，stats.norm.pdf(x) 似乎比 stats.norm().pdf(x) 更自然。但是，前一种方法有一些固有的缺点；例如，SciPy 的所有 125 个连续分布都必须在导入时实例化，每次调用方法时都必须验证分布参数，并且方法的文档必须 a) 为每个分布的每个方法单独生成，或者 b) 省略每个分布特有的形状参数。为了解决这些和其他缺点，gh-15928 提出了一种新的、独立的基础设施，该基础设施基于后一种（随机变量）方法。本转换指南记录了 rv_continuous 和 rv_continuous_frozen 的用户如何迁移到新的基础设施。

注意：新的基础设施对于某些用例可能还不太方便，尤其是在将分布参数拟合到数据时。如果旧的基础设施适合他们的需求，欢迎用户继续使用；目前它没有被弃用。

基础

在新的基础设施中，分布族是按照 CamelCase 约定命名的类。它们必须在使用前实例化，参数作为仅关键字参数传递。分布族类的实例可以被认为是随机变量，在数学中通常使用大写字母表示。

from scipy import stats
mu, sigma = 0, 1
X = stats.Normal(mu=mu, sigma=sigma)
X

Normal(mu=0.0, sigma=1.0)

实例化后，形状参数可以作为属性读取（但不能写入）。

X.mu, X.sigma

(0.0, 1.0)

scipy.stats.Normal 的文档包含指向其每个方法详细文档的链接。（例如，与 scipy.stats.norm 的文档进行比较。）

注意：尽管截至 SciPy 1.15，只有 Normal 和 Uniform 具有呈现的文档，并且可以直接从 stats 导入，但几乎所有其他旧分布都可以通过 make_distribution（稍后讨论）与新接口一起使用。

然后通过仅传递参数（如果有），而不传递分布参数来调用方法。

X.mean()

0.0

X.pdf(x)

0.24197072451914337

对于像这样的简单调用（例如，参数是有效的浮点数），对新随机变量的方法的调用通常比对旧分布方法的类似调用更快。

%timeit X.pdf(x)

2.47 µs ± 5.18 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

%timeit frozen.pdf(x)

21.8 µs ± 161 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

%timeit dist.pdf(x, loc=mu, scale=sigma)

22 µs ± 21.2 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

请注意，X.pdf 和 frozen.pdf 的调用方式是相同的，而对 dist.pdf 的调用非常相似 - 唯一的区别在于对 dist.pdf 的调用包括形状参数，而在新的基础设施中，形状参数仅在实例化随机变量时提供。

除了 pdf 和 mean 之外，新基础设施的几个其他方法本质上与旧方法相同。

• logpdf（概率密度函数的对数）

• cdf（累积分布函数）

• logcdf（累积分布函数的对数）

• entropy（微分熵）

• 中位数

• 支持

其他方法有新的名称，但在其他方面是直接替换。

• sf（生存函数） \(\rightarrow\) ccdf（互补累积分布函数）

• logsf \(\rightarrow\) logccdf

• ppf（百分点函数） \(\rightarrow\) icdf（逆累积分布函数）

• isf（逆生存函数） \(\rightarrow\) iccdf（逆互补累积分布函数）

• std \(\rightarrow\) standard_deviation

• var \(\rightarrow\) variance

新的基础设施有几个与上述相同类型的新方法。

• ilogcdf（累积分布函数对数的逆）

• ilogccdf（互补累积分布函数对数的逆）

• logentropy（熵的对数）

• mode（分布的众数）

• 偏度

• kurtosis（非超额峰度；请参阅下面的“标准化矩”）

为了方便起见，它有一个新的 plot 方法

import matplotlib.pyplot as plt
X.plot()
plt.show()

旧基础设施的大部分剩余方法（rvs, moment, stats, interval, fit, nnlf, fit_loc_scale 和 expect）具有类似的功能，但需要注意一些事项。在描述替代方案之前，我们简要介绍如何处理非正态分布的随机变量：几乎所有旧的分布对象都可以通过 scipy.stats.make_distribution 转换为新的分布类，并且可以通过将形状参数作为关键字参数传递来实例化新的分布类。例如，考虑 Weibull 分布。我们可以创建一个新的类，它是分布族的抽象，如下所示

Weibull = stats.make_distribution(stats.weibull_min)
print(Weibull.__doc__[:288])  # help(Weibull) prints this and (too much) more

This class represents `scipy.stats.weibull_min` as a subclass of `ContinuousDistribution`.
The `repr`/`str` of class instances is `Weibull`.
The PDF of the distribution is defined for :math:`x` in [0.0, \infty].
This class accepts one parameterization:
`c` for :math:`c \in (0, \infty)`.

根据 weibull_min 的文档，形状参数表示为 c，因此我们可以通过将 c 作为关键字参数传递给新的 Weibull 类来实例化一个随机变量。

c = 2.
X = Weibull(c=c)
X.plot()
plt.show()

以前，所有分布都继承了 loc 和 scale 分布参数。这些不再被接受。相反，随机变量可以使用算术运算符进行平移和缩放。

Y = 2*X + 1
Y.plot()  # note the change in the abscissae
plt.show()

提供了一个单独的分布 weibull_max，它是关于原点的 weibull_min 的反射。现在，这只是 -X。

Y = -X
Y.plot()
plt.show()

矩

旧的基础设施提供了一个 moments 方法来计算原始矩。当仅指定阶数时，新的 moment 方法是直接替换。

stats.weibull_min.moment(1, c), X.moment(1)  # first raw moment of the Weibull distribution with shape c

(0.8862269254527579, 0.8862269254527579)

然而，旧的基础设施也有一个 stats 方法，它提供了分布的各种统计信息。以下统计信息与指示的字符相关联。

• 均值（关于原点的第一个原始矩，'m'）

• 方差（第二个中心矩，'v'）

• 偏度（第三个标准化矩，'s'）

• 超额峰度（第四个标准化矩减去 3，'k'）

例如

stats.weibull_min.stats(c, moments='mvsk')

(0.8862269254527579,
 0.21460183660255183,
 0.6311106578189344,
 0.24508930068764556)

现在，可以通过将适当的参数传递给新的 moment 方法来计算任何 order 和 kind（原始、中心和标准化）的矩。

(X.moment(order=1, kind='raw'), 
 X.moment(order=2, kind='central'), 
 X.moment(order=3, kind='standardized'), 
 X.moment(order=4, kind='standardized'))

(0.8862269254527579,
 0.21460183660255183,
 0.6311106578189344,
 3.2450893006876314)

请注意 峰度 定义的差异。以前，提供的是“超额”（又名“Fisher”）峰度。按照惯例（而不是自然的数学定义），这是标准化四阶矩移动一个常数值（3），使正态分布的值为 0.0。

stats.norm.stats(moments='k')

0.0

新的 moment 和 kurtosis 方法默认不遵循此惯例；它们根据数学定义报告标准化四阶矩。这也称为“非超额”（或“Pearson”）峰度，对于标准正态分布，它的值为 3.0。

stats.Normal().moment(4, kind='standardized') 

3.0

为了方便起见，kurtosis 方法提供了一个 convention 参数来选择两者之一。

stats.Normal().kurtosis(convention='excess'), stats.Normal().kurtosis(convention='non-excess')

(0.0, 3.0)

随机变量

旧基础设施的 rvs 方法已被 sample 替换，但是当 1) 需要控制随机状态或 2) 将数组用于形状、位置或尺度参数时，应注意两个差异。

首先，参数 random_state 根据 SPEC 7 被 rng 替换。过去控制随机状态的一种模式是使用 numpy.random.seed 或将整数种子传递给 random_state 参数。整数种子被转换为 numpy.random.RandomState 的实例，因此在两种情况下，给定整数种子的行为将是相同的

import numpy as np
np.random.seed(1)
dist = stats.norm
dist.rvs(), dist.rvs(random_state=1)

(1.6243453636632417, 1.6243453636632417)

在新基础设施中，将 numpy.Generator 的实例传递给 sample 的 rng 参数以控制随机状态。请注意，整数种子现在被转换为 numpy.random.Generator 的实例，而不是 numpy.random.RandomState 的实例。

X = stats.Normal()
rng = np.random.default_rng(1)  # instantiate a numpy.random.Generator
X.sample(rng=rng), X.sample(rng=1)

(0.345584192064786, 0.345584192064786)

其次，参数 shape 替换参数 size。

dist.rvs(size=(3, 4))

array([[-0.61175641, -0.52817175, -1.07296862,  0.86540763],
       [-2.3015387 ,  1.74481176, -0.7612069 ,  0.3190391 ],
       [-0.24937038,  1.46210794, -2.06014071, -0.3224172 ]])

X.sample(shape=(3, 4))

array([[ 1.09030882,  1.00534614,  2.31221283, -1.19646394],
       [ 0.68798893,  0.75769842,  1.52220417,  0.3066763 ],
       [ 1.15948958,  1.51430283, -1.01990881, -0.14666494]])

除了名称上的差异之外，当涉及数组形状参数时，行为上也存在差异。以前，分布参数数组的形状必须包含在 size 中。

dist.rvs(size=(3, 4, 2), loc=[0, 1]).shape  # `loc` has shape (2,)

(3, 4, 2)

现在，参数数组的形状被认为是随机变量对象本身的属性。指定数组形状参数的形状是多余的，因此在指定样本的 shape 时不包括它。

Y = stats.Normal(mu = [0, 1])
Y.sample(shape=(3, 4)).shape  # the sample has shape (3, 4); each element is of shape (2,)

(3, 4, 2)

概率质量（又名“置信”）区间

旧的基础设施提供了一个 interval 方法，默认情况下，该方法将返回一个包含分布指定百分比概率质量的对称（关于中位数）区间。

a = 0.95
dist.interval(confidence=a)

(-1.959963984540054, 1.959963984540054)

现在，调用具有每个尾部所需概率质量的逆 CDF 和互补 CDF 方法。

p = 1 - a
X.icdf(p/2), X.iccdf(p/2)

(-1.959963984540054, 1.959963984540054)

似然函数

旧的基础设施提供了一个函数来计算负对数似然函数，错误地命名为 nnlf（而不是 nllf）。它接受分布的参数作为元组，观测值作为数组。

mu = 0
sigma = 1
data = stats.norm.rvs(size=100, loc=mu, scale=sigma)
stats.norm.nnlf((mu, sigma), data)

131.07116392314364

现在，只需根据其数学定义计算负对数似然。

X = stats.Normal(mu=mu, sigma=sigma)
-X.logpdf(data).sum()

131.07116392314364

期望值

旧基础设施的 expect 方法估计由分布的 PDF 加权的任意函数的定积分。例如，分布的四阶矩由下式给出

def f(x): return x**4
stats.norm.expect(f, lb=-np.inf, ub=np.inf)

3.000000000000001

这与源代码所做的相比，并没有增加多少便利性：使用 scipy.integrate.quad 执行数值积分。

from scipy import integrate
def f(x): return x**4 * stats.norm.pdf(x)
integrate.quad(f, a=-np.inf, b=np.inf)  # integral estimate, estimate of the error

(3.0000000000000053, 1.2043244026618655e-08)

当然，这也可以使用新的基础设施轻松完成。

def f(x): return x**4 * X.pdf(x)
integrate.quad(f, a=-np.inf, b=np.inf)  # integral estimate, estimate of the error

(3.0000000000000053, 1.2043244043338238e-08)

conditional 参数只是将结果按区间内包含的概率质量的倒数进行缩放。

a, b = -1, 3
def f(x): return x**4
stats.norm.expect(f, lb=a, ub=b, conditional=True)

1.657813862143119

直接使用 quad，即

def f(x): return x**4 * stats.norm.pdf(x)
prob = stats.norm.cdf(b) - stats.norm.cdf(a)
integrate.quad(f, a=a, b=b)[0] / prob

1.6578138621431193

以及新的随机变量

integrate.quad(f, a=a, b=b)[0] / X.cdf(a, b)

1.657813862143119

请注意，新随机变量的 cdf 方法接受两个参数来计算指定区间内的概率质量。在许多情况下，当概率质量非常小时，这可以避免减法误差（“灾难性抵消”）。

拟合

位置/尺度估计

旧的基础设施提供了一个从数据估计分布的位置和尺度参数的函数。

rng = np.random.default_rng(91392794601852341588152500276639671056)
dist = stats.weibull_min
c, loc, scale = 0.5, 0, 4.
data = dist.rvs(size=100, c=c, loc=loc, scale=scale, random_state=rng)
dist.fit_loc_scale(data, c)
# compare against 0 (default location) and 4 (specified scale)

(1.8268432669973347, 1.9640905288934327)

根据源代码，很容易复制该方法。

X = Weibull(c=c)
def fit_loc_scale(X, data):
    m, v = X.mean(), X.variance()
    m_, v_ = data.mean(), data.var()
    scale = np.sqrt(v_ / v)
    loc = m_ - scale*m
    return loc, scale

fit_loc_scale(X, data)

(1.8268432669973347, 1.9640905288934327)

请注意，在此示例中，估计值非常差，并且启发式算法的性能不佳是决定不提供替代方案的原因。

最大似然估计

旧基础设施的最后一个方法是 fit，它根据分布的样本估计基础分布的位置、尺度和形状参数。

c_, loc_, scale_ = stats.weibull_min.fit(data)
c_, loc_, scale_

(0.5960699537061247, 1.8950104942593064e-05, 4.015875612107546)

此方法的便利性既是祝福也是诅咒。

当它起作用时，这种简单性非常受欢迎。对于某些分布，已经手动编程了参数的最大似然估计 (MLE) 的解析表达式，对于这些分布，fit 方法非常可靠。

然而，对于大多数分布，拟合是通过数值最小化负对数似然函数来完成的。这并不能保证成功——既有MLE的固有局限性，也有数值优化技术的现状的限制——在这些情况下，fit 的用户会陷入困境。但是，只要稍加理解，就能大大提高成功的可能性，因此，我们在这里介绍一个使用新基础设施进行拟合的迷你教程。

首先，请注意，MLE 是将分布拟合到数据的几种潜在策略之一。它所做的只是找到使负对数似然函数 (NLLF) 最小的分布参数值。请注意，对于生成数据的分布，NLLF 为

stats.weibull_min.nnlf((c, loc, scale), data)

246.9503871545133

使用 fit 估计的参数的 NLLF 更低

stats.weibull_min.nnlf((c_, loc_, scale_), data)

230.35385152369787

因此，无论这是否满足用户对“良好拟合”的概念，或者参数是否在合理范围内，fit 都认为它的工作已经完成。除了允许用户指定分布参数的猜测外，它几乎没有对过程的控制。（*注意：从技术上讲，可以使用 optimizer 参数进行更强的控制，但这没有很好的文档记录，并且它的使用会抵消 fit 方法提供的任何便利性。*）

另一个问题是，对于某些分布，MLE 本质上表现不佳。例如，在这个例子中，我们可以通过将位置设置为数据的最小值来使 NLLF 尽可能低——即使对于形状和尺度参数的虚假值也是如此。

stats.weibull_min.nnlf((0.1, np.min(data), 10), data)

-inf

与这些问题相伴的是，fit 默认情况下估计所有分布参数，包括位置。在上面的例子中，我们可能只对更常见的双参数 Weibull 分布感兴趣，因为实际位置是零。fit *可以* 通过指定位置为固定参数来适应这种情况。

# c_, loc_, scale_ = stats.weibull_min.fit(data, loc=0)  # careful! this provides loc=0 as a *guess*
c_, loc_, scale_ = stats.weibull_min.fit(data, floc=0)
c_, loc_, scale_

(0.5834128280886433, 0, 3.829972848420729)

stats.weibull_min.nnlf((c_, loc_, scale_), data)

244.9617429249505

但是，通过提供一个据称“只需工作”的 fit 方法，会鼓励用户忽略一些重要的细节。

相反，我们建议，使用 SciPy 提供的通用优化工具将分布拟合到数据几乎同样简单。但是，这允许用户根据其首选的目标函数执行拟合，在出现意外情况时提供选项，并帮助用户避免最常见的陷阱。

例如，假设我们期望的目标函数确实是负对数似然。将其定义为仅形状参数 c 和尺度 的函数很简单，将位置保留为默认值零。

def nllf(x):
    c, scale = x  # pass (only) `c` and `scale` as elements of `x`
    X = Weibull(c=c) * scale
    return -X.logpdf(data).sum()

nllf((c_, scale_))

244.96174292495053

为了执行最小化，我们可以使用 scipy.optimize.minimize。

from scipy import optimize
eps = 1e-10  # numerical tolerance to avoid invalid parameters
lb = [eps, eps]  # lower bound on `c` and `scale`
ub = [10, 10]  # upper bound on `c` and `scale`
x0 = [1, 1]  # guess to get optimization started
bounds = optimize.Bounds(lb, ub)
res_mle = optimize.minimize(nllf, x0, bounds=bounds)
res_mle

  message: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
  success: True
   status: 0
      fun: 244.96174292277675
        x: [ 5.834e-01  3.830e+00]
      nit: 13
      jac: [ 5.684e-06  0.000e+00]
     nfev: 45
     njev: 15
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>

目标函数的值基本相同，并且 PDF 无法区分。

import matplotlib.pyplot as plt
x = np.linspace(eps, 15, 300)

c_, scale_ = res_mle.x
X = Weibull(c=c_)*scale_
plt.plot(x, X.pdf(x), '-', label='numerical optimization')

c_, _, scale_ = stats.weibull_min.fit(data, floc=0)
Y = Weibull(c=c_)*scale_
plt.plot(x, Y.pdf(x), '--', label='`weibull_min.fit`')

plt.hist(data, bins=np.linspace(0, 20, 30), density=True, alpha=0.1)
plt.xlabel('x')
plt.ylabel('pdf(x)')
plt.legend()
plt.ylim(0, 0.5)
plt.show()

但是，使用这种方法，可以很容易地更改拟合过程以满足我们的需求，从而实现除 rv_continuous.fit 提供的那些之外的估计程序。例如，我们可以通过将目标函数更改为如下的 **n**egative **l**og **p**roduct **s**pacing 来执行 最大间距估计 (MSE)。

def nlps(x):
    c, scale = x
    X = Weibull(c=c) * scale
    x = np.sort(np.concatenate((data, X.support())))  # Append the endpoints of the support to the data
    return -X.logcdf(x[:-1], x[1:]).sum().real  # Minimize the sum of the logs the probability mass between points

res_mps = optimize.minimize(nlps, x0, bounds=bounds)
res_mps.x

array([0.55917917, 3.85552195])

对于 L 矩的方法（试图匹配分布和样本 L 矩）

def lmoment_residual(x):
    c, scale = x
    X = Weibull(c=c) * scale
    E11 = stats.order_statistic(X, r=1, n=1).mean()
    E12, E22 = stats.order_statistic(X, r=[1, 2], n=2).mean()
    lmoments_X = [E11, 0.5*(E22 - E12)]  # the first two l-moments of the distribution
    lmoments_x = stats.lmoment(data, order=[1, 2])  # first two l-moments of the data
    return np.linalg.norm(lmoments_x - lmoments_X)  # Minimize the norm of the difference

x0 = [0.4, 3]  # This method is a bit sensitive to the initial guess
res_lmom = optimize.minimize(lmoment_residual, x0, bounds=bounds)
res_lmom.x, res_lmom.fun  # residual should be ~0

(array([0.60943275, 3.90234501]), 5.058079590921588e-07)

所有三种方法的图看起来相似，并且似乎描述了数据，这让我们有信心我们获得了良好的拟合。

import matplotlib.pyplot as plt
x = np.linspace(eps, 10, 300)

for res, label in [(res_mle, "MLE"), (res_mps, "MPS"), (res_lmom, "L-moments")]:
    c_, scale_ = res.x
    X = Weibull(c=c_)*scale_
    plt.plot(x, X.pdf(x), '-', label=label)

plt.hist(data, bins=np.linspace(0, 10, 30), density=True, alpha=0.1)
plt.xlabel('x')
plt.ylabel('pdf(x)')
plt.legend()
plt.ylim(0, 0.3)
plt.show()

估计并不总是涉及将分布拟合到数据。例如，在 gh-12134 中，用户需要计算 Weibull 分布的形状和尺度参数，以实现给定的均值和标准差。这很容易适应相同的框架。

moments_0 = np.asarray([8, 20])  # desired mean and standard deviation
def f(x):
    c, scale = x
    X = Weibull(c=c) * scale
    moments_X = X.mean(), X.standard_deviation()
    return np.linalg.norm(moments_X - moments_0)

bounds.lb = np.asarray([0.1, 0.1])  # the Weibull distribution is problematic for very small c
res = optimize.minimize(f, x0, bounds=bounds, method='slsqp')  # easily change the minimization method
res

 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 4.595040620463497e-07
       x: [ 4.627e-01  3.430e+00]
     nit: 21
     jac: [ 1.322e+02 -2.328e-01]
    nfev: 100
    njev: 21

新的和高级功能

转换

数学转换可以应用于随机变量。例如，随机变量和标量之间的许多基本算术运算（+、-、*、/、**）都可以工作。

例如，我们可以看到，Weibull 分布的 RV 的倒数是根据 scipy.stats.invweibull 分布的。（“逆”分布通常是随机变量倒数的基础分布。）

c = 10.6

X = Weibull(c=10.6)  
Y = 1 / X  # compare to `invweibull`
Y.plot()

x = np.linspace(0.8, 2.05, 300)
plt.plot(x, stats.invweibull(c=c).pdf(x), '--')
plt.legend(['`1 / X`', '`invweibull`'])
plt.title("Inverse Weibull PDF")
plt.show()

scipy.stats.chi2 描述了正态分布随机变量的平方和。

X = stats.Normal()
Y = X**2  # compare to chi2
Y.plot(t=('x', eps, 5));

x = np.linspace(eps, 5, 300)
plt.plot(x, stats.chi2(df=1).pdf(x), '--')
plt.ylim(0, 3)
plt.legend(['`X**2`', '`chi2`'])
plt.title("Chi-squared PDF")
plt.show()

scipy.stats.foldcauchy 描述了柯西分布随机变量的绝对值。（“折叠”分布是随机变量绝对值的基础分布。）

Cauchy = stats.make_distribution(stats.cauchy)
c = 4.72

X = Cauchy() + c  
Y = abs(X)  # compare to `foldcauchy`
Y.plot(t=('x', 0, 60))

x = np.linspace(0, 60, 300)
plt.plot(x, stats.foldcauchy(c=c).pdf(x), '--')
plt.legend(['`abs(X)`', '`foldcauchy`'])
plt.title("Folded Cauchy PDF")
plt.show()

scipy.stats.lognormal 描述了正态分布随机变量的指数。之所以这样命名，是因为所得随机变量的对数是正态分布的。（通常，“对数”分布通常是随机变量*指数*的基础分布。）

u, s = 1, 0.5

X = stats.Normal()*s + u
Y = stats.exp(X)  # compare to `lognorm`
Y.plot(t=('x', eps, 9))

x = np.linspace(eps, 9, 300)
plt.plot(x, stats.lognorm(s=s, scale=np.exp(u)).pdf(x), '--')
plt.legend(['`exp(X)`', '`lognorm`'])
plt.title("Log-Normal PDF")
plt.show()

scipy.stats.loggamma 描述了伽玛分布随机变量的对数。请注意，按照典型的惯例，此分布的适当名称应为 exp-gamma，因为 RV 的指数是伽玛分布的。

Gamma = stats.make_distribution(stats.gamma)
a = 0.414

X = Gamma(a=a)  
Y = stats.log(X)  # compare to `loggamma`
Y.plot()

x = np.linspace(-17.5, 2, 300)
plt.plot(x, stats.loggamma(c=a).pdf(x), '--')
plt.legend(['`log(X)`', '`loggamma`'])
plt.title("Exp-Gamma PDF")
plt.show()

scipy.stats.truncnorm 是截断正态随机变量的基础分布。请注意，截断转换可以在平移和缩放正态分布随机变量之前或之后应用，这比 scipy.stats.truncnorm（它本质上在截断*之后*平移和缩放）更容易获得期望的结果。

a, b = 0.1, 2

X = stats.Normal()  
Y = stats.truncate(X, a, b)  # compare to `truncnorm`
Y.plot()

x = np.linspace(a, b, 300)
plt.plot(x, stats.truncnorm(a, b).pdf(x), '--')
plt.legend(['`truncate(X, a, b)`', '`truncnorm`'])
plt.title("Truncated Normal PDF")
plt.show()

scipy.stats.dgamma 是两个伽玛分布的 RV 的混合，其中一个围绕原点反射。（通常，“双”分布是 RV 的混合的基础分布，其中一个被反射。）

a = 1.1
X = Gamma(a=a)
Y = stats.Mixture((X, -X), weights=[0.5, 0.5])
# plot method not available for mixtures

x = np.linspace(-4, 4, 300)
plt.plot(x, Y.pdf(x))
plt.plot(x, stats.dgamma(a=a).pdf(x), '--')
plt.legend(['`Mixture(X, -X)`', '`dgamma`'])
plt.title("Double Gammma PDF")
plt.show()

限制

虽然支持随机变量与 Python 标量或 NumPy 数组之间的大多数算术转换，但有一些限制。

• 仅当参数为正整数时，才实现将随机变量提升为参数的幂。

• 仅当参数为 1 以外的正标量时，才实现将参数提升为随机变量的幂。

• 仅当支持完全是非负数或非正数时，才实现除以随机变量。

• 仅当支持是非负数时，才实现随机变量的对数。（负值的对数是虚数。）

两个随机变量之间的算术运算*尚未*支持。请注意，此类运算在数学上要复杂得多；例如，两个随机变量之和的 PDF 涉及两个 PDF 的卷积。

此外，虽然转换是可组合的，但 a) 截断和 b) 平移/缩放都只能执行一次。例如，Y = 2 * stats.exp(X + 1) 将引发错误，因为这将需要在指数运算之前平移和在指数运算之后缩放；这被视为“平移/缩放”两次。但是，这里（在许多情况下）可以进行数学简化来避免该问题：Y = (2*math.exp(2)) * stats.exp(X) 是等效的，并且只需要一个缩放操作。

尽管转换相当稳健，但它们都依赖于通用的实现，这可能会导致数值困难。如果您担心转换结果的准确性，请考虑将生成的 PDF 与随机样本的直方图进行比较。

X = stats.Normal()
Y = X**3
x = np.linspace(-5, 5, 300)
plt.plot(x, Y.pdf(x), label='pdf')
plt.hist(X.sample(100000)**3, density=True, bins=np.linspace(-5, 5, 100), alpha=0.5);
plt.ylim(0, 2)
plt.show()

准蒙特卡洛采样

随机变量可以从统计分布生成准随机、低差异的样本。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
X = stats.Normal()

rng = np.random.default_rng(7824387278234)
qrng = stats.qmc.Sobol(1, rng=rng)  # instantiate a QMCEngine

bins = np.linspace(-3.5, 3.5, 31)
plt.hist(X.sample(512, rng=qrng), bins, alpha=0.5, label='quasi-random')
plt.hist(X.sample(512, rng=rng), bins, alpha=0.5, label='pseudo-random')
plt.title('Histogram of normally-distributed sample')
plt.legend()
plt.show()

请注意，在生成多个样本时（例如，len(shape) > 1），沿着零轴的每个切片都是独立的、低差异的序列。以下生成两个独立的 QMC 样本，每个长度为 512。

samples = X.sample((512, 2), rng=qrng)
plt.hist(samples[:, 0], bins, alpha=0.5, label='sample 0')
plt.hist(samples[:, 1], bins, alpha=0.5, label='sample 1')
plt.title('Histograms of normally-distributed samples')
plt.legend()
plt.show()

如果我们生成 512 个独立的 QMC 序列，每个长度为 2，则结果会大相径庭。

samples = X.sample((2, 512), rng=qrng)
plt.hist(samples[0], bins, alpha=0.5, label='sample 0')
plt.hist(samples[1], bins, alpha=0.5, label='sample 1')
plt.title('Histograms of normally-distributed samples')
plt.legend()
plt.show()

精度

对于某些分布，例如 scipy.stats.norm，几乎所有分布的方法都经过自定义，以确保准确的计算。对于其他分布，例如 scipy.stats.gausshyper，定义的只有 PDF，其他方法必须基于 PDF 以数值方式计算。例如，高斯超几何分布的生存函数（互补 CDF）是通过将 PDF 从 0（支持的左端）数值积分到 x 来获得 CDF，然后取补数来计算的。

a, b, c, z = 1.5, 2.5, 2, 0
x = 0.5
frozen = stats.gausshyper(a=a, b=b, c=c, z=z)
frozen.sf(x)

0.28779340921080643

frozen.sf(x) == 1 - integrate.quad(frozen.pdf, 0, x)[0]

True

但是，另一种方法是将 PDF 从 x 数值积分到 1（支持的右端）。

integrate.quad(frozen.pdf, x, 1)[0]

0.28779340919669805

这些是不同的但同样有效的方法，因此，假设 PDF 是准确的，则两个结果不可能以相同的方式不准确。因此，我们可以通过比较它们来估计结果的准确性。

res1 = frozen.sf(x)
res2 = integrate.quad(frozen.pdf, x, 1)[0]
abs((res1 - res2) / res1)

4.902259981810363e-11

新的基础设施知道计算大多数量的几种不同方式。例如，下图说明了各种分布函数之间的关系。

它遵循决策树来选择最准确的估计量的方法。这些用于确定“最佳”方法的特定决策树可能会在 SciPy 的版本之间发生变化，但是用于计算互补 CDF 的示例可能类似于

和用于计算矩的示例

但是，您可以使用 method 参数覆盖它使用的方法。

GaussHyper = stats.make_distribution(stats.gausshyper)
X = GaussHyper(a=a, b=b, c=c, z=z)

对于 ccdf，method='quadrature' 在右尾积分 PDF。

X.ccdf(x, method='quadrature') == integrate.tanhsinh(X.pdf, x, 1).integral

True

method='complement' 取 CDF 的补数。

X.ccdf(x, method='complement') == (1 - integrate.tanhsinh(X.pdf, 0, x).integral)

True

method='logexp' 取 CCDF 对数的指数。

X.ccdf(x, method='logexp') == np.exp(integrate.tanhsinh(lambda x: X.logpdf(x), x, 1, log=True).integral)

True

如有疑问，请考虑尝试不同的 method 来计算给定的量。

性能

考虑一下涉及高斯超几何分布的计算的性能。对于标量形状和参数，新旧基础设施的性能相当。新基础设施预计不会快得多；尽管它减少了参数验证的“开销”，但在标量量的数值积分方面并没有明显更快，而后者在这里占主导地位的执行时间。

%timeit X.cdf(x)  # new infrastructure

402 µs ± 868 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%timeit stats.gausshyper.cdf(x, a, b, c, z)  # old infrastructure

392 µs ± 955 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

np.isclose(X.cdf(x), stats.gausshyper.cdf(x, a, b, c, z))

True

但当涉及数组时，新基础设施的速度要快得多。这是因为新基础设施的底层积分器（scipy.integrate.tanhsinh）和求根器（scipy.optimize.elementwise.find_root）是原生向量化的，而旧基础设施使用的传统例程（scipy.integrate.quad 和 scipy.optimize.brentq）则不是。

x = np.linspace(0, 1, 1000)

%timeit X.cdf(x)  # new infrastructure

2.47 ms ± 47.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit stats.gausshyper.cdf(x, a, b, c, z)  # old infrastructure

330 ms ± 150 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

对于通过数值积分或反演（求根）计算底层函数的数组计算，可以预期有类似的性能提升。

可视化

我们可以使用方便的方法 plot 轻松可视化随机变量的底层分布函数。

import matplotlib.pyplot as plt
ax = X.plot()
plt.show()

X.plot(y='cdf')
plt.show()

plot 方法非常灵活，其签名灵感来源于 《The Grammar of Graphics》。例如，对于 \(x\) 在 \([-10, 10]\) 上的情况，绘制 PDF 与 CDF 的对比图。

X.plot('cdf', 'pdf', t=('x', -10, 10))
plt.show()

顺序统计量

如上面的“拟合”部分所示，现在支持来自分布的随机样本的顺序统计量的分布。例如，我们可以绘制样本大小为 4 的正态分布的顺序统计量的概率密度函数。

n = 4
r = np.arange(1, n+1)
X = stats.Normal()
Y = stats.order_statistic(X, r=r, n=n)
Y.plot()
plt.show()

计算这些顺序统计量的期望值

Y.mean()

array([-1.02937537, -0.29701138,  0.29701138,  1.02937537])

结论

尽管改变可能让人感到不适，但新的基础设施为改进 SciPy 概率分布用户体验铺平了道路。当然，它并不完整！在编写本指南时，计划在未来的版本中添加以下功能

• 额外的内置分布（除了 Normal 和 Uniform 之外）

• 离散分布

• 循环分布，包括将任意连续分布环绕在圆上的能力

• 除了 NumPy 之外，还支持其他与 Python Array API 兼容的后端（例如 CuPy、PyTorch、JAX）

与往常一样，SciPy 的开发是由 SciPy 社区为 SciPy 社区完成的。我们始终欢迎关于其他功能的建议和帮助实现它们的帮助。

我们希望本指南能够激发您尝试新基础设施的优势并简化这一过程。