SciPy 1.12 中文文档（三）
稀疏数组（scipy.sparse）
简介
scipy.sparse 及其子模块提供了用于处理稀疏数组的工具。稀疏数组是只有数组中少数位置包含任何数据的数组，大多数位置被视为“空”。稀疏数组很有用，因为它们允许用于线性代数（scipy.sparse.linalg）或基于图的计算（scipy.sparse.csgraph）的算法更简单、更快速或内存消耗较少，但是它们通常在像切片、重塑或赋值等操作上不太灵活。本指南将介绍scipy.sparse中稀疏数组的基础知识，解释稀疏数据结构的独特之处，并引导用户查看用户指南中解释稀疏线性代数和图方法的其他部分。
入门稀疏数组
稀疏数组是一种特殊类型的数组，其中数组中只有少数位置包含数据。这允许使用压缩表示数据，仅记录存在数据的位置。有许多不同的稀疏数组格式，每种格式在压缩和功能之间进行不同的权衡。首先，让我们构建一个非常简单的稀疏数组，即坐标（COO）数组 (coo_array) 并将其与密集数组进行比较：
>>> import scipy as sp
>>> import numpy
>>> dense = numpy.array([[1, 0, 0, 2], [0, 4, 1, 0], [0, 0, 5, 0]])
>>> sparse = sp.sparse.coo_array(dense)
>>> dense
array([[1, 0, 0, 2],
[0, 4, 1, 0],
[0, 0, 5, 0]])
>>> sparse
<3x4 sparse array of type '<class 'numpy.int64'>'
with 5 stored elements in COOrdinate format>
注意，在我们的密集数组中，有五个非零值。例如，2 在位置 0,3，4 在位置 1,1。所有其他值都为零。稀疏数组显式记录这五个值（参见 COOrdinate format 中的 5 个存储元素），然后将所有其余的零表示为隐式值。
大多数稀疏数组方法的工作方式与密集数组方法类似：
>>> sparse.max()
5
>>> dense.max()
5
>>> sparse.argmax()
10
>>> dense.argmax()
10
>>> sparse.mean()
1.0833333333333333
>>> dense.mean()
1.0833333333333333
稀疏数组还具有一些“额外”的属性，例如 .nnz，它返回存储值的数量：
>>> sparse.nnz
5
大多数减少操作，例如 .mean()、.sum() 或 .max()，在应用到稀疏数组的轴上时将返回一个 numpy 数组：
>>> sparse.mean(axis=1)
array([0.75, 1.25, 1.25])
这是因为稀疏数组上的减少操作通常是密集的。
理解稀疏数组格式
不同类型的稀疏数组具有不同的功能。例如，COO 数组不能被索引或切片：
>>> dense[2, 2]
5
>>> sparse[2, 2]
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: 'coo_array' object is not subscriptable
但是，其他格式，例如压缩稀疏行（CSR）csr_array 支持切片和元素索引：
>>> sparse.tocsr()[2, 2]
5
有时，scipy.sparse会返回与输入稀疏矩阵格式不同的稀疏矩阵格式。例如，COO 格式的两个稀疏数组的点积将是 CSR 格式数组：
>>> sparse @ sparse.T
<3x3 sparse array of type '<class 'numpy.int64'>'
with 5 stored elements in Compressed Sparse Row format>
这种改变是因为scipy.sparse会改变输入稀疏数组的格式，以使用最有效的计算方法。
scipy.sparse模块包含以下格式，每种格式都有自己独特的优势和劣势：
-
块状稀疏行（BSR）数组
scipy.sparse.bsr_array，在数组的数据部分以连续的块出现时最合适。 -
坐标（COO）数组
scipy.sparse.coo_array提供了一种简单的构建稀疏数组和原地修改它们的方法。COO 也可以快速转换为其他格式，如 CSR、CSC 或 BSR。 -
压缩稀疏行（CSR）数组
scipy.sparse.csr_array，最适用于快速算术运算、向量乘积和按行切片。 -
压缩稀疏列（CSC）数组
scipy.sparse.csc_array最适用于快速算术运算、向量乘积和按列切片。 -
对角线（DIA）数组
scipy.sparse.dia_array，对于有效存储和快速算术运算很有用，只要数据主要沿着数组的对角线出现。 -
键值字典（DOK）数组
scipy.sparse.dok_array，对于快速构建和单元素访问很有用。 -
列表列表（LIL）数组
scipy.sparse.lil_array，对于快速构建和修改稀疏数组很有用。
每种稀疏数组格式的优势和劣势的更多信息可以在它们的文档中找到。
所有scipy.sparse数组格式都可以直接从numpy.ndarray构建。然而，某些稀疏格式也可以以不同的方式构建。每个稀疏数组格式都有不同的优势，并且这些优势在每个类中都有文档记录。例如，构建稀疏数组最常见的方法之一是从单独的row、column和data值构建稀疏数组。对于我们之前的数组：
>>> dense
array([[1, 0, 0, 2],
[0, 4, 1, 0],
[0, 0, 5, 0]])
row、column和data数组描述了稀疏数组中条目的行、列和值：
>>> row = [0,0,1,1,2]
>>> col = [0,3,1,2,2]
>>> data = [1,2,4,1,5]
使用这些，我们现在可以定义一个稀疏数组而不需要首先构建一个密集数组：
>>> csr = sp.sparse.csr_array((data, (row, col)))
>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
with 5 stored elements in Compressed Sparse Row format>
不同的类有不同的构造函数，但是scipy.sparse.csr_array、scipy.sparse.csc_array和scipy.sparse.coo_array允许使用这种构造方式。
稀疏数组，隐式零和重复
稀疏数组很有用，因为它们表示大部分值是隐式的，而不是存储一个实际的占位值。在scipy.sparse中，用于表示“无数据”的值是隐式零。当需要显式零时，这可能会令人困惑。例如，在图方法中，我们经常需要能够区分（A）连接节点 i 和 j 的权重为零的链接和（B）i 和 j 之间没有链接。稀疏矩阵可以做到这一点，只要我们记住显式和隐式零即可。
例如，在我们之前的csr数组中，我们可以通过将它包含在data列表中来包含一个显式零。让我们把底行和最后一列数组中的最后一个条目视为显式零：
>>> row = [0,0,1,1,2,2]
>>> col = [0,3,1,2,2,3]
>>> data = [1,2,4,1,5,0]
那么，我们的稀疏数组将有六个存储元素，而不是五个：
>>> csr = sp.sparse.csr_array((data, (row, col)))
>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
with 6 stored elements in Compressed Sparse Row format>
“额外”的元素就是我们的显式零。当转换回密集数组时，两者仍然是相同的，因为密集数组将所有东西都显式地表示：
>>> csr.todense()
array([[1, 0, 0, 2],
[0, 4, 1, 0],
[0, 0, 5, 0]])
>>> dense
array([[1, 0, 0, 2],
[0, 4, 1, 0],
[0, 0, 5, 0]])
但是，对于稀疏算术、线性代数和图方法，位置2,3处的值将被视为显式零。要去除此显式零，我们可以使用csr.eliminate_zeros()方法。这个方法在稀疏数组中原地操作，并移除任何零值存储元素：
>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
with 6 stored elements in Compressed Sparse Row format>
>>> csr.eliminate_zeros()
>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
with 5 stored elements in Compressed Sparse Row format>
在csr.eliminate_zeros()之前，有六个存储元素。之后，只有五个存储元素。
另一个复杂性点源于在构建稀疏数组时处理 重复项 的方式。当我们在构建稀疏数组时在 row,col 处有两个或更多条目时，就会发生 重复项。例如，我们可能用重复值在 1,1 处表示先前的数组：
>>> row = [0,0,1,1,1,2]
>>> col = [0,3,1,1,2,2]
>>> data = [1,2,1,3,1,5]
在这种情况下，我们可以看到有 两个 data 值对应于我们最终数组中的 1,1 位置。scipy.sparse 将单独存储这些值：
>>> dupes = sp.sparse.coo_array((data, (row, col)))
>>> dupes
<3x4 sparse array of type '<class 'numpy.int64'>'
with 6 stored elements in COOrdinate format>
请注意，这个稀疏数组中有六个存储的元素，尽管只有五个唯一的数据位置。当这些数组转换回密集数组时，重复值将被求和。因此，在位置 1,1 处，密集数组将包含重复存储条目的总和，即 1 + 3：
>>> dupes.todense()
array([[1, 0, 0, 2],
[0, 4, 1, 0],
[0, 0, 5, 0]])
要删除稀疏数组本身中的重复值，从而减少存储元素的数量，可以使用 .sum_duplicates() 方法：
>>> dupes.sum_duplicates()
>>> dupes
<3x4 sparse array of type '<class 'numpy.int64'>'
with 5 stored elements in COOrdinate format>
现在我们的稀疏数组中只有五个存储的元素，且与本指南中一直使用的数组相同：
>>> dupes.todense()
array([[1, 0, 0, 2],
[0, 4, 1, 0],
[0, 0, 5, 0]])
规范格式
几种稀疏数组格式具有“规范格式”，以实现更高效的操作。通常这些格式包括像增加限制这样的额外条件：
-
任何值都没有重复条目
-
已排序的索引
具有规范形式的类包括：coo_array，csr_array，csc_array 和 bsr_array。详细信息请参阅这些类的文档字符串，了解每种规范表示的细节。
要检查这些类的实例是否处于规范形式，请使用 .has_canonical_format 属性：
>>> coo = sp.sparse.coo_array(([1, 1, 1], ([0, 2, 1], [0, 1, 2])))
>>> coo.has_canonical_format
False
要将实例转换为规范形式，请使用 .sum_duplicates() 方法：
>>> coo.sum_duplicates()
>>> coo.has_canonical_format
True
稀疏数组的下一步操作
当处理大型、几乎为空的数组时，稀疏数组类型最为有用。特别是在这些情况下，稀疏线性代数和稀疏图方法的效率显著提高。
使用 ARPACK 解决稀疏特征值问题
介绍
ARPACK [1] 是一个 Fortran 包，提供快速查找大稀疏矩阵少数特征值/特征向量的例程。为了找到这些解，它只需左乘问题矩阵。此操作通过反向通信接口执行。这种结构的结果是 ARPACK 能够找到任何线性函数映射向量到向量的特征值和特征向量。
ARPACK 提供的所有功能都包含在两个高级接口中 scipy.sparse.linalg.eigs 和 scipy.sparse.linalg.eigsh。 eigs 提供实现接口，用于查找实数或复数非对称方阵的特征值/特征向量，而 eigsh 提供接口，用于实对称或复共轭矩阵。
基本功能
ARPACK 可以解决形如
[A \mathbf{x} = \lambda \mathbf{x}]
或一般的特征值问题形式
[A \mathbf{x} = \lambda M \mathbf{x}.]
ARPACK 的优势在于它可以计算特定子集的特征值/特征向量对。这是通过关键字which实现的。以下which值可用：
-
which = 'LM': 最大幅度的特征值 (eigs,eigsh), 即复数的欧几里德范数中的最大特征值. -
which = 'SM': 最小幅度的特征值 (eigs,eigsh), 即复数的欧几里德范数中的最小特征值. -
which = 'LR': 最大实部的特征值 (eigs). -
which = 'SR': 最小实部的特征值 (eigs). -
which = 'LI': 最大虚部的特征值 (eigs). -
which = 'SI': 最小虚部的特征值 (eigs). -
which = 'LA': 最大代数值的特征值 (eigsh), 即包含任何负号的最大特征值. -
which = 'SA': 最小代数值的特征值 (eigsh), 即包含任何负号的最小特征值. -
which = 'BE': 光谱两端的特征值 (eigsh).
注意，ARPACK 通常更擅长找到极端特征值，即具有较大幅度的特征值。特别是，使用which = 'SM'可能导致执行时间缓慢和/或异常结果。更好的方法是使用转移反演模式。
转移反演模式
Shift-invert mode 依赖于以下观察。对于广义特征值问题
[A \mathbf{x} = \lambda M \mathbf{x},]
可以证明
[(A - \sigma M)^{-1} M \mathbf{x} = \nu \mathbf{x},]
其中
[\nu = \frac{1}{\lambda - \sigma}.]
举例
假设您想要查找大矩阵的最小和最大特征值以及相应的特征向量。ARPACK 可以处理多种输入形式：如密集矩阵，例如numpy.ndarray 实例，稀疏矩阵，例如scipy.sparse.csr_matrix，或者从scipy.sparse.linalg.LinearOperator 派生的一般线性操作员。为了简单起见，在本例中，我们将构建一个对称的正定矩阵。
>>> import numpy as np
>>> from scipy.linalg import eig, eigh
>>> from scipy.sparse.linalg import eigs, eigsh
>>> np.set_printoptions(suppress=True)
>>> rng = np.random.default_rng()
>>>
>>> X = rng.random((100, 100)) - 0.5
>>> X = np.dot(X, X.T) # create a symmetric matrix
现在我们有一个对称矩阵 X，用来测试这些程序。首先，使用 eigh 计算标准特征值分解：
>>> evals_all, evecs_all = eigh(X)
随着 X 的维度增长，这个例程变得非常慢。特别是，如果只需要少量特征向量和特征值，ARPACK 可能是一个更好的选择。首先，计算 X 的最大特征值 (which = 'LM') 并将其与已知结果进行比较：
>>> evals_large, evecs_large = eigsh(X, 3, which='LM')
>>> print(evals_all[-3:])
[29.22435321 30.05590784 30.58591252]
>>> print(evals_large)
[29.22435321 30.05590784 30.58591252]
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1\. 0\. 0.], # may vary (signs)
[ 0\. 1\. 0.],
[-0\. 0\. -1.]])
结果如预期。ARPACK 恢复了所需的特征值，并且它们与先前已知的结果相匹配。此外，特征向量是正交的，这是我们预期的。现在，让我们尝试解最小幅度特征值的问题：
>>> evals_small, evecs_small = eigsh(X, 3, which='SM')
Traceback (most recent call last): # may vary (convergence)
...
scipy.sparse.linalg._eigen.arpack.arpack.ArpackNoConvergence:
ARPACK error -1: No convergence (1001 iterations, 0/3 eigenvectors converged)
糟糕。我们看到，如上所述，ARPACK 并不太擅长找到小特征值。可以通过几种方法来解决这个问题。我们可以增加容差 (tol) 以加快收敛速度：
>>> evals_small, evecs_small = eigsh(X, 3, which='SM', tol=1E-2)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 0.99999999 0.00000024 -0.00000049], # may vary (signs)
[-0.00000023 0.99999999 0.00000056],
[ 0.00000031 -0.00000037 0.99999852]])
这种方法有效，但结果的精度会降低。另一种选择是将最大迭代次数 (maxiter) 从 1000 增加到 5000：
>>> evals_small, evecs_small = eigsh(X, 3, which='SM', maxiter=5000)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1\. 0\. 0.], # may vary (signs)
[-0\. 1\. 0.],
[ 0\. 0\. -1.]])
我们得到了预期的结果，但计算时间要长得多。幸运的是，ARPACK 包含一个模式，可以快速确定非外部特征值：shift-invert mode。如上所述，这种模式涉及将特征值问题转换为具有不同特征值的等价问题。在这种情况下，我们希望找到接近零的特征值，因此我们选择 sigma = 0。然后转换后的特征值将满足 (\nu = 1/(\lambda - \sigma) = 1/\lambda)，因此我们的小特征值 (\lambda) 变为大特征值 (\nu)。
>>> evals_small, evecs_small = eigsh(X, 3, sigma=0, which='LM')
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1\. 0\. 0.], # may vary (signs)
[ 0\. -1\. -0.],
[-0\. -0\. 1.]])
我们得到了我们希望的结果，计算时间大大减少。请注意，从 (\nu \to \lambda) 的转换完全在后台进行。用户不必担心细节。
移位-反转模式提供的不仅仅是获取少量小特征值的快速方法。比如，您希望找到内部特征值和特征向量，例如那些接近(\lambda = 1)的。只需设置sigma = 1，ARPACK 将处理其余部分：
>>> evals_mid, evecs_mid = eigsh(X, 3, sigma=1, which='LM')
>>> i_sort = np.argsort(abs(1. / (1 - evals_all)))[-3:]
>>> evals_all[i_sort]
array([0.94164107, 1.05464515, 0.99090277])
>>> evals_mid
array([0.94164107, 0.99090277, 1.05464515])
>>> print(np.dot(evecs_mid.T, evecs_all[:,i_sort]))
array([[-0\. 1\. 0.], # may vary (signs)
[-0\. -0\. 1.],
[ 1\. 0\. 0.]]
特征值的顺序不同，但它们都在那里。请注意，移位-反转模式需要内部解决矩阵的逆问题。这由eigsh和eigs自动处理，但用户也可以指定该操作。有关详细信息，请参阅scipy.sparse.linalg.eigsh和scipy.sparse.linalg.eigs的文档字符串。
使用 LinearOperator
现在我们考虑一个情况，您想避免创建密集矩阵，而是使用scipy.sparse.linalg.LinearOperator。我们的第一个线性算子应用于输入向量和用户提供给算子本身的向量(\mathbf{d})之间的逐元素乘法。这个算子模拟了一个对角矩阵，其主对角线上的元素是(\mathbf{d})，它的主要优点在于前向和伴随操作都是简单的逐元素乘法，而不是矩阵-向量乘法。对于对角矩阵，我们期望的特征值等于沿主对角线的元素，即(\mathbf{d})。使用eigsh得到的特征值和特征向量与应用于密集矩阵时使用eigh得到的进行比较：
>>> from scipy.sparse.linalg import LinearOperator
>>> class Diagonal(LinearOperator):
... def __init__(self, diag, dtype='float32'):
... self.diag = diag
... self.shape = (len(self.diag), len(self.diag))
... self.dtype = np.dtype(dtype)
... def _matvec(self, x):
... return self.diag*x
... def _rmatvec(self, x):
... return self.diag*x
>>> N = 100
>>> rng = np.random.default_rng()
>>> d = rng.normal(0, 1, N).astype(np.float64)
>>> D = np.diag(d)
>>> Dop = Diagonal(d, dtype=np.float64)
>>> evals_all, evecs_all = eigh(D)
>>> evals_large, evecs_large = eigsh(Dop, 3, which='LA', maxiter=1e3)
>>> evals_all[-3:]
array([1.53092498, 1.77243671, 2.00582508])
>>> evals_large
array([1.53092498, 1.77243671, 2.00582508])
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1\. 0\. 0.], # may vary (signs)
[-0\. -1\. 0.],
[ 0\. 0\. -1.]]
在这种情况下，我们创建了一个快速且简单的Diagonal算子。外部库PyLops提供了与Diagonal算子类似的功能，以及其他几个算子。
最后，我们考虑一个线性算子，模仿一阶导数模板的应用。在这种情况下，该算子等效于一个实非对称矩阵。再次，我们将估计的特征值和特征向量与将相同的一阶导数应用于输入信号的密集矩阵进行比较：
>>> class FirstDerivative(LinearOperator):
... def __init__(self, N, dtype='float32'):
... self.N = N
... self.shape = (self.N, self.N)
... self.dtype = np.dtype(dtype)
... def _matvec(self, x):
... y = np.zeros(self.N, self.dtype)
... y[1:-1] = (0.5*x[2:]-0.5*x[0:-2])
... return y
... def _rmatvec(self, x):
... y = np.zeros(self.N, self.dtype)
... y[0:-2] = y[0:-2] - (0.5*x[1:-1])
... y[2:] = y[2:] + (0.5*x[1:-1])
... return y
>>> N = 21
>>> D = np.diag(0.5*np.ones(N-1), k=1) - np.diag(0.5*np.ones(N-1), k=-1)
>>> D[0] = D[-1] = 0 # take away edge effects
>>> Dop = FirstDerivative(N, dtype=np.float64)
>>> evals_all, evecs_all = eig(D)
>>> evals_large, evecs_large = eigs(Dop, 4, which='LI')
>>> evals_all_imag = evals_all.imag
>>> isort_imag = np.argsort(np.abs(evals_all_imag))
>>> evals_all_imag = evals_all_imag[isort_imag]
>>> evals_large_imag = evals_large.imag
>>> isort_imag = np.argsort(np.abs(evals_large_imag))
>>> evals_large_imag = evals_large_imag[isort_imag]
>>> evals_all_imag[-4:]
array([-0.95105652, 0.95105652, -0.98768834, 0.98768834])
>>> evals_large_imag
array([0.95105652, -0.95105652, 0.98768834, -0.98768834]) # may vary
注意，这个算子的特征值都是虚数。此外，scipy.sparse.linalg.eigs函数的关键字which='LI'会产生具有最大绝对虚部（正负都有）的特征值。同样，在PyLops库中有一个更高级的一阶导数算子的实现，名为FirstDerivative算子。
参考资料
压缩稀疏图例程（scipy.sparse.csgraph）
示例：单词阶梯
单词阶梯是由刘易斯·卡罗尔发明的一种文字游戏，玩家通过逐个改变一个字母来找出单词之间的路径。例如，可以这样连接“ape”和“man”：
[{\rm ape \to apt \to ait \to bit \to big \to bag \to mag \to man}]
请注意，每一步都涉及更改单词中的一个字母。这只是从“ape”到“man”的一条可能路径，但是否是最短路径呢？如果我们希望找到两个给定单词之间的最短单词阶梯路径，稀疏图子模块可以帮助。
首先，我们需要一个有效单词的列表。许多操作系统都内置了这样的列表。例如，在 Linux 上，可以在以下位置之一找到单词列表：
/usr/share/dict
/var/lib/dict
另一个获取单词的简单来源是各种互联网上的 Scrabble 单词列表（使用您喜爱的搜索引擎进行搜索）。我们首先要创建这个列表。系统单词列表由一个每行一个单词的文件组成。以下内容应修改以使用您现有的单词列表：
>>> word_list = open('/usr/share/dict/words').readlines()
>>> word_list = map(str.strip, word_list)
我们想要查看长度为 3 的单词，所以让我们只选择正确长度的单词。我们还将消除以大写字母开头（专有名词）或包含非字母数字字符（如撇号和连字符）的单词。最后，我们会确保为后续比较转换为小写：
>>> word_list = [word for word in word_list if len(word) == 3]
>>> word_list = [word for word in word_list if word[0].islower()]
>>> word_list = [word for word in word_list if word.isalpha()]
>>> word_list = list(map(str.lower, word_list))
>>> len(word_list)
586 # may vary
现在我们有一个包含 586 个有效的三个字母单词的列表（具体数字可能根据特定列表而变化）。这些单词中的每一个将成为我们图中的一个节点，并且我们将创建连接每对仅相差一个字母的单词节点的边。
有有效的方法来做到这一点，也有低效的方法。为了尽可能高效地完成这个任务，我们将使用一些复杂的 numpy 数组操作：
>>> import numpy as np
>>> word_list = np.asarray(word_list)
>>> word_list.dtype # these are unicode characters in Python 3
dtype('<U3')
>>> word_list.sort() # sort for quick searching later
我们有一个数组，其中每个条目都是三个 Unicode 字符长。我们希望找到所有只有一个字符不同的配对。我们将从将每个单词转换为三维向量开始：
>>> word_bytes = np.ndarray((word_list.size, word_list.itemsize),
... dtype='uint8',
... buffer=word_list.data)
>>> # each unicode character is four bytes long. We only need first byte
>>> # we know that there are three characters in each word
>>> word_bytes = word_bytes[:, ::word_list.itemsize//3]
>>> word_bytes.shape
(586, 3) # may vary
现在，我们将使用汉明距离来确定哪些单词对之间存在连接。汉明距离衡量两个向量之间不同条目的比例：任何两个汉明距离等于(1/N)的单词，其中(N)是字母数，将在单词阶梯中连接起来：
>>> from scipy.spatial.distance import pdist, squareform
>>> from scipy.sparse import csr_matrix
>>> hamming_dist = pdist(word_bytes, metric='hamming')
>>> # there are three characters in each word
>>> graph = csr_matrix(squareform(hamming_dist < 1.5 / 3))
在比较距离时，我们不使用相等性，因为这对于浮点值来说可能不稳定。不等式产生了期望的结果，只要单词列表中没有两个条目完全相同。现在，我们的图已经设置好了，我们将使用最短路径搜索来找到图中任意两个单词之间的路径：
>>> i1 = word_list.searchsorted('ape')
>>> i2 = word_list.searchsorted('man')
>>> word_list[i1]
'ape'
>>> word_list[i2]
'man'
我们需要检查这些是否匹配，因为如果单词不在列表中，情况就不同。现在，我们只需要找到图中这两个索引之间的最短路径。我们将使用Dijkstra 算法，因为它允许我们仅为一个节点找到路径：
>>> from scipy.sparse.csgraph import dijkstra
>>> distances, predecessors = dijkstra(graph, indices=i1,
... return_predecessors=True)
>>> print(distances[i2])
5.0 # may vary
因此，我们看到“猿”和“人”之间的最短路径只包含五步。我们可以使用算法返回的前导来重构这条路径：
>>> path = []
>>> i = i2
>>> while i != i1:
... path.append(word_list[i])
... i = predecessors[i]
>>> path.append(word_list[i1])
>>> print(path[::-1])
['ape', 'apt', 'opt', 'oat', 'mat', 'man'] # may vary
比我们最初的例子少了三个链接：从“猿”到“人”的路径只有五步。
使用模块中的其他工具，我们可以回答其他问题。例如，有没有在单词梯子中没有链接的三个字母单词？这是关于图中连通分量的一个问题：
>>> from scipy.sparse.csgraph import connected_components
>>> N_components, component_list = connected_components(graph)
>>> print(N_components)
15 # may vary
在这个特定的三个字母单词样本中，有 15 个连通分量：即，有 15 个不同的单词集合，这些集合之间没有路径。每个集合中有多少个单词？我们可以从连通分量的列表中学到这些信息：
>>> [np.sum(component_list == i) for i in range(N_components)]
[571, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] # may vary
有一个大的连通集和 14 个较小的连通集。让我们来看看较小连通集中的单词：
>>> [list(word_list[np.nonzero(component_list == i)]) for i in range(1, N_components)]
[['aha'], # may vary
['chi'],
['ebb'],
['ems', 'emu'],
['gnu'],
['ism'],
['khz'],
['nth'],
['ova'],
['qua'],
['ugh'],
['ups'],
['urn'],
['use']]
这些都是不通过单词梯子与其他单词连接的三个字母单词。
我们可能还对哪些单词之间的分离最大感到好奇。哪两个单词需要最多的链接才能连接起来？我们可以通过计算所有最短路径的矩阵来确定这一点。请注意，按照惯例，两个非连接点之间的距离被报告为无穷大，因此在找到最大值之前我们需要将这些移除：
>>> distances, predecessors = dijkstra(graph, return_predecessors=True)
>>> max_distance = np.max(distances[~np.isinf(distances)])
>>> print(max_distance)
13.0 # may vary
因此，至少有一对单词需要 13 步才能从一个单词到另一个单词！让我们确定这些是哪些：
>>> i1, i2 = np.nonzero(distances == max_distance)
>>> list(zip(word_list[i1], word_list[i2]))
[('imp', 'ohm'), # may vary
('imp', 'ohs'),
('ohm', 'imp'),
('ohm', 'ump'),
('ohs', 'imp'),
('ohs', 'ump'),
('ump', 'ohm'),
('ump', 'ohs')]
我们看到有两对单词彼此之间的最大分离：一方面是‘imp’和‘ump’，另一方面是‘ohm’和‘ohs’。我们可以以与上述相同的方式找到连接列表：
>>> path = []
>>> i = i2[0]
>>> while i != i1[0]:
... path.append(word_list[i])
... i = predecessors[i1[0], i]
>>> path.append(word_list[i1[0]])
>>> print(path[::-1])
['imp', 'amp', 'asp', 'ass', 'ads', 'add', 'aid', 'mid', 'mod', 'moo', 'too', 'tho', 'oho', 'ohm'] # may vary
这给我们展示了我们所期望看到的路径。
单词梯子只是 scipy 稀疏矩阵快速图算法的一个潜在应用。图论在数学、数据分析和机器学习的许多领域中都有出现。稀疏图工具足够灵活，可以处理许多这些情况。
空间数据结构和算法（scipy.spatial）
这告诉我们，这个三角形有三角形＃0 作为邻居，但没有其他邻居。此外，它告诉我们，邻居 0 位于三角形的顶点 1 的对面：
Qhull 还可以对更高维度点集（例如，在 3D 中划分成四面体）执行简单化到简单形。
凸包
出现了两个新的三角形。但我们看到它们是退化的，面积为零。
我们可以可视化它：
>>> from scipy.spatial import Delaunay
>>> import numpy as np
>>> points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])
>>> tri = Delaunay(points)
>>> import matplotlib.pyplot as plt
>>> plt.triplot(points[:,0], points[:,1], tri.simplices)
>>> plt.plot(points[:,0], points[:,1], 'o')
然而，Qhull 提供了“QJ”选项，指示它随机扰动输入数据直到去除退化情况：
>>> for j, p in enumerate(points):
... plt.text(p[0]-0.03, p[1]+0.03, j, ha='right') # label the points
>>> for j, s in enumerate(tri.simplices):
... p = points[s].mean(axis=0)
... plt.text(p[0], p[1], '#%d' % j, ha='center') # label triangles
>>> plt.xlim(-0.5, 1.5); plt.ylim(-0.5, 1.5)
>>> plt.show()
Delaunay 三角化是将一组点分割成一组不重叠的三角形的过程，使得任何点都不在任何三角形的外接圆内。在实践中，这种三角化倾向于避免具有小角度的三角形。
可以使用scipy.spatial计算 Delaunay 三角化，如下所示：
>>> i = 1
>>> tri.simplices[i,:]
array([3, 1, 0], dtype=int32)
>>> points[tri.simplices[i,:]]
array([[ 1\. , 1\. ],
[ 0\. , 1.1],
[ 0\. , 0\. ]])
注意，这种退化情况不仅可能因为重复点而发生，甚至在乍一看似乎行为良好的点集中，也可能因更复杂的几何原因而发生。
>>> tri.neighbors[i]
array([-1, 0, -1], dtype=int32)
并且添加一些进一步的装饰：
>>> points[tri.simplices[i, 1]]
array([ 0\. , 1.1])
此外，还可以找到相邻的三角形：
的确，从图中我们看到这种情况。
观察到点＃4 是一个重复的点，并未出现在三角化的顶点中。这种情况已记录：
需要注意，并非所有的点都必然出现在三角化的顶点中，这是由于在形成三角化过程中的数值精度问题。
>>> points = np.array([[0, 0], [0, 1], [1, 0], [1, 1], [1, 1]])
>>> tri = Delaunay(points)
>>> np.unique(tri.simplices.ravel())
array([0, 1, 2, 3], dtype=int32)
这意味着点 4 位于三角形 0 和顶点 3 附近，但未包含在三角化中。
>>> tri.coplanar
array([[4, 0, 3]], dtype=int32)
Delaunay 三角化
scipy.spatial可以通过利用Qhull库来计算一组点的三角化、Voronoi 图和凸包。
考虑到上述的重复点：
>>> tri = Delaunay(points, qhull_options="QJ Pp")
>>> points[tri.simplices]
array([[[1, 0],
[1, 1],
[0, 0]],
[[1, 1],
[1, 1],
[1, 0]],
[[1, 1],
[0, 1],
[0, 0]],
[[0, 1],
[1, 1],
[1, 1]]])
共面点
此外，它包含了KDTree用于最近邻点查询的实现，以及在各种度量中进行距离计算的实用程序。
凸包是包含给定点集中所有点的最小凸对象。
可以通过scipy.spatial中的 Qhull 包装器计算如下：
>>> from scipy.spatial import ConvexHull
>>> rng = np.random.default_rng()
>>> points = rng.random((30, 2)) # 30 random points in 2-D
>>> hull = ConvexHull(points)
凸包表示为 N 个 1 维简单形式的集合，在二维中意味着线段。存储方案与上面讨论的 Delaunay 三角剖分中的简单形式完全相同。
我们可以说明上述结果：
>>> import matplotlib.pyplot as plt
>>> plt.plot(points[:,0], points[:,1], 'o')
>>> for simplex in hull.simplices:
... plt.plot(points[simplex,0], points[simplex,1], 'k-')
>>> plt.show()
使用scipy.spatial.convex_hull_plot_2d也可以实现同样的效果。
Voronoi 图
Voronoi 图是将空间分割为给定一组点的最近邻域的子集。
使用scipy.spatial有两种方法来处理此对象。首先，可以使用KDTree回答“这个点最接近哪些点”的问题，并以此定义区域：
>>> from scipy.spatial import KDTree
>>> points = np.array([[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2],
... [2, 0], [2, 1], [2, 2]])
>>> tree = KDTree(points)
>>> tree.query([0.1, 0.1])
(0.14142135623730953, 0)
因此，点(0.1, 0.1)属于区域0。以颜色显示：
>>> x = np.linspace(-0.5, 2.5, 31)
>>> y = np.linspace(-0.5, 2.5, 33)
>>> xx, yy = np.meshgrid(x, y)
>>> xy = np.c_[xx.ravel(), yy.ravel()]
>>> import matplotlib.pyplot as plt
>>> dx_half, dy_half = np.diff(x[:2])[0] / 2., np.diff(y[:2])[0] / 2.
>>> x_edges = np.concatenate((x - dx_half, [x[-1] + dx_half]))
>>> y_edges = np.concatenate((y - dy_half, [y[-1] + dy_half]))
>>> plt.pcolormesh(x_edges, y_edges, tree.query(xy)[1].reshape(33, 31), shading='flat')
>>> plt.plot(points[:,0], points[:,1], 'ko')
>>> plt.show()
然而，这并不会将 Voronoi 图作为几何对象给出。
用scipy.spatial中的 Qhull 包装器再次以线段和点的形式表示：
>>> from scipy.spatial import Voronoi
>>> vor = Voronoi(points)
>>> vor.vertices
array([[0.5, 0.5],
[0.5, 1.5],
[1.5, 0.5],
[1.5, 1.5]])
Voronoi 顶点表示形成 Voronoi 区域多边形边缘的点集。在这种情况下，有 9 个不同的区域：
>>> vor.regions
[[], [-1, 0], [-1, 1], [1, -1, 0], [3, -1, 2], [-1, 3], [-1, 2], [0, 1, 3, 2], [2, -1, 0], [3, -1, 1]]
负值-1再次表示无穷远点。事实上，只有一个区域 [0, 1, 3, 2] 是有界的。请注意，由于与上述 Delaunay 三角剖分中类似的数值精度问题，Voronoi 区域可能少于输入点。
区域之间的脊（二维中的线）被描述为凸包部分的简单形式集合：
>>> vor.ridge_vertices
[[-1, 0], [-1, 0], [-1, 1], [-1, 1], [0, 1], [-1, 3], [-1, 2], [2, 3], [-1, 3], [-1, 2], [1, 3], [0, 2]]
这些数字表示组成线段的 Voronoi 顶点的索引。-1再次表示无穷远点 —— 只有 12 条线中的 4 条是有界线段，其他的延伸到无穷远。
Voronoi 脊是在输入点之间绘制的线段的垂直线。每个脊对应的两个点也记录在案：
>>> vor.ridge_points
array([[0, 3],
[0, 1],
[2, 5],
[2, 1],
[1, 4],
[7, 8],
[7, 6],
[7, 4],
[8, 5],
[6, 3],
[4, 5],
[4, 3]], dtype=int32)
这些信息综合起来足以构建完整的图表。
我们可以按如下方式绘制它。首先是点和 Voronoi 顶点：
>>> plt.plot(points[:, 0], points[:, 1], 'o')
>>> plt.plot(vor.vertices[:, 0], vor.vertices[:, 1], '*')
>>> plt.xlim(-1, 3); plt.ylim(-1, 3)
绘制有限线段的方式与凸包类似，但现在我们必须防范无限边：
>>> for simplex in vor.ridge_vertices:
... simplex = np.asarray(simplex)
... if np.all(simplex >= 0):
... plt.plot(vor.vertices[simplex, 0], vor.vertices[simplex, 1], 'k-')
延伸到无限的脊需要更多的注意：
>>> center = points.mean(axis=0)
>>> for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
... simplex = np.asarray(simplex)
... if np.any(simplex < 0):
... i = simplex[simplex >= 0][0] # finite end Voronoi vertex
... t = points[pointidx[1]] - points[pointidx[0]] # tangent
... t = t / np.linalg.norm(t)
... n = np.array([-t[1], t[0]]) # normal
... midpoint = points[pointidx].mean(axis=0)
... far_point = vor.vertices[i] + np.sign(np.dot(midpoint - center, n)) * n * 100
... plt.plot([vor.vertices[i,0], far_point[0]],
... [vor.vertices[i,1], far_point[1]], 'k--')
>>> plt.show()
使用scipy.spatial.voronoi_plot_2d也可以创建此图。
Voronoi 图可以用来创建有趣的生成艺术。尝试调整mandala函数的设置，创作属于你自己的作品！
>>> import numpy as np
>>> from scipy import spatial
>>> import matplotlib.pyplot as plt
>>> def mandala(n_iter, n_points, radius):
... """Creates a mandala figure using Voronoi tessellations.
...
... Parameters
... ----------
... n_iter : int
... Number of iterations, i.e. how many times the equidistant points will
... be generated.
... n_points : int
... Number of points to draw per iteration.
... radius : scalar
... The radial expansion factor.
...
... Returns
... -------
... fig : matplotlib.Figure instance
...
... Notes
... -----
... This code is adapted from the work of Audrey Roy Greenfeld [1]_ and Carlos
... Focil-Espinosa [2]_, who created beautiful mandalas with Python code. That
... code in turn was based on Antonio Sánchez Chinchón's R code [3]_.
...
... References
... ----------
... .. [1] https://www.codemakesmehappy.com/2019/09/voronoi-mandalas.html
...
... .. [2] https://github.com/CarlosFocil/mandalapy
...
... .. [3] https://github.com/aschinchon/mandalas
...
... """
... fig = plt.figure(figsize=(10, 10))
... ax = fig.add_subplot(111)
... ax.set_axis_off()
... ax.set_aspect('equal', adjustable='box')
...
... angles = np.linspace(0, 2*np.pi * (1 - 1/n_points), num=n_points) + np.pi/2
... # Starting from a single center point, add points iteratively
... xy = np.array([[0, 0]])
... for k in range(n_iter):
... t1 = np.array([])
... t2 = np.array([])
... # Add `n_points` new points around each existing point in this iteration
... for i in range(xy.shape[0]):
... t1 = np.append(t1, xy[i, 0] + radius**k * np.cos(angles))
... t2 = np.append(t2, xy[i, 1] + radius**k * np.sin(angles))
...
... xy = np.column_stack((t1, t2))
...
... # Create the Mandala figure via a Voronoi plot
... spatial.voronoi_plot_2d(spatial.Voronoi(xy), ax=ax)
...
... return fig
>>> # Modify the following parameters in order to get different figures
>>> n_iter = 3
>>> n_points = 6
>>> radius = 4
>>> fig = mandala(n_iter, n_points, radius)
>>> plt.show()
统计（scipy.stats）
简介
在本教程中，我们讨论了许多scipy.stats的功能，但肯定不是所有功能。这里的意图是为用户提供对此包的工作知识。我们参考参考手册获取更多详细信息。
注意：本文档正在进行中。
-
离散统计分布
-
连续统计分布
-
SciPy 中的通用非均匀随机数采样
-
重采样和蒙特卡罗方法
随机变量
已经实现了两个通用分布类，用于封装连续随机变量和离散随机变量。使用这些类已经实现了超过 80 个连续随机变量（RVs）和 10 个离散随机变量。此外，用户可以轻松地添加新的例程和分布（如果您创建了一个，请贡献出来）。
所有统计函数位于子包scipy.stats中，可以使用info(stats)获取这些函数的相当完整的列表。还可以从 stats 子包的文档字符串中获取可用的随机变量列表。
在以下讨论中，我们主要关注连续 RVs。几乎所有内容也适用于离散变量，但我们在这里指出了一些差异：离散分布的特定点。
在下面的代码示例中，我们假设scipy.stats包已经导入为
>>> from scipy import stats
在某些情况下，我们假设已经导入了个别对象
>>> from scipy.stats import norm
获取帮助
首先，所有分布都伴随有帮助函数。要仅获得一些基本信息，我们打印相关的文档字符串：print(stats.norm.__doc__)。
要找到支持区间，即分布的上限和下限，请调用：
>>> print('bounds of distribution lower: %s, upper: %s' % norm.support())
bounds of distribution lower: -inf, upper: inf
我们可以用dir(norm)列出分布的所有方法和属性。事实证明，一些方法是私有的，尽管它们没有以下划线开头（它们的名称不是这样命名的），例如veccdf，只能用于内部计算（当试图使用它们时会发出警告，并且将在某个时候删除）。
要获得真正的主要方法，我们列出冻结分布的方法。（我们将在下面解释冻结分布的含义）。
>>> rv = norm()
>>> dir(rv) # reformatted
['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__',
'__format__', '__ge__', '__getattribute__', '__gt__', '__hash__',
'__init__', '__le__', '__lt__', '__module__', '__ne__', '__new__',
'__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__',
'__str__', '__subclasshook__', '__weakref__', 'a', 'args', 'b', 'cdf',
'dist', 'entropy', 'expect', 'interval', 'isf', 'kwds', 'logcdf',
'logpdf', 'logpmf', 'logsf', 'mean', 'median', 'moment', 'pdf', 'pmf',
'ppf', 'random_state', 'rvs', 'sf', 'stats', 'std', 'var']
最后，我们可以通过内省获得可用分布的列表：
>>> dist_continu = [d for d in dir(stats) if
... isinstance(getattr(stats, d), stats.rv_continuous)]
>>> dist_discrete = [d for d in dir(stats) if
... isinstance(getattr(stats, d), stats.rv_discrete)]
>>> print('number of continuous distributions: %d' % len(dist_continu))
number of continuous distributions: 108
>>> print('number of discrete distributions: %d' % len(dist_discrete))
number of discrete distributions: 20
常见方法
连续 RVs 的主要公共方法包括：
-
rvs：随机变量
-
pdf：概率密度函数
-
cdf：累积分布函数
-
sf：生存函数（1-CDF）
-
ppf：百分点函数（CDF 的逆）
-
isf：逆生存函数（SF 的逆）
-
统计：返回均值、方差、（Fisher 的）偏度或（Fisher 的）峰度
-
moment：分布的非中心矩
让我们以正态随机变量为例。
>>> norm.cdf(0)
0.5
要在多个点计算cdf，可以传递一个列表或 NumPy 数组。
>>> norm.cdf([-1., 0, 1])
array([ 0.15865525, 0.5, 0.84134475])
>>> import numpy as np
>>> norm.cdf(np.array([-1., 0, 1]))
array([ 0.15865525, 0.5, 0.84134475])
因此，基本方法，如pdf、cdf等，都是矢量化的。
其他通常有用的方法也受到支持：
>>> norm.mean(), norm.std(), norm.var()
(0.0, 1.0, 1.0)
>>> norm.stats(moments="mv")
(array(0.0), array(1.0))
要找到分布的中位数，我们可以使用百分点函数ppf，它是cdf的逆函数：
>>> norm.ppf(0.5)
0.0
要生成一系列随机变量，使用size关键字参数：
>>> norm.rvs(size=3)
array([-0.35687759, 1.34347647, -0.11710531]) # random
不要认为norm.rvs(5)会生成 5 个变量：
>>> norm.rvs(5)
5.471435163732493 # random
在这里，5而没有关键字被解释为第一个可能的关键字参数loc，它是所有连续分布采用的一对关键字参数中的第一个。这将我们带到下一小节的主题。
随机数生成
抽取随机数依赖于numpy.random包中的生成器。在上述示例中，特定的随机数流在多次运行中无法重现。要实现可重复性，可以显式地种子一个随机数生成器。在 NumPy 中，生成器是numpy.random.Generator的实例。以下是创建生成器的标准方式：
>>> from numpy.random import default_rng
>>> rng = default_rng()
并且可以通过以下方式固定种子：
>>> # do NOT copy this value
>>> rng = default_rng(301439351238479871608357552876690613766)
警告
不要使用此数或常见值，如 0. 仅使用一小组种子来实例化更大的状态空间意味着存在一些不可能达到的初始状态。如果每个人都使用这样的值，会造成一些偏差。获取种子的好方法是使用numpy.random.SeedSequence:
>>> from numpy.random import SeedSequence
>>> print(SeedSequence().entropy)
301439351238479871608357552876690613766 # random
分布中的random_state参数接受一个numpy.random.Generator类的实例或一个整数，然后用于种子化内部的Generator对象：
>>> norm.rvs(size=5, random_state=rng)
array([ 0.47143516, -1.19097569, 1.43270697, -0.3126519 , -0.72058873]) # random
欲了解更多信息，请参阅NumPy 文档。
要了解更多关于 SciPy 中实现的随机数采样器，请参阅非均匀随机数采样教程和准蒙特卡洛教程
移位和缩放
所有连续分布都接受loc和scale作为关键参数，以调整分布的位置和尺度，例如，对于标准正态分布，位置是均值，尺度是标准差。
>>> norm.stats(loc=3, scale=4, moments="mv")
(array(3.0), array(16.0))
在许多情况下，随机变量X的标准化分布是通过变换(X - loc) / scale来获得的。 默认值为loc = 0和scale = 1。
智能使用loc和scale可以帮助以许多方式修改标准分布。 为了进一步说明缩放，指数分布随机变量的cdf是
[F(x) = 1 - \exp(-\lambda x)]
通过应用上述的缩放规则，可以看到通过取scale = 1./lambda，我们得到了适当的比例。
>>> from scipy.stats import expon
>>> expon.mean(scale=3.)
3.0
注意
需要形状参数(a)的分布可能需要不仅仅是简单地应用loc和/或scale来实现所需的形式。 例如，给定长度为(R)的恒定向量，每个分量受独立的 N(0, (\sigma²))扰动影响的 2-D 向量长度分布为 rice((R/\sigma), scale= (\sigma))。 第一个参数是一个需要与(x)一起缩放的形状参数。
均匀分布也很有趣：
>>> from scipy.stats import uniform
>>> uniform.cdf([0, 1, 2, 3, 4, 5], loc=1, scale=4)
array([ 0\. , 0\. , 0.25, 0.5 , 0.75, 1\. ])
最后，请回顾前段提到的问题，即我们需要理解norm.rvs(5)的含义。 结果表明，像这样调用分布时，第一个参数即 5 被传递以设置loc参数。 让我们看看：
>>> np.mean(norm.rvs(5, size=500))
5.0098355106969992 # random
因此，要解释上一节示例的输出：norm.rvs(5)生成具有均值loc=5的单个正态分布随机变量，因为默认的size=1。
建议通过将值作为关键词而不是作为参数明确设置loc和scale参数。 当调用同一随机变量的多个方法时，可以使用冻结分布的技术来最小化重复，如下所述。
形状参数
虽然一般连续随机变量可以通过loc和scale参数进行移位和缩放，但某些分布需要额外的形状参数。 例如，具有密度的伽玛分布
[\gamma(x, a) = \frac{\lambda (\lambda x)^{a-1}}{\Gamma(a)} e^{-\lambda x};,]
需要形状参数(a)的分布。 注意，通过将scale关键字设置为(1/\lambda)可以获得(\lambda)。
让我们检查伽玛分布的形状参数的数量和名称。 （我们从上面知道应该是 1。）
>>> from scipy.stats import gamma
>>> gamma.numargs
1
>>> gamma.shapes
'a'
现在，我们将形状变量的值设为 1，以获得指数分布，以便轻松比较我们是否得到预期的结果。
>>> gamma(1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))
注意，我们也可以将形状参数指定为关键词：
>>> gamma(a=1, scale=2.).stats(moments="mv")
(array(2.0), array(4.0))
冻结分布
反复传递loc和scale关键字可能会变得非常麻烦。 使用冻结随机变量的概念来解决此类问题。
>>> rv = gamma(1, scale=2.)
通过使用rv，我们不再需要在每次方法调用中包含尺度或形状参数了。因此，分布可以通过两种方式之一使用，要么将所有分布参数传递给每个方法调用（例如之前所做的），要么冻结分布实例的参数。让我们来检查一下：
>>> rv.mean(), rv.std()
(2.0, 2.0)
这确实是我们应该得到的结果。
广播
基本方法pdf等符合通常的 numpy 广播规则。例如，我们可以计算不同概率和自由度的 t 分布上尾的临界值。
>>> stats.t.isf([0.1, 0.05, 0.01], [[10], [11]])
array([[ 1.37218364, 1.81246112, 2.76376946],
[ 1.36343032, 1.79588482, 2.71807918]])
在这里，第一行包含 10 自由度的临界值，第二行包含 11 自由度（d.o.f.）。因此，广播规则使得两次调用isf得到相同的结果：
>>> stats.t.isf([0.1, 0.05, 0.01], 10)
array([ 1.37218364, 1.81246112, 2.76376946])
>>> stats.t.isf([0.1, 0.05, 0.01], 11)
array([ 1.36343032, 1.79588482, 2.71807918])
如果概率数组，即[0.1, 0.05, 0.01]和自由度数组，即[10, 11, 12]具有相同的数组形状，则使用逐元素匹配。例如，我们可以通过调用以下方式获得 10 自由度的 10%尾部，11 自由度的 5%尾部和 12 自由度的 1%尾部。
>>> stats.t.isf([0.1, 0.05, 0.01], [10, 11, 12])
array([ 1.37218364, 1.79588482, 2.68099799])
离散分布的特定点
离散分布与连续分布大多数具有相同的基本方法。但是，pdf被概率质量函数pmf替代，没有可用的估计方法（例如拟合），而scale不是有效的关键字参数。仍然可以使用位置参数关键字loc来移动分布。
计算累积分布函数（cdf）需要额外注意。在连续分布的情况下，累积分布函数在界限（a，b）内在大多数标准情况下是严格单调递增的，因此具有唯一的反函数。然而，离散分布的累积分布函数是一个步骤函数，因此反函数即百分位点函数需要不同的定义：
ppf(q) = min{x : cdf(x) >= q, x integer}
获取更多信息，请参阅此处的文档 here。
我们可以以超几何分布为例
>>> from scipy.stats import hypergeom
>>> [M, n, N] = [20, 7, 12]
如果我们在某些整数点使用累积分布函数（cdf），然后在这些 cdf 值上评估百分位点函数（ppf），我们可以得到初始的整数值，例如
>>> x = np.arange(4) * 2
>>> x
array([0, 2, 4, 6])
>>> prb = hypergeom.cdf(x, M, n, N)
>>> prb
array([ 1.03199174e-04, 5.21155831e-02, 6.08359133e-01,
9.89783282e-01])
>>> hypergeom.ppf(prb, M, n, N)
array([ 0., 2., 4., 6.])
如果我们使用的值不是累积分布函数（cdf）步骤函数的拐点，则会得到下一个更高的整数值：
>>> hypergeom.ppf(prb + 1e-8, M, n, N)
array([ 1., 3., 5., 7.])
>>> hypergeom.ppf(prb - 1e-8, M, n, N)
array([ 0., 2., 4., 6.])
拟合分布
未冻结分布的主要附加方法与分布参数的估计有关：
-
fit：包括位置的分布参数的最大似然估计
和尺度
-
fit_loc_scale：在给定形状参数时估计位置和尺度
-
nnlf：负对数似然函数
-
expect：计算函数对 pdf 或 pmf 的期望
性能问题和注意事项
就个体方法的性能而言，速度方面有很大的差异，这些方法的结果是通过显式计算或与特定分布无关的通用算法之一获得的。
显式计算，一方面要求对于给定的分布直接指定方法，可以是通过解析公式或scipy.special或numpy.random中的特殊函数进行计算。这些通常是相对快速的计算。
另一方面，如果分布未指定任何显式计算，则使用通用方法。要定义一个分布，只需要 pdf 或 cdf 中的一个；所有其他方法可以使用数值积分和根查找导出。然而，这些间接方法可能非常慢。例如，rgh = stats.gausshyper.rvs(0.5, 2, 2, 2, size=100)以非常间接的方式创建随机变量，在我的计算机上为 100 个随机变量大约需要 19 秒，而从标准正态分布或 t 分布中生成一百万个随机变量仅需要超过一秒。
剩余问题
scipy.stats中的分布最近已经得到了修正和改进，并增加了相当多的测试套件；然而，还有一些问题存在：
-
已经对分布在一些参数范围内进行了测试；然而，在某些角落范围内，可能仍然存在一些不正确的结果。
-
fit中的最大似然估计不能使用所有分布的默认起始参数，并且用户需要提供良好的起始参数。此外，对于某些分布来说，使用最大似然估计可能本质上并不是最佳选择。
构建特定的分布
下一个示例展示了如何构建自己的分布。更多示例展示了分布的使用和一些统计测试。
制作连续分布，即子类化rv_continuous
制作连续分布是相当简单的。
>>> from scipy import stats
>>> class deterministic_gen(stats.rv_continuous):
... def _cdf(self, x):
... return np.where(x < 0, 0., 1.)
... def _stats(self):
... return 0., 0., 0., 0.
>>> deterministic = deterministic_gen(name="deterministic")
>>> deterministic.cdf(np.arange(-3, 3, 0.5))
array([ 0., 0., 0., 0., 0., 0., 1., 1., 1., 1., 1., 1.])
有趣的是，pdf现在自动计算：
>>> deterministic.pdf(np.arange(-3, 3, 0.5))
array([ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
5.83333333e+04, 4.16333634e-12, 4.16333634e-12,
4.16333634e-12, 4.16333634e-12, 4.16333634e-12])
注意在性能问题和警告中提到的性能问题。未指定的常见方法的计算可能变得非常慢，因为只调用通用方法，这些方法本质上不能使用关于分布的任何具体信息。因此，作为一个警示例子：
>>> from scipy.integrate import quad
>>> quad(deterministic.pdf, -1e-1, 1e-1)
(4.163336342344337e-13, 0.0)
但这是不正确的：该 pdf 的积分应该为 1. 让我们缩小积分区间：
>>> quad(deterministic.pdf, -1e-3, 1e-3) # warning removed
(1.000076872229173, 0.0010625571718182458)
看起来更好了。然而，问题的根源在于确定性分布类定义中未指定 pdf。
子类化rv_discrete
在接下来的内容中，我们使用stats.rv_discrete生成一个离散分布，该分布具有以整数为中心的间隔的截断正态的概率。
一般信息
从help(stats.rv_discrete)的文档字符串中，
“您可以通过将 (xk, pk) 元组传递给 rv_discrete 初始化方法（通过 values= 关键字），来构造一个任意的离散随机变量，其中 P{X=xk} = pk，该元组仅描述了具有非零概率的 X 值 (xk)。”
此外，这种方法需要满足一些进一步的要求：
-
关键字 name 是必需的。
-
分布的支持点 xk 必须是整数。
-
需要指定有效数字（小数点位数）。
实际上，如果不满足最后两个要求，则可能会引发异常或生成的数字可能不正确。
一个例子
让我们开始工作。首先：
>>> npoints = 20 # number of integer support points of the distribution minus 1
>>> npointsh = npoints // 2
>>> npointsf = float(npoints)
>>> nbound = 4 # bounds for the truncated normal
>>> normbound = (1+1/npointsf) * nbound # actual bounds of truncated normal
>>> grid = np.arange(-npointsh, npointsh+2, 1) # integer grid
>>> gridlimitsnorm = (grid-0.5) / npointsh * nbound # bin limits for the truncnorm
>>> gridlimits = grid - 0.5 # used later in the analysis
>>> grid = grid[:-1]
>>> probs = np.diff(stats.truncnorm.cdf(gridlimitsnorm, -normbound, normbound))
>>> gridint = grid
最后，我们可以派生自 rv_discrete：
>>> normdiscrete = stats.rv_discrete(values=(gridint,
... np.round(probs, decimals=7)), name='normdiscrete')
现在我们已经定义了分布，我们可以访问所有离散分布的常用方法。
>>> print('mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f' %
... normdiscrete.stats(moments='mvsk'))
mean = -0.0000, variance = 6.3302, skew = 0.0000, kurtosis = -0.0076
>>> nd_std = np.sqrt(normdiscrete.stats(moments='v'))
测试实现
让我们生成一个随机样本并比较观察频率与概率。
>>> n_sample = 500
>>> rvs = normdiscrete.rvs(size=n_sample)
>>> f, l = np.histogram(rvs, bins=gridlimits)
>>> sfreq = np.vstack([gridint, f, probs*n_sample]).T
>>> print(sfreq)
[[-1.00000000e+01 0.00000000e+00 2.95019349e-02] # random
[-9.00000000e+00 0.00000000e+00 1.32294142e-01]
[-8.00000000e+00 0.00000000e+00 5.06497902e-01]
[-7.00000000e+00 2.00000000e+00 1.65568919e+00]
[-6.00000000e+00 1.00000000e+00 4.62125309e+00]
[-5.00000000e+00 9.00000000e+00 1.10137298e+01]
[-4.00000000e+00 2.60000000e+01 2.24137683e+01]
[-3.00000000e+00 3.70000000e+01 3.89503370e+01]
[-2.00000000e+00 5.10000000e+01 5.78004747e+01]
[-1.00000000e+00 7.10000000e+01 7.32455414e+01]
[ 0.00000000e+00 7.40000000e+01 7.92618251e+01]
[ 1.00000000e+00 8.90000000e+01 7.32455414e+01]
[ 2.00000000e+00 5.50000000e+01 5.78004747e+01]
[ 3.00000000e+00 5.00000000e+01 3.89503370e+01]
[ 4.00000000e+00 1.70000000e+01 2.24137683e+01]
[ 5.00000000e+00 1.10000000e+01 1.10137298e+01]
[ 6.00000000e+00 4.00000000e+00 4.62125309e+00]
[ 7.00000000e+00 3.00000000e+00 1.65568919e+00]
[ 8.00000000e+00 0.00000000e+00 5.06497902e-01]
[ 9.00000000e+00 0.00000000e+00 1.32294142e-01]
[ 1.00000000e+01 0.00000000e+00 2.95019349e-02]]
接下来，我们可以测试我们的样本是否由我们的离散正态分布生成。这也验证了随机数是否生成正确。
卡方检验要求每个箱中至少有一定数量的观测值。我们将尾箱合并成更大的箱，以确保它们包含足够的观测值。
>>> f2 = np.hstack([f[:5].sum(), f[5:-5], f[-5:].sum()])
>>> p2 = np.hstack([probs[:5].sum(), probs[5:-5], probs[-5:].sum()])
>>> ch2, pval = stats.chisquare(f2, p2*n_sample)
>>> print('chisquare for normdiscrete: chi2 = %6.3f pvalue = %6.4f' % (ch2, pval))
chisquare for normdiscrete: chi2 = 12.466 pvalue = 0.4090 # random
在这种情况下，p 值很高，因此我们可以非常确信我们的随机样本确实是由该分布生成的。
分析一个样本
首先，我们创建一些随机变量。我们设置一个种子，以便每次运行时获得相同的结果以进行查看。例如，我们从学生 t 分布中取一个样本：
>>> x = stats.t.rvs(10, size=1000)
在这里，我们设置了 t 分布的所需形状参数，它在统计学上对应于自由度，设置 size=1000 意味着我们的样本包含 1000 个独立抽取的（伪）随机数。由于我们未指定关键字参数 loc 和 scale，它们被设置为它们的默认值零和一。
描述统计
x 是一个 numpy 数组，我们可以直接访问所有数组方法，例如，
>>> print(x.min()) # equivalent to np.min(x)
-3.78975572422 # random
>>> print(x.max()) # equivalent to np.max(x)
5.26327732981 # random
>>> print(x.mean()) # equivalent to np.mean(x)
0.0140610663985 # random
>>> print(x.var()) # equivalent to np.var(x))
1.28899386208 # random
样本属性如何与其理论对应物比较？
>>> m, v, s, k = stats.t.stats(10, moments='mvsk')
>>> n, (smin, smax), sm, sv, ss, sk = stats.describe(x)
>>> sstr = '%-14s mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'
>>> print(sstr % ('distribution:', m, v, s ,k))
distribution: mean = 0.0000, variance = 1.2500, skew = 0.0000, kurtosis = 1.0000 # random
>>> print(sstr % ('sample:', sm, sv, ss, sk))
sample: mean = 0.0141, variance = 1.2903, skew = 0.2165, kurtosis = 1.0556 # random
注意：stats.describe 使用无偏估计器计算方差，而 np.var 使用的是有偏估计器。
对于我们的样本，样本统计数据与它们的理论对应物略有不同。
T 检验和 KS 检验
我们可以使用 t 检验来测试我们的样本平均值是否在统计上显著地不同于理论期望。
>>> print('t-statistic = %6.3f pvalue = %6.4f' % stats.ttest_1samp(x, m))
t-statistic = 0.391 pvalue = 0.6955 # random
P 值为 0.7，这意味着在例如 10%的 alpha 误差下，我们无法拒绝样本均值等于零的假设，即标准 t 分布的期望。
作为练习，我们也可以直接计算我们的 t 检验，而不使用提供的函数，这应该会给我们相同的答案，确实如此：
>>> tt = (sm-m)/np.sqrt(sv/float(n)) # t-statistic for mean
>>> pval = stats.t.sf(np.abs(tt), n-1)*2 # two-sided pvalue = Prob(abs(t)>tt)
>>> print('t-statistic = %6.3f pvalue = %6.4f' % (tt, pval))
t-statistic = 0.391 pvalue = 0.6955 # random
Kolmogorov-Smirnov 检验可用于测试样本是否来自标准 t 分布的假设。
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 't', (10,)))
KS-statistic D = 0.016 pvalue = 0.9571 # random
再次，P 值足够高，我们无法拒绝随机样本确实按照 t 分布分布的假设。在实际应用中，我们不知道底层分布是什么。如果我们对样本执行 Kolmogorov-Smirnov 检验，以检验其是否符合标准正态分布，则同样不能拒绝我们的样本是由正态分布生成的假设，因为在本例中，P 值几乎为 40%。
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 'norm'))
KS-statistic D = 0.028 pvalue = 0.3918 # random
然而，标准正态分布的方差为 1，而我们的样本方差为 1.29。如果我们对样本进行标准化，并将其与正态分布进行检验，则 P 值再次足够大，我们无法拒绝样本来自正态分布的假设。
>>> d, pval = stats.kstest((x-x.mean())/x.std(), 'norm')
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % (d, pval))
KS-statistic D = 0.032 pvalue = 0.2397 # random
注意：Kolmogorov-Smirnov 检验假定我们针对具有给定参数的分布进行检验，由于在最后一种情况下，我们估计了均值和方差，这一假设被违反，并且基于此的 P 值的检验统计分布不正确。
分布的尾部
最后，我们可以检查分布的上尾部。我们可以使用百分位点函数 PPF，它是 CDF 函数的反函数，来获取临界值，或者更直接地，我们可以使用生存函数的反函数。
>>> crit01, crit05, crit10 = stats.t.ppf([1-0.01, 1-0.05, 1-0.10], 10)
>>> print('critical values from ppf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % (crit01, crit05, crit10))
critical values from ppf at 1%, 5% and 10% 2.7638 1.8125 1.3722
>>> print('critical values from isf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % tuple(stats.t.isf([0.01,0.05,0.10],10)))
critical values from isf at 1%, 5% and 10% 2.7638 1.8125 1.3722
>>> freq01 = np.sum(x>crit01) / float(n) * 100
>>> freq05 = np.sum(x>crit05) / float(n) * 100
>>> freq10 = np.sum(x>crit10) / float(n) * 100
>>> print('sample %%-frequency at 1%%, 5%% and 10%% tail %8.4f %8.4f %8.4f' % (freq01, freq05, freq10))
sample %-frequency at 1%, 5% and 10% tail 1.4000 5.8000 10.5000 # random
在所有三种情况下，我们的样本在尾部的权重比底层分布更大。我们可以简要检查一个更大的样本，看看是否能更接近匹配。在这种情况下，经验频率与理论概率非常接近，但如果我们重复几次，波动仍然相当大。
>>> freq05l = np.sum(stats.t.rvs(10, size=10000) > crit05) / 10000.0 * 100
>>> print('larger sample %%-frequency at 5%% tail %8.4f' % freq05l)
larger sample %-frequency at 5% tail 4.8000 # random
我们还可以将其与正态分布的尾部进行比较，正态分布在尾部的权重较小：
>>> print('tail prob. of normal at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' %
... tuple(stats.norm.sf([crit01, crit05, crit10])*100))
tail prob. of normal at 1%, 5% and 10% 0.2857 3.4957 8.5003
卡方检验可用于测试有限数量的箱子，观察到的频率是否显著不同于假设分布的概率。
>>> quantiles = [0.0, 0.01, 0.05, 0.1, 1-0.10, 1-0.05, 1-0.01, 1.0]
>>> crit = stats.t.ppf(quantiles, 10)
>>> crit
array([ -inf, -2.76376946, -1.81246112, -1.37218364, 1.37218364,
1.81246112, 2.76376946, inf])
>>> n_sample = x.size
>>> freqcount = np.histogram(x, bins=crit)[0]
>>> tprob = np.diff(quantiles)
>>> nprob = np.diff(stats.norm.cdf(crit))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t: chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t: chi2 = 2.30 pvalue = 0.8901 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 64.60 pvalue = 0.0000 # random
我们看到标准正态分布明显被拒绝，而标准 t 分布则无法被拒绝。由于我们样本的方差与两个标准分布都不同，我们可以再次考虑估计尺度和位置来重新进行测试。
分布的拟合方法可用于估计分布的参数，并且使用估计分布的概率重复进行测试。
>>> tdof, tloc, tscale = stats.t.fit(x)
>>> nloc, nscale = stats.norm.fit(x)
>>> tprob = np.diff(stats.t.cdf(crit, tdof, loc=tloc, scale=tscale))
>>> nprob = np.diff(stats.norm.cdf(crit, loc=nloc, scale=nscale))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t: chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t: chi2 = 1.58 pvalue = 0.9542 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 11.08 pvalue = 0.0858 # random
考虑到估计的参数，我们仍然可以拒绝我们的样本来自正态分布的假设（在 5% 的水平上），但同样地，我们无法拒绝 t 分布，其 p 值为 0.95。
正态分布的特殊测试
由于正态分布是统计学中最常见的分布，有几个额外的函数可用于测试样本是否可能来自正态分布。
首先，我们可以测试我们样本的偏度和峰度是否显著不同于正态分布：
>>> print('normal skewtest teststat = %6.3f pvalue = %6.4f' % stats.skewtest(x))
normal skewtest teststat = 2.785 pvalue = 0.0054 # random
>>> print('normal kurtosistest teststat = %6.3f pvalue = %6.4f' % stats.kurtosistest(x))
normal kurtosistest teststat = 4.757 pvalue = 0.0000 # random
这两个测试被合并到正态性测试中
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(x))
normaltest teststat = 30.379 pvalue = 0.0000 # random
在所有三个测试中，p 值非常低，我们可以拒绝我们的样本具有正态分布的假设。
由于我们样本的偏度和峰度基于中心矩，如果我们测试标准化样本，我们得到完全相同的结果：
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest((x-x.mean())/x.std()))
normaltest teststat = 30.379 pvalue = 0.0000 # random
由于正态性被强烈拒绝，我们可以检查正态测试是否对其他情况给出合理的结果：
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.t.rvs(10, size=100)))
normaltest teststat = 4.698 pvalue = 0.0955 # random
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.norm.rvs(size=1000)))
normaltest teststat = 0.613 pvalue = 0.7361 # random
当测试小样本的 t 分布观测值和大样本的正态分布观测值的正态性时，在两种情况下我们都无法拒绝样本来自正态分布的原假设。在第一种情况下，这是因为测试不足以区分小样本中的 t 分布和正态分布随机变量。
比较两个样本
在接下来的内容中，我们得到了两个样本，可以来自相同或不同的分布，我们想要测试这些样本是否具有相同的统计特性。
比较均值
测试样本的均值是否相同：
>>> rvs1 = stats.norm.rvs(loc=5, scale=10, size=500)
>>> rvs2 = stats.norm.rvs(loc=5, scale=10, size=500)
>>> stats.ttest_ind(rvs1, rvs2)
Ttest_indResult(statistic=-0.5489036175088705, pvalue=0.5831943748663959) # random
测试具有不同均值的样本：
>>> rvs3 = stats.norm.rvs(loc=8, scale=10, size=500)
>>> stats.ttest_ind(rvs1, rvs3)
Ttest_indResult(statistic=-4.533414290175026, pvalue=6.507128186389019e-06) # random
Kolmogorov-Smirnov 测试两个样本的 ks_2samp
对于从同一分布中抽取的样本，我们无法拒绝原假设，因为 p 值很高。
>>> stats.ks_2samp(rvs1, rvs2)
KstestResult(statistic=0.026, pvalue=0.9959527565364388) # random
在第二个例子中，即具有不同位置（即均值）的情况下，我们可以拒绝原假设，因为 p 值低于 1%。
>>> stats.ks_2samp(rvs1, rvs3)
KstestResult(statistic=0.114, pvalue=0.00299005061044668) # random
核密度估计
统计学中的一个常见任务是从一组数据样本中估计随机变量的概率密度函数（PDF）。这个任务被称为密度估计。最著名的工具是直方图。直方图是一种用于可视化的有用工具（主要是因为每个人都能理解它），但是它没有有效地使用可用的数据。核密度估计（KDE）是同一任务的更有效工具。gaussian_kde 估计器可用于估计单变量和多变量数据的 PDF。如果数据是单峰的话，它的效果最好。
单变量估计
我们从最少量的数据开始，以了解gaussian_kde的工作原理以及带宽选择的不同选项。从 PDF 中采样的数据显示为图的底部蓝色虚线（称为拉格图）：
>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> x1 = np.array([-7, -5, 1, 4, 5], dtype=np.float64)
>>> kde1 = stats.gaussian_kde(x1)
>>> kde2 = stats.gaussian_kde(x1, bw_method='silverman')
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20) # rug plot
>>> x_eval = np.linspace(-10, 10, num=200)
>>> ax.plot(x_eval, kde1(x_eval), 'k-', label="Scott's Rule")
>>> ax.plot(x_eval, kde2(x_eval), 'r-', label="Silverman's Rule")
>>> plt.show()
我们发现Scott's Rule和Silverman's Rule之间的差异非常小，而且在有限的数据量下，带宽选择可能有点过宽。我们可以定义自己的带宽函数来获得更少平滑的结果。
>>> def my_kde_bandwidth(obj, fac=1./5):
... """We use Scott's Rule, multiplied by a constant factor."""
... return np.power(obj.n, -1./(obj.d+4)) * fac
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(x1, np.zeros(x1.shape), 'b+', ms=20) # rug plot
>>> kde3 = stats.gaussian_kde(x1, bw_method=my_kde_bandwidth)
>>> ax.plot(x_eval, kde3(x_eval), 'g-', label="With smaller BW")
>>> plt.show()
我们发现，如果将带宽设置得非常窄，得到的概率密度函数（PDF）估计值就简单地是围绕每个数据点的高斯函数之和。
现在我们来看一个更现实的例子，并比较两种可用带宽选择规则之间的差异。这些规则已知对（接近）正态分布很有效，但即使对于非常强烈非正态的单峰分布，它们也能工作得相当好。作为一个非正态分布，我们采用自由度为 5 的学生 T 分布。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng()
x1 = rng.normal(size=200) # random data, normal distribution
xs = np.linspace(x1.min()-1, x1.max()+1, 200)
kde1 = stats.gaussian_kde(x1)
kde2 = stats.gaussian_kde(x1, bw_method='silverman')
fig = plt.figure(figsize=(8, 6))
ax1 = fig.add_subplot(211)
ax1.plot(x1, np.zeros(x1.shape), 'b+', ms=12) # rug plot
ax1.plot(xs, kde1(xs), 'k-', label="Scott's Rule")
ax1.plot(xs, kde2(xs), 'b-', label="Silverman's Rule")
ax1.plot(xs, stats.norm.pdf(xs), 'r--', label="True PDF")
ax1.set_xlabel('x')
ax1.set_ylabel('Density')
ax1.set_title("Normal (top) and Student's T$_{df=5}$ (bottom) distributions")
ax1.legend(loc=1)
x2 = stats.t.rvs(5, size=200, random_state=rng) # random data, T distribution
xs = np.linspace(x2.min() - 1, x2.max() + 1, 200)
kde3 = stats.gaussian_kde(x2)
kde4 = stats.gaussian_kde(x2, bw_method='silverman')
ax2 = fig.add_subplot(212)
ax2.plot(x2, np.zeros(x2.shape), 'b+', ms=12) # rug plot
ax2.plot(xs, kde3(xs), 'k-', label="Scott's Rule")
ax2.plot(xs, kde4(xs), 'b-', label="Silverman's Rule")
ax2.plot(xs, stats.t.pdf(xs, 5), 'r--', label="True PDF")
ax2.set_xlabel('x')
ax2.set_ylabel('Density')
plt.show()
现在我们来看一个双峰分布，其中一个特征的高斯较宽，另一个特征的高斯较窄。我们预期这将是一个更难近似的密度，因为需要精确解析每个特征所需的不同带宽。
>>> from functools import partial
>>> loc1, scale1, size1 = (-2, 1, 175)
>>> loc2, scale2, size2 = (2, 0.2, 50)
>>> x2 = np.concatenate([np.random.normal(loc=loc1, scale=scale1, size=size1),
... np.random.normal(loc=loc2, scale=scale2, size=size2)])
>>> x_eval = np.linspace(x2.min() - 1, x2.max() + 1, 500)
>>> kde = stats.gaussian_kde(x2)
>>> kde2 = stats.gaussian_kde(x2, bw_method='silverman')
>>> kde3 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.2))
>>> kde4 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.5))
>>> pdf = stats.norm.pdf
>>> bimodal_pdf = pdf(x_eval, loc=loc1, scale=scale1) * float(size1) / x2.size + \
... pdf(x_eval, loc=loc2, scale=scale2) * float(size2) / x2.size
>>> fig = plt.figure(figsize=(8, 6))
>>> ax = fig.add_subplot(111)
>>> ax.plot(x2, np.zeros(x2.shape), 'b+', ms=12)
>>> ax.plot(x_eval, kde(x_eval), 'k-', label="Scott's Rule")
>>> ax.plot(x_eval, kde2(x_eval), 'b-', label="Silverman's Rule")
>>> ax.plot(x_eval, kde3(x_eval), 'g-', label="Scott * 0.2")
>>> ax.plot(x_eval, kde4(x_eval), 'c-', label="Scott * 0.5")
>>> ax.plot(x_eval, bimodal_pdf, 'r--', label="Actual PDF")
>>> ax.set_xlim([x_eval.min(), x_eval.max()])
>>> ax.legend(loc=2)
>>> ax.set_xlabel('x')
>>> ax.set_ylabel('Density')
>>> plt.show()
正如预期的那样，由于双峰分布两个特征的不同特征尺寸，KDE 与真实 PDF 的接近程度并不如我们希望的那样。通过将默认带宽减半（Scott * 0.5），我们可以稍微改善一些，而使用比默认带宽小 5 倍的因子则不足以平滑。然而，在这种情况下，我们真正需要的是非均匀（自适应）带宽。
多变量估计
使用 gaussian_kde 我们可以进行多变量以及单变量估计。这里展示了双变量情况。首先，我们生成一些随机数据，模拟两个变量之间的相关性。
>>> def measure(n):
... """Measurement model, return two coupled measurements."""
... m1 = np.random.normal(size=n)
... m2 = np.random.normal(scale=0.5, size=n)
... return m1+m2, m1-m2
>>> m1, m2 = measure(2000)
>>> xmin = m1.min()
>>> xmax = m1.max()
>>> ymin = m2.min()
>>> ymax = m2.max()
然后我们将 KDE 应用于数据：
>>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
>>> positions = np.vstack([X.ravel(), Y.ravel()])
>>> values = np.vstack([m1, m2])
>>> kernel = stats.gaussian_kde(values)
>>> Z = np.reshape(kernel.evaluate(positions).T, X.shape)
最后，我们将估计的双变量分布作为色彩映射图绘制出来，并在顶部绘制个别数据点。
>>> fig = plt.figure(figsize=(8, 6))
>>> ax = fig.add_subplot(111)
>>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
... extent=[xmin, xmax, ymin, ymax])
>>> ax.plot(m1, m2, 'k.', markersize=2)
>>> ax.set_xlim([xmin, xmax])
>>> ax.set_ylim([ymin, ymax])
>>> plt.show()
多尺度图相关性（MGC）
使用 multiscale_graphcorr，我们可以对高维和非线性数据进行独立性测试。在开始之前，让我们导入一些有用的包：
>>> import numpy as np
>>> import matplotlib.pyplot as plt; plt.style.use('classic')
>>> from scipy.stats import multiscale_graphcorr
让我们使用自定义绘图函数来绘制数据关系：
>>> def mgc_plot(x, y, sim_name, mgc_dict=None, only_viz=False,
... only_mgc=False):
... """Plot sim and MGC-plot"""
... if not only_mgc:
... # simulation
... plt.figure(figsize=(8, 8))
... ax = plt.gca()
... ax.set_title(sim_name + " Simulation", fontsize=20)
... ax.scatter(x, y)
... ax.set_xlabel('X', fontsize=15)
... ax.set_ylabel('Y', fontsize=15)
... ax.axis('equal')
... ax.tick_params(axis="x", labelsize=15)
... ax.tick_params(axis="y", labelsize=15)
... plt.show()
... if not only_viz:
... # local correlation map
... plt.figure(figsize=(8,8))
... ax = plt.gca()
... mgc_map = mgc_dict["mgc_map"]
... # draw heatmap
... ax.set_title("Local Correlation Map", fontsize=20)
... im = ax.imshow(mgc_map, cmap='YlGnBu')
... # colorbar
... cbar = ax.figure.colorbar(im, ax=ax)
... cbar.ax.set_ylabel("", rotation=-90, va="bottom")
... ax.invert_yaxis()
... # Turn spines off and create white grid.
... for edge, spine in ax.spines.items():
... spine.set_visible(False)
... # optimal scale
... opt_scale = mgc_dict["opt_scale"]
... ax.scatter(opt_scale[0], opt_scale[1],
... marker='X', s=200, color='red')
... # other formatting
... ax.tick_params(bottom="off", left="off")
... ax.set_xlabel('#Neighbors for X', fontsize=15)
... ax.set_ylabel('#Neighbors for Y', fontsize=15)
... ax.tick_params(axis="x", labelsize=15)
... ax.tick_params(axis="y", labelsize=15)
... ax.set_xlim(0, 100)
... ax.set_ylim(0, 100)
... plt.show()
让我们先看一些线性数据：
>>> rng = np.random.default_rng()
>>> x = np.linspace(-1, 1, num=100)
>>> y = x + 0.3 * rng.random(x.size)
可以在下方绘制模拟关系：
>>> mgc_plot(x, y, "Linear", only_viz=True)
现在，我们可以看到测试统计量、p 值和 MGC 地图在下方进行了可视化。最优尺度显示为地图上的红色“x”：
>>> stat, pvalue, mgc_dict = multiscale_graphcorr(x, y)
>>> print("MGC test statistic: ", round(stat, 1))
MGC test statistic: 1.0
>>> print("P-value: ", round(pvalue, 1))
P-value: 0.0
>>> mgc_plot(x, y, "Linear", mgc_dict, only_mgc=True)
从这里可以清楚地看出，MGC 能够确定输入数据矩阵之间的关系，因为 p 值非常低，而 MGC 检验统计量相对较高。MGC 地图指示出强烈的线性关系。直观上来说，这是因为拥有更多的邻居将有助于识别 (x) 和 (y) 之间的线性关系。在这种情况下，最优尺度等同于全局尺度，用地图上的红点标记。
对于非线性数据集，同样可以执行相同操作。以下的 (x) 和 (y) 数组来自非线性模拟：
>>> unif = np.array(rng.uniform(0, 5, size=100))
>>> x = unif * np.cos(np.pi * unif)
>>> y = unif * np.sin(np.pi * unif) + 0.4 * rng.random(x.size)
可以在下方绘制模拟关系：
>>> mgc_plot(x, y, "Spiral", only_viz=True)
现在，我们可以看到测试统计量、p 值和 MGC 地图在下方进行了可视化。最优尺度显示为地图上的红色“x”：
>>> stat, pvalue, mgc_dict = multiscale_graphcorr(x, y)
>>> print("MGC test statistic: ", round(stat, 1))
MGC test statistic: 0.2 # random
>>> print("P-value: ", round(pvalue, 1))
P-value: 0.0
>>> mgc_plot(x, y, "Spiral", mgc_dict, only_mgc=True)
从这里可以清楚地看出，MGC 能够再次确定关系，因为 p 值非常低，而 MGC 检验统计量相对较高。MGC 地图指示出强烈的非线性关系。在这种情况下，最优尺度等同于局部尺度，用地图上的红点标记。
拟准蒙特卡罗
在讨论准蒙特卡洛（QMC）之前，先快速介绍一下蒙特卡洛（MC）。MC 方法，或 MC 实验，是一类广泛的计算算法，依赖于重复随机抽样来获得数值结果。其基本概念是利用随机性来解决在原则上可能是确定性的问题。它们通常用于物理和数学问题，并且在难以或不可能使用其他方法时最为有用。MC 方法主要用于三类问题：优化、数值积分和从概率分布中生成抽样。
具有特定属性的随机数生成比听起来更复杂。简单的 MC 方法旨在采样独立同分布的点。但生成多组随机点可能会产生完全不同的结果。
在上述图中的两种情况中，点是随机生成的，不了解之前绘制的点。显然，空间的某些区域未被探索——这可能在模拟中造成问题，因为特定的点集可能会触发完全不同的行为。
蒙特卡洛（MC）的一个巨大优势是它具有已知的收敛性质。让我们来看看在 5 维空间中平方和的均值：
[f(\mathbf{x}) = \left( \sum_{j=1}^{5}x_j \right)²,]
其中，(x_j \sim \mathcal{U}(0,1))。它具有已知的均值，(\mu = 5/3+5(5-1)/4)。使用 MC 抽样，我们可以数值计算这个均值，并且近似误差遵循理论速率(O(n^{-1/2}))。
尽管收敛性是确保的，实践者倾向于希望有一个更确定性的探索过程。使用普通的 MC，可以使用种子来获得可重复的过程。但是固定种子会破坏收敛性质：一个给定的种子可以对一类问题有效，但对另一类问题无效。
为了以确定性方式穿过空间，通常使用跨越所有参数维度的常规网格，也称为饱和设计。让我们考虑单位超立方体，所有边界从 0 到 1。现在，点之间的距离为 0.1，填充单位间隔所需的点数将是 10。在二维超立方体中，相同的间距需要 100 个点，而在三维中则需要 1,000 个点。随着维度数量的增加，填充空间所需的实验数量将呈指数增长。这种指数增长被称为“维数灾难”。
>>> import numpy as np
>>> disc = 10
>>> x1 = np.linspace(0, 1, disc)
>>> x2 = np.linspace(0, 1, disc)
>>> x3 = np.linspace(0, 1, disc)
>>> x1, x2, x3 = np.meshgrid(x1, x2, x3)
为了缓解这个问题，设计了 QMC 方法。它们是确定性的，对空间覆盖良好，其中一些可以继续并保持良好的特性。与 MC 方法的主要区别在于，这些点不是独立同分布的，而是知道之前的点。因此，有些方法也被称为序列。
此图展示了 2 组 256 个点。左侧的设计是普通的 MC，而右侧的设计是使用Sobol’ 方法的 QMC 设计。我们清楚地看到 QMC 版本更均匀。点样本更接近边界，且聚集和间隙较少。
评估均匀性的一种方法是使用称为差异的度量。这里*Sobol’*点的差异比粗糙的 MC 方法更好。
回到平均值的计算，QMC 方法对于误差的收敛速率也更好。它们可以实现(O(n^{-1}))的速率，对于非常平滑的函数甚至可以实现更好的速率。此图表明*Sobol’*方法具有(O(n^{-1}))的速率：
我们参考scipy.stats.qmc的文档以获取更多数学细节。
计算差异
让我们考虑两组点集。从下图可以明显看出，左侧的设计覆盖了更多的空间，而右侧的设计则较少。可以使用称为discrepancy的度量来量化。差异越低，样本越均匀。
>>> import numpy as np
>>> from scipy.stats import qmc
>>> space_1 = np.array([[1, 3], [2, 6], [3, 2], [4, 5], [5, 1], [6, 4]])
>>> space_2 = np.array([[1, 5], [2, 4], [3, 3], [4, 2], [5, 1], [6, 6]])
>>> l_bounds = [0.5, 0.5]
>>> u_bounds = [6.5, 6.5]
>>> space_1 = qmc.scale(space_1, l_bounds, u_bounds, reverse=True)
>>> space_2 = qmc.scale(space_2, l_bounds, u_bounds, reverse=True)
>>> qmc.discrepancy(space_1)
0.008142039609053464
>>> qmc.discrepancy(space_2)
0.010456854423869011
使用 QMC 引擎
实现了几个 QMC 抽样器/引擎。这里我们看看两种最常用的 QMC 方法：Sobol 和 Halton 序列。
"""Sobol' and Halton sequences."""
from scipy.stats import qmc
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng()
n_sample = 256
dim = 2
sample = {}
# Sobol'
engine = qmc.Sobol(d=dim, seed=rng)
sample["Sobol'"] = engine.random(n_sample)
# Halton
engine = qmc.Halton(d=dim, seed=rng)
sample["Halton"] = engine.random(n_sample)
fig, axs = plt.subplots(1, 2, figsize=(8, 4))
for i, kind in enumerate(sample):
axs[i].scatter(sample[kind][:, 0], sample[kind][:, 1])
axs[i].set_aspect('equal')
axs[i].set_xlabel(r'$x_1$')
axs[i].set_ylabel(r'$x_2$')
axs[i].set_title(f'{kind}—$C² = ${qmc.discrepancy(sample[kind]):.2}')
plt.tight_layout()
plt.show()
警告
QMC 方法需要特别小心，用户必须阅读文档以避免常见的陷阱。例如，Sobol’ 方法要求点数是 2 的幂次方。此外，稀疏化、燃烧或其他点选择可能会破坏序列的性质，导致点集比 MC 方法更差。
QMC 引擎是状态感知的。这意味着可以继续序列、跳过某些点或重置它。让我们从Halton获取 5 个点。然后请求第二组 5 个点：
>>> from scipy.stats import qmc
>>> engine = qmc.Halton(d=2)
>>> engine.random(5)
array([[0.22166437, 0.07980522], # random
[0.72166437, 0.93165708],
[0.47166437, 0.41313856],
[0.97166437, 0.19091633],
[0.01853937, 0.74647189]])
>>> engine.random(5)
array([[0.51853937, 0.52424967], # random
[0.26853937, 0.30202745],
[0.76853937, 0.857583 ],
[0.14353937, 0.63536078],
[0.64353937, 0.01807683]])
现在我们重置序列。请求 5 个点会得到相同的第一个 5 个点：
>>> engine.reset()
>>> engine.random(5)
array([[0.22166437, 0.07980522], # random
[0.72166437, 0.93165708],
[0.47166437, 0.41313856],
[0.97166437, 0.19091633],
[0.01853937, 0.74647189]])
现在我们推进序列以获取相同的第二组 5 个点：
>>> engine.reset()
>>> engine.fast_forward(5)
>>> engine.random(5)
array([[0.51853937, 0.52424967], # random
[0.26853937, 0.30202745],
[0.76853937, 0.857583 ],
[0.14353937, 0.63536078],
[0.64353937, 0.01807683]])
注意
默认情况下，Sobol和Halton 都是乱序的。收敛性能更好，并且防止在高维度中出现点的边缘或明显的模式。没有实际理由不使用乱序版本。
制作 QMC 引擎，即继承QMCEngine
要创建自己的QMCEngine，必须定义几种方法。以下是一个包装numpy.random.Generator的示例。
>>> import numpy as np
>>> from scipy.stats import qmc
>>> class RandomEngine(qmc.QMCEngine):
... def __init__(self, d, seed=None):
... super().__init__(d=d, seed=seed)
... self.rng = np.random.default_rng(self.rng_seed)
...
...
... def _random(self, n=1, *, workers=1):
... return self.rng.random((n, self.d))
...
...
... def reset(self):
... self.rng = np.random.default_rng(self.rng_seed)
... self.num_generated = 0
... return self
...
...
... def fast_forward(self, n):
... self.random(n)
... return self
然后像任何其他 QMC 引擎一样使用它：
>>> engine = RandomEngine(2)
>>> engine.random(5)
array([[0.22733602, 0.31675834], # random
[0.79736546, 0.67625467],
[0.39110955, 0.33281393],
[0.59830875, 0.18673419],
[0.67275604, 0.94180287]])
>>> engine.reset()
>>> engine.random(5)
array([[0.22733602, 0.31675834], # random
[0.79736546, 0.67625467],
[0.39110955, 0.33281393],
[0.59830875, 0.18673419],
[0.67275604, 0.94180287]])
使用 QMC 的指南
-
QMC 有规则！确保阅读文档，否则可能不会比 MC 有任何好处。
-
如果需要精确(2^m)个点，请使用
Sobol。 -
Halton允许采样或跳过任意数量的点。这是以比*Sobol'*更慢的收敛速率为代价的。 -
永远不要移除序列的第一个点。这将破坏其属性。
-
乱序总是更好的。
-
如果使用基于 LHS 的方法，不能添加点而不丢失 LHS 属性。（有一些方法可以做到，但这些方法尚未实现。）
