随机变量转换指南#
背景#
在 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)
np.float64(0.24197072451914337)
较不常用的方法是调用分布对象的 __call__ 方法，无论原始类是什么，它都会返回一个 rv_continuous_frozen 的实例。
frozen = stats.norm()
type(frozen)
scipy.stats._distn_infrastructure.rv_continuous_frozen
这个新对象的方法仅接受自变量，不接受分布参数。
frozen.pdf(x)
np.float64(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=np.float64(0.0), sigma=np.float64(1.0))
实例化后，形状参数可以作为属性读取（但不能写入）。
X.mu, X.sigma
(np.float64(0.0), np.float64(1.0))
scipy.stats.Normal 的文档包含其每个方法的详细文档链接。（相比之下，可以查看 scipy.stats.norm 的文档。）
注：虽然截至 SciPy 1.15，只有 Normal 和 Uniform 拥有渲染好的文档并能直接从 stats 导入，但几乎所有其他旧分布都可以通过 make_distribution（稍后讨论）在新接口中使用。
然后通过仅传递自变量（如果有）来调用方法，而不传递分布参数。
X.mean()
np.float64(0.0)
X.pdf(x)
np.float64(0.24197072451914337)
对于这类简单的调用（例如自变量是有效的浮点数），调用新随机变量的方法通常会比调用旧分布的同类方法更快。
%timeit X.pdf(x)
2.26 μs ± 3.3 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
%timeit frozen.pdf(x)
17.7 μs ± 80.6 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
%timeit dist.pdf(x, loc=mu, scale=sigma)
17.7 μs ± 20.6 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
请注意，调用 X.pdf 和 frozen.pdf 的方式是完全相同的，而调用 dist.pdf 也非常相似——唯一的区别是调用 dist.pdf 时包含形状参数，而在新架构中，形状参数仅在随机变量实例化时提供。
除了 pdf 和 mean 之外，新架构中的其他几个方法与旧方法基本相同。
logpdf（概率密度函数的对数）cdf（累积分布函数）logcdf（累积分布函数的对数）entropy（微分熵）median（中位数）support (支撑集)
其他方法有新的名称，但在其他方面是直接替代品。
sf（生存函数）\(\rightarrow\)ccdf（互补累积分布函数）logsf\(\rightarrow\)logccdfppf（百分点函数）\(\rightarrow\)icdf（逆累积分布函数）isf（逆生存函数）\(\rightarrow\)iccdf（逆互补累积分布函数）std\(\rightarrow\)standard_deviation（标准差）var\(\rightarrow\)variance（方差）
新架构还增加了一些与上述方法类似的新方法。
ilogcdf（累积分布函数对数的逆函数）ilogccdf（互补累积分布函数对数的逆函数）logentropy（熵的对数）mode（分布的众数）skewness（偏度）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 `<class 'scipy.stats._distribution_infrastructure.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 p
根据 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
(np.float64(0.8862269254527579), np.float64(0.8862269254527579))
然而，旧架构还有一个 stats 方法，它提供了分布的各种统计量。以下统计量与指示的字符相关联。
均值（关于原点的一阶原点矩，
'm'）方差（二阶中心矩，
'v'）偏度（三阶标准化矩，
's'）超额峰度（四阶标准化矩减去 3，
'k'）
例如
stats.weibull_min.stats(c, moments='mvsk')
(np.float64(0.8862269254527579),
np.float64(0.21460183660255183),
np.float64(0.6311106578189344),
np.float64(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'))
(np.float64(0.8862269254527579),
np.float64(0.21460183660255183),
np.float64(0.6311106578189344),
np.float64(3.2450893006876314))
请注意 峰度 (kurtosis) 定义的差异。以前提供的是“超额”（又称“Fisher”）峰度。根据惯例（而非自然数学定义），这是标准化的四阶矩平移一个常数值 (3)，使得正态分布的值为 0.0。
stats.norm.stats(moments='k')
np.float64(0.0)
新的 moment 和 kurtosis 方法默认不遵循这一惯例；它们根据数学定义报告标准化四阶矩。这也被称为“非超额”（或“Pearson”）峰度，对于标准正态分布，其值为 3.0。
stats.Normal().moment(4, kind='standardized')
np.float64(3.0)
为了方便起见，kurtosis 方法提供了一个 convention 参数来在两者之间进行选择。
stats.Normal().kurtosis(convention='excess'), stats.Normal().kurtosis(convention='non-excess')
(np.float64(0.0), np.float64(3.0))
随机样本#
旧架构的 rvs 方法已被 sample 取代，但在 1) 需要控制随机状态或 2) 形状、位置或缩放参数使用数组时，应注意两个区别。
首先，根据 SPEC 7，参数 random_state 被 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)
(np.float64(1.6243453636632417), np.float64(1.6243453636632417))
在新架构中，向 sample 的 rng 参数传递 numpy.random.Generator 实例来控制随机状态。请注意，整数种子现在被转换为 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)
(np.float64(0.345584192064786), np.float64(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([[-0.74566906, -0.99838065, 0.87051894, -1.72507176],
[-0.09593492, -1.25679529, -0.07644145, -0.61029409],
[ 0.62681897, 1.68828187, 0.17886689, 0.34121354]])
除了名称上的差异外，在涉及数组形状参数时，行为也有所不同。以前，分布参数数组的形状必须包含在 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)
(np.float64(-1.959963984540054), np.float64(1.959963984540054))
现在，请使用所需的每侧尾部概率质量来调用逆 CDF 和互补 CDF 方法。
p = 1 - a
X.icdf(p/2), X.iccdf(p/2)
(np.float64(-1.959963984540054), np.float64(1.959963984540054))
似然函数#
旧架构提供了一个计算负对数似然函数的函数，错误地命名为 nnlf（应为 nllf）。它接受分布参数作为元组，并将观测值作为数组。
mu = 0
sigma = 1
data = stats.norm.rvs(size=100, loc=mu, scale=sigma)
stats.norm.nnlf((mu, sigma), data)
np.float64(131.07116392314364)
现在，只需根据其数学定义计算负对数似然即可。
X = stats.Normal(mu=mu, sigma=sigma)
-X.logpdf(data).sum()
np.float64(131.07116392314364)
期望值#
旧架构的 expect 方法通过分布的 PDF 加权来估计任意函数的定积分。例如，分布的四阶矩由下式给出
def f(x): return x**4
stats.norm.expect(f, lb=-np.inf, ub=np.inf)
np.float64(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)
np.float64(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
np.float64(1.6578138621431193)
以及使用新的随机变量
integrate.quad(f, a=a, b=b)[0] / X.cdf(a, b)
np.float64(1.6578138621431193)
请注意，新随机变量的 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)
(np.float64(1.8268432669973347), np.float64(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)
(np.float64(1.8268432669973347), np.float64(1.9640905288934327))
请注意，此示例中的估计结果非常糟糕，启发式方法的糟糕性能也是决定不提供替代方案的原因之一。
最大似然估计#
旧架构的最后一个方法是 fit，它在给定分布样本的情况下估计底层分布的位置、尺度和形状参数。
c_, loc_, scale_ = stats.weibull_min.fit(data)
c_, loc_, scale_
(np.float64(0.5960699537061247),
np.float64(1.8950104942593064e-05),
np.float64(4.015875612107546))
这种方法的便利性既是福音也是诅咒。
当它奏效时，其简洁性备受赞赏。对于某些分布，参数的最大似然估计 (MLE) 的解析表达式已经过手动编程，对于这些分布，fit 方法非常可靠。
然而，对于大多数分布，拟合是通过负对数似然函数的数值最小化完成的。这不能保证成功——既因为 MLE 固有的局限性，也因为数值优化的技术现状——在这些情况下，fit 的用户会陷入困境。然而，少量的理解对于确保成功大有裨益，因此在这里我们展示一个使用新架构进行拟合的微型教程。
首先，请注意 MLE 是将分布拟合到数据的几种潜在策略之一。它只不过是找到使负对数似然 (NLLF) 最小化的分布参数值。请注意，对于生成数据的分布，NLLF 是
stats.weibull_min.nnlf((c, loc, scale), data)
np.float64(246.9503871545133)
使用 fit 估计的参数的 NLLF 更低
stats.weibull_min.nnlf((c_, loc_, scale_), data)
np.float64(230.35385152369787)
因此，fit 认为它的工作完成了，无论这是否满足用户对“良好拟合”的想法，或者参数是否在合理的范围内。除了允许用户指定分布参数的猜测值外，它对过程几乎没有控制权。（注：技术上，optimizer 参数可以用于更好的控制，但这没有很好的文档记录，而且它的使用抵消了 fit 方法提供的任何便利。）
另一个问题是，MLE 对于某些分布而言天生表现不佳。例如在这个例子中，只需将位置设置为数据的最小值，我们就可以让 NLLF 变得任意低——即使对于形状和尺度参数的荒谬值也是如此。
stats.weibull_min.nnlf((0.1, np.min(data), 10), data)
np.float64(-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_
(np.float64(0.5834128280886433), 0, np.float64(3.829972848420729))
stats.weibull_min.nnlf((c_, loc_, scale_), data)
np.float64(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_))
np.float64(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.96174292277612
x: [ 5.834e-01 3.830e+00]
nit: 13
jac: [ 2.842e-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 未提供的估计过程。例如，我们可以通过将客观函数更改为负对数乘积间距来执行最大间距估计 (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.55917929, 3.85553545])
对于 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]), np.float64(5.036804210160496e-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.5947808997191174e-07
x: [ 4.627e-01 3.430e+00]
nit: 21
jac: [ 1.322e+02 -2.329e-01]
nfev: 100
njev: 21
multipliers: []
新功能和高级功能#
转换#
数学变换可以应用于随机变量。例如，随机变量与标量之间的许多基本算术运算（+, -, *, /, **）都可以工作。
例如，我们可以看到 Weibull 分布随机变量的倒数遵循 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()
/home/treddy/github_projects/scipy/doc/build/env/lib/python3.11/site-packages/scipy/stats/_distribution_infrastructure.py:1974: RuntimeWarning: divide by zero encountered in divide
h=lambda u: 1 / u, dh=lambda u: 1 / u ** 2)
/home/treddy/github_projects/scipy/doc/build/env/lib/python3.11/site-packages/scipy/stats/_continuous_distns.py:2760: RuntimeWarning: invalid value encountered in multiply
return c*pow(x, c-1)*np.exp(-pow(x, c))
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 分布随机变量的绝对值。（“折叠”分布是随机变量绝对值所遵循的分布。）
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 描述了 Gamma 分布随机变量的对数。请注意，根据典型惯例，该分布的适当名称应该是 exp-gamma，因为随机变量的指数是 gamma 分布的。
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 是两个 gamma 分布随机变量的混合，其中一个关于原点对称。（通常，“双”分布是随机变量及其对称混合后所遵循的分布。）
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），沿第 0 轴的每个切片都是独立的低差异序列。以下生成两个独立的 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 进行数值计算。例如，Gauss 超几何分布的生存函数（互补 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)
np.float64(0.28779340921080643)
frozen.sf(x) == 1 - integrate.quad(frozen.pdf, 0, x)[0]
np.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)
np.float64(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
np.True_
method='complement' 取 CDF 的补集。
X.ccdf(x, method='complement') == (1 - integrate.tanhsinh(X.pdf, 0, x).integral)
np.True_
method='logexp' 对 CCDF 的对数取指数。
X.ccdf(x, method='logexp') == np.exp(integrate.tanhsinh(lambda x: X.logpdf(x), x, 1, log=True).integral)
np.True_
如有疑问，请考虑尝试不同的 method 来计算给定量。
性能#
考虑涉及 Gauss 超几何分布的计算性能。对于标量形状和自变量，新旧架构的性能相当。新架构预计不会快很多；虽然它减少了参数验证的“开销”，但在标量的数值积分方面并没有显著加快，而后者在这里占据了执行时间的大部分。
%timeit X.cdf(x) # new infrastructure
389 μs ± 1.33 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
%timeit stats.gausshyper.cdf(x, a, b, c, z) # old infrastructure
336 μs ± 6.88 μs 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))
np.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.35 ms ± 1.51 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit stats.gausshyper.cdf(x, a, b, c, z) # old infrastructure
294 ms ± 247 μ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)。例如，在 \([-10, 10]\) 上使用自变量 \(x\)，绘制 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 社区完成的。我们随时欢迎对其他功能的建议以及对实现这些功能的帮助。
我们希望本指南能阐明尝试新架构的优势，并简化这一过程。