SciPy-1-12-中文文档-二-

SciPy 1.12 中文文档（二）

原文：docs.scipy.org/doc/scipy-1.12.0/index.html

压缩稀疏图例程（scipy.sparse.csgraph）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/csgraph.html

示例：单词阶梯

单词阶梯是由刘易斯·卡罗尔发明的一种文字游戏，玩家通过逐个改变一个字母来找出单词之间的路径。例如，可以这样连接“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）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/spatial.html

这告诉我们，这个三角形有三角形＃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）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/stats.html

简介

在本教程中，我们讨论了许多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}))的速率：

我们参考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 属性。（有一些方法可以做到，但这些方法尚未实现。）

多维图像处理（scipy.ndimage）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/ndimage.html

简介

图像处理和分析通常被视为对值的二维数组进行操作。然而，在一些领域中，必须分析更高维度的图像。医学成像和生物成像是这方面的典型例子。由于其固有的多维特性，numpy非常适合于这类应用。scipy.ndimage包提供了许多通用的图像处理和分析函数，设计用于处理任意维度的数组。目前这些包括：线性和非线性滤波函数，二值形态学，B 样条插值以及对象测量。 ## 所有函数共享的属性

所有函数都具有一些共同的特性。特别是，所有函数都允许使用output参数来指定输出数组。通过这个参数，你可以指定一个数组，在操作中会就地改变这个数组并存储结果。在这种情况下，不会返回结果。通常情况下，使用output参数更加高效，因为可以利用现有数组来存储结果。

返回的数组类型取决于操作的类型，但大多数情况下等于输入的类型。然而，如果使用output参数，则结果的类型将等于指定的输出参数的类型。如果没有给出输出参数，则仍然可以指定输出结果的类型。这可以通过简单地将所需的numpy类型对象分配给输出参数来实现。例如：

>>> from scipy.ndimage import correlate
>>> import numpy as np
>>> correlate(np.arange(10), [1, 2.5])
array([ 0,  2,  6,  9, 13, 16, 20, 23, 27, 30])
>>> correlate(np.arange(10), [1, 2.5], output=np.float64)
array([  0\. ,   2.5,   6\. ,   9.5,  13\. ,  16.5,  20\. ,  23.5,  27\. ,  30.5]) 
```  ## 过滤函数

本节描述的函数都执行某种类型的空间过滤输入数组的操作：输出中的元素是相应输入元素邻域值的某种函数。我们称这些元素邻域为滤波核，通常是矩形的形状，但也可以具有任意的足迹。下面描述的许多函数允许通过*footprint*参数定义核的足迹。例如，可以如下定义十字形的核：

```py
>>> footprint = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
>>> footprint
array([[0, 1, 0],
 [1, 1, 1],
 [0, 1, 0]]) 

通常，核的原点位于核形状的维度的中心，通过将核形状的维度除以 2 来计算。例如，长度为 3 的一维核的原点位于第二个元素。例如，与长度为 3 的由 1 组成的滤波器的相关性：

>>> from scipy.ndimage import correlate1d
>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> correlate1d(a, [1, 1, 1])
array([0, 0, 1, 1, 1, 0, 0]) 

有时，为了选择内核的不同起点更为方便。因此，大多数函数支持 origin 参数，该参数指定了滤波器相对于其中心的起点。例如：

>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> correlate1d(a, [1, 1, 1], origin = -1)
array([0, 1, 1, 1, 0, 0, 0]) 

这会使结果向左移动。这个特性通常不经常需要，但可能会很有用，特别是对于具有偶数大小的滤波器。一个很好的例子是计算后向和前向差分：

>>> a = [0, 0, 1, 1, 1, 0, 0]
>>> correlate1d(a, [-1, 1])               # backward difference
array([ 0,  0,  1,  0,  0, -1,  0])
>>> correlate1d(a, [-1, 1], origin = -1)  # forward difference
array([ 0,  1,  0,  0, -1,  0,  0]) 

我们也可以按以下方式计算前向差分：

>>> correlate1d(a, [0, -1, 1])
array([ 0,  1,  0,  0, -1,  0,  0]) 

然而，相比于较大的内核，使用 origin 参数更为高效。对于多维内核，origin 可以是一个数字，此时假定各轴上的起点相等，或者是一个序列，分别给出每个轴上的起点。

由于输出元素是输入元素邻域内元素的函数，所以数组的边界需要适当处理，即在边界外提供值。这是通过假设根据特定边界条件扩展数组来完成的。在下述描述的函数中，可以使用 mode 参数来选择边界条件，该参数必须是一个字符串，表示边界条件的名称。当前支持以下边界条件：

模式	描述	示例
“nearest”	使用边界处的值	[1 2 3]->[1 1 2 3 3]
“wrap”	周期性地复制数组	[1 2 3]->[3 1 2 3 1]
“reflect”	在边界处反射数组	[1 2 3]->[1 1 2 3 3]
“mirror”	在边界处镜像数组	[1 2 3]->[2 1 2 3 2]
“constant”	使用常量值，默认为 0.0	[1 2 3]->[0 1 2 3 0]

为了保持与插值例程的一致性，还支持以下同义词：

模式	描述
“grid-constant”	等同于 “constant”*
“grid-mirror”	等同于 “reflect”
“grid-wrap”	等同于 “wrap”

• “grid-constant” 和 “constant” 在滤波操作中等效，但在插值函数中具有不同的行为。为了 API 的一致性，滤波函数接受任一名称。

“constant” 模式是特殊的，因为它需要额外的参数来指定应使用的常量值。

注意，mirror 和 reflect 模式的区别仅在于边界处的样本是否在反射时重复。对于 mirror 模式，对称点正好位于最后一个样本，因此该值不重复。这种模式也被称为整样本对称，因为对称点位于最后一个样本上。类似地，reflect 常被称为半样本对称，因为对称点位于数组边界的半样本之外。

注意

实现这种边界条件的最简单方法是将数据复制到一个更大的数组中，并根据边界条件在边界处扩展数据。对于大数组和大滤波器核，这将非常消耗内存，因此下面描述的函数使用一种不需要分配大临时缓冲区的不同方法。

相关和卷积

• 函数correlate1d沿指定轴计算 1-D 相关。数组沿指定轴的行与给定的权重相关。权重参数必须是一个 1-D 数字序列。

• 函数correlate实现输入数组与给定核的多维相关。

• 函数convolve1d沿指定轴计算 1-D 卷积。数组沿指定轴的行与给定的权重进行卷积。权重参数必须是一个 1-D 数字序列。

• 函数convolve实现输入数组与给定核的多维卷积。

  注意

  卷积本质上是在镜像核之后进行相关。因此，origin参数的行为与相关的情况不同：结果向相反方向移动。

平滑滤波器

• 函数gaussian_filter1d实现 1-D 高斯滤波器。高斯滤波器的标准差通过参数sigma传递。将order = 0 设置为与高斯核的卷积。顺序为 1、2 或 3 对应于与高斯的一、二或三阶导数的卷积。高阶导数未实现。

• gaussian_filter 函数实现了多维高斯滤波器。高斯滤波器沿着每个轴的标准差通过参数sigma作为数字序列传递。如果sigma不是一个序列而是一个单一数字，则滤波器的标准差在所有方向上都相等。滤波器的顺序可以分别为每个轴指定。顺序为 0 对应于与高斯核的卷积。顺序为 1、2 或 3 对应于与高斯的一阶、二阶或三阶导数的卷积。高阶导数未实现。order参数必须是一个数字，以指定所有轴的相同顺序，或者是一个数字序列，以指定每个轴的不同顺序。下面的示例显示了在具有不同sigma值的测试数据上应用的滤波器。order参数保持为 0。

  

  注意

  多维滤波器被实现为一系列 1-D 高斯滤波器。中间数组以与输出相同的数据类型存储。因此，对于精度较低的输出类型，结果可能不精确，因为中间结果可能以不足的精度存储。可以通过指定更精确的输出类型来防止这种情况。

• uniform_filter1d 函数沿着给定轴计算给定size的 1-D 均匀滤波器。

• uniform_filter 实现了多维均匀滤波器。均匀滤波器的大小由size参数作为整数序列给出，每个轴一个。如果size不是一个序列，而是一个单一数字，则假定所有轴上的大小相等。

  注意

  多维滤波器被实现为一系列 1-D 均匀滤波器。中间数组以与输出相同的数据类型存储。因此，对于精度较低的输出类型，结果可能不精确，因为中间结果可能以不足的精度存储。可以通过指定更精确的输出类型来防止这种情况。

基于顺序统计的滤波器

• minimum_filter1d 函数沿着给定轴计算给定size的 1-D 最小值滤波器。

• maximum_filter1d 函数沿着给定轴计算给定size的 1-D 最大值滤波器。

• minimum_filter 函数用于计算多维最小值滤波器。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• maximum_filter 函数用于计算多维最大值滤波器。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• rank_filter 函数用于计算多维等级滤波器。rank 可以小于零，例如rank = -1 表示最大元素。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• percentile_filter 函数用于计算多维百分位数滤波器。percentile 可以小于零，例如percentile = -20 等同于 percentile = 80。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• median_filter 函数用于计算多维中值滤波器。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

导数

滤波器可以通过几种方法构造。函数gaussian_filter1d，描述于平滑滤波器，可用于沿指定轴使用order参数计算导数。其他导数滤波器包括 Prewitt 和 Sobel 滤波器：

• 函数prewitt计算沿给定轴的导数。

• 函数sobel计算沿给定轴的导数。

Laplace 滤波器通过所有轴上的二阶导数之和计算。因此，可以使用不同的二阶导数函数构造不同的 Laplace 滤波器。因此，我们提供了一个通用函数，该函数接受函数参数以计算给定方向上的二阶导数。

• 函数generic_laplace使用通过derivative2传递的函数计算二阶导数来计算拉普拉斯滤波器。函数derivative2应具有以下签名

  derivative2(input, axis, output, mode, cval, *extra_arguments, **extra_keywords) 
  

  应计算沿尺寸axis的二阶导数。如果output不为None，则应将其用于输出并返回None，否则应返回结果。mode，cval具有通常的意义。

  extra_arguments和extra_keywords参数可用于传递每次调用derivative2时传递给其的额外参数元组和命名参数字典。

  例如

  >>> def d2(input, axis, output, mode, cval):
  ...     return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0)
  ...
  >>> a = np.zeros((5, 5))
  >>> a[2, 2] = 1
  >>> from scipy.ndimage import generic_laplace
  >>> generic_laplace(a, d2)
  array([[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]) 
  

  为了演示extra_arguments参数的使用，我们可以执行

  >>> def d2(input, axis, output, mode, cval, weights):
  ...     return correlate1d(input, weights, axis, output, mode, cval, 0,)
  ...
  >>> a = np.zeros((5, 5))
  >>> a[2, 2] = 1
  >>> generic_laplace(a, d2, extra_arguments = ([1, -2, 1],))
  array([[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]) 
  

  或者

  >>> generic_laplace(a, d2, extra_keywords = {'weights': [1, -2, 1]})
  array([[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]) 
  

通过提供适当的二阶导数函数，以下两个函数使用generic_laplace实现：

• 函数laplace使用离散差分计算二阶导数的拉普拉斯（即与[1, -2, 1]卷积）。

• 函数gaussian_laplace使用gaussian_filter计算二阶导数的拉普拉斯滤波器。通过参数sigma作为数字序列传递高斯滤波器沿每个轴的标准差。如果sigma不是序列而是单个数字，则滤波器的标准差在所有方向上是相等的。

梯度幅度定义为所有方向梯度平方和的平方根。与通用拉普拉斯函数类似，有一个generic_gradient_magnitude函数用于计算数组的梯度幅度。

• 函数generic_gradient_magnitude使用通过derivative计算一阶导数的函数来计算梯度幅度。函数derivative应具有以下签名

  derivative(input, axis, output, mode, cval, *extra_arguments, **extra_keywords) 
  

  它应计算维度axis上的导数。如果output不是None，则应用于输出并返回None，否则应返回结果。mode，cval具有通常的意义。

  extra_arguments和extra_keywords参数可用于传递额外参数的元组和传递给derivative的命名参数字典。

  例如，sobel函数符合所需的签名

  >>> a = np.zeros((5, 5))
  >>> a[2, 2] = 1
  >>> from scipy.ndimage import sobel, generic_gradient_magnitude
  >>> generic_gradient_magnitude(a, sobel)
  array([[ 0\.        ,  0\.        ,  0\.        ,  0\.        ,  0\.        ],
   [ 0\.        ,  1.41421356,  2\.        ,  1.41421356,  0\.        ],
   [ 0\.        ,  2\.        ,  0\.        ,  2\.        ,  0\.        ],
   [ 0\.        ,  1.41421356,  2\.        ,  1.41421356,  0\.        ],
   [ 0\.        ,  0\.        ,  0\.        ,  0\.        ,  0\.        ]]) 
  

  查看generic_laplace的文档以获取使用extra_arguments和extra_keywords参数的示例。

函数sobel和prewitt函数符合所需的签名，因此可以直接与generic_gradient_magnitude一起使用。

• 函数gaussian_gradient_magnitude使用gaussian_filter来计算梯度幅度。高斯滤波器沿每个轴的标准偏差通过参数sigma作为数字序列传递。如果sigma不是序列而是一个单一数字，则滤波器的标准偏差沿所有方向均相等。

通用滤波函数

要实现滤波函数，可以使用通用函数，它们接受一个实现滤波操作的可调用对象。这些通用函数处理输入和输出数组的迭代，以及边界条件的实现等详细信息。只需提供一个实现实际滤波工作的回调函数的可调用对象。回调函数也可以用 C 语言编写，并使用PyCapsule传递（有关更多信息，请参见在 C 中扩展 scipy.ndimage）。

• generic_filter1d函数实现了一个通用的一维滤波函数，其中实际的滤波操作必须作为 Python 函数（或其他可调用对象）提供。generic_filter1d函数迭代数组的行，并在每行调用function。传递给function的参数是numpy.float64类型的一维数组。第一个包含当前行的值。根据filter_size和origin参数，在开头和结尾进行扩展。第二个数组应该就地修改，以提供行的输出值。例如，考虑沿一个维度的相关性：

  >>> a = np.arange(12).reshape(3,4)
  >>> correlate1d(a, [1, 2, 3])
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

  使用generic_filter1d函数可以实现相同的操作，如下所示：

  >>> def fnc(iline, oline):
  ...     oline[...] = iline[:-2] + 2 * iline[1:-1] + 3 * iline[2:]
  ...
  >>> from scipy.ndimage import generic_filter1d
  >>> generic_filter1d(a, fnc, 3)
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

  在这里，默认情况下，内核的起源被假定为长度为 3 的滤波器的中间。因此，在调用函数之前，每个输入行都在开始和结束时扩展了一个值。

  可以选择定义并传递额外参数给滤波器函数。extra_arguments和extra_keywords参数可用于传递一组额外参数的元组和/或传递给每次调用的命名参数的字典。例如，我们可以将我们的滤波器的参数作为一个参数传递

  >>> def fnc(iline, oline, a, b):
  ...     oline[...] = iline[:-2] + a * iline[1:-1] + b * iline[2:]
  ...
  >>> generic_filter1d(a, fnc, 3, extra_arguments = (2, 3))
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

  或

  >>> generic_filter1d(a, fnc, 3, extra_keywords = {'a':2, 'b':3})
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

• generic_filter函数实现了一个通用的滤波器函数，其中实际的过滤操作必须作为 Python 函数（或其他可调用对象）提供。generic_filter函数迭代数组并在每个元素上调用function。function的参数是一个 1-D 数组，类型为numpy.float64，包含在滤波器足迹内的当前元素周围的值。函数应返回一个可以转换为双精度数的单个值。例如，考虑一个相关性：

  >>> a = np.arange(12).reshape(3,4)
  >>> correlate(a, [[1, 0], [0, 3]])
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

  可以使用generic_filter执行相同的操作，如下所示：

  >>> def fnc(buffer):
  ...     return (buffer * np.array([1, 3])).sum()
  ...
  >>> from scipy.ndimage import generic_filter
  >>> generic_filter(a, fnc, footprint = [[1, 0], [0, 1]])
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

  在这里，指定了一个仅包含两个元素的核足迹。因此，过滤器函数接收到一个长度等于两个的缓冲区，该缓冲区与适当的权重相乘并求和得到结果。

  当调用generic_filter函数时，必须提供矩形核的大小或核的足迹。如果提供了size参数，则必须是大小序列或单个数字，其中情况下假定过滤器沿每个轴的大小相等。如果提供了footprint参数，则必须是一个数组，通过其中的非零元素定义核的形状。

  可以选择定义并传递额外参数给滤波器函数。extra_arguments和extra_keywords参数可用于传递一组额外参数的元组和/或传递给每次调用的命名参数的字典。例如，我们可以将我们的滤波器的参数作为一个参数传递

  >>> def fnc(buffer, weights):
  ...     weights = np.asarray(weights)
  ...     return (buffer * weights).sum()
  ...
  >>> generic_filter(a, fnc, footprint = [[1, 0], [0, 1]], extra_arguments = ([1, 3],))
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

  或

  >>> generic_filter(a, fnc, footprint = [[1, 0], [0, 1]], extra_keywords= {'weights': [1, 3]})
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

这些函数按照从最后一个轴开始的行或元素进行迭代，即最后一个索引变化最快。对于重要的情况，这种迭代顺序是有保证的，以便根据空间位置调整滤波器。以下是使用实现滤波器并在迭代时跟踪当前坐标的类的示例。它执行与上述generic_filter相同的滤波操作，但还打印当前坐标：

>>> a = np.arange(12).reshape(3,4)
>>>
>>> class fnc_class:
...     def __init__(self, shape):
...         # store the shape:
...         self.shape = shape
...         # initialize the coordinates:
...         self.coordinates = [0] * len(shape)
...
...     def filter(self, buffer):
...         result = (buffer * np.array([1, 3])).sum()
...         print(self.coordinates)
...         # calculate the next coordinates:
...         axes = list(range(len(self.shape)))
...         axes.reverse()
...         for jj in axes:
...             if self.coordinates[jj] < self.shape[jj] - 1:
...                 self.coordinates[jj] += 1
...                 break
...             else:
...                 self.coordinates[jj] = 0
...         return result
...
>>> fnc = fnc_class(shape = (3,4))
>>> generic_filter(a, fnc.filter, footprint = [[1, 0], [0, 1]])
[0, 0]
[0, 1]
[0, 2]
[0, 3]
[1, 0]
[1, 1]
[1, 2]
[1, 3]
[2, 0]
[2, 1]
[2, 2]
[2, 3]
array([[ 0,  3,  7, 11],
 [12, 15, 19, 23],
 [28, 31, 35, 39]]) 

对于generic_filter1d函数，相同的方法适用，只是这个函数不会迭代正在被过滤的轴。然后generic_filter1d的示例如下：

>>> a = np.arange(12).reshape(3,4)
>>>
>>> class fnc1d_class:
...     def __init__(self, shape, axis = -1):
...         # store the filter axis:
...         self.axis = axis
...         # store the shape:
...         self.shape = shape
...         # initialize the coordinates:
...         self.coordinates = [0] * len(shape)
...
...     def filter(self, iline, oline):
...         oline[...] = iline[:-2] + 2 * iline[1:-1] + 3 * iline[2:]
...         print(self.coordinates)
...         # calculate the next coordinates:
...         axes = list(range(len(self.shape)))
...         # skip the filter axis:
...         del axes[self.axis]
...         axes.reverse()
...         for jj in axes:
...             if self.coordinates[jj] < self.shape[jj] - 1:
...                 self.coordinates[jj] += 1
...                 break
...             else:
...                 self.coordinates[jj] = 0
...
>>> fnc = fnc1d_class(shape = (3,4))
>>> generic_filter1d(a, fnc.filter, 3)
[0, 0]
[1, 0]
[2, 0]
array([[ 3,  8, 14, 17],
 [27, 32, 38, 41],
 [51, 56, 62, 65]]) 

傅里叶域滤波器

本节描述的函数在傅里叶域执行滤波操作。因此，此类函数的输入数组应与逆傅里叶变换函数兼容，例如numpy.fft模块中的函数。因此，我们需要处理可能是实部或复数傅里叶变换结果的数组。在实傅里叶变换的情况下，仅存储对称复数变换的一半。此外，需要知道在实 fft 变换之前被转换的轴的长度。本节描述的函数提供了一个n参数，在实变换的情况下，必须等于变换之前的实变换轴的长度。如果此参数小于零，则假定输入数组是复数傅里叶变换的结果。axis参数可用于指示执行实变换的轴。

• fourier_shift函数将输入数组乘以给定移位操作的多维傅里叶变换。shift参数是每个维度的移位序列或所有维度的单个值。

• fourier_gaussian函数将输入数组乘以具有给定标准差sigma的高斯滤波器的多维傅里叶变换。sigma参数是每个维度的值序列或所有维度的单个值。

• fourier_uniform函数将输入数组乘以具有给定大小size的均匀滤波器的多维傅里叶变换。size参数是每个维度的值序列或所有维度的单个值。

• fourier_ellipsoid 函数将输入数组与给定尺寸size的椭球形滤波器的多维傅里叶变换相乘。size参数是每个维度的值序列或所有维度的单个值。此函数仅实现维度为 1、2 和 3 的情况。## 插值函数

本节描述基于 B 样条理论的各种插值函数。关于 B 样条的良好介绍可参见[1]，图像插值的详细算法见[5]。

样条预滤波器

使用大于 1 阶的样条进行插值需要预滤波步骤。在 section [Interpolation functions]中描述的插值函数通过调用spline_filter应用预滤波，但可以通过将prefilter关键字设置为 False 来指示不执行此操作。如果在同一数组上执行多次插值操作，则这很有用。在这种情况下，只需进行一次预滤波并使用预滤波后的数组作为插值函数的输入更有效。以下两个函数实现预滤波：

• spline_filter1d 函数沿给定轴计算 1-D 样条滤波器。可选择提供输出数组。样条的阶数必须大于 1 且小于 6。

• spline_filter 函数计算多维样条滤波器。

  注意

  多维滤波器是通过一系列 1-D 样条滤波器实现的。中间数组与输出相同数据类型存储。因此，如果要求有限精度的输出，由于中间结果可能存储不足精度，结果可能不准确。可以通过指定高精度的输出类型来避免这种情况。

插值边界处理

所有的插值函数都使用样条插值来实现输入数组的某种几何变换。这要求将输出坐标映射到输入坐标，因此，可能需要处理超出边界的输入值。对于多维滤波函数中描述的相同问题，解决方法与过滤函数相同。因此，这些函数都支持一个 mode 参数，用于确定如何处理边界，并支持一个 cval 参数，在使用“constant”模式时给出一个常数值。下面展示了所有模式的行为，包括在非整数位置的行为。请注意，不同模式下的边界处理方式不同；reflect（又称grid-mirror）和grid-wrap涉及关于距离图像样本之间一半位置的对称或重复（虚线垂直线），而mirror和wrap将图像视为其范围恰好结束于第一个和最后一个样本点，而不是过了 0.5 个样本点。

图像样本的坐标落在每个轴上从 0 到shape[i] - 1的整数采样位置范围内，其中 i 是轴的索引。下图展示了在形状为(7, 7)的图像中位置为(3.7, 3.3)处的插值。对于阶数为 n 的插值，每个轴上涉及 n + 1 个样本。填充的圆圈显示了插值中涉及到的采样位置，在红色 x 的位置插值的图中。

插值函数

• geometric_transform 函数将任意几何变换应用于输入。给定的 mapping 函数在输出中的每个点被调用，以找到对应的输入坐标。mapping 必须是一个可调用对象，接受一个长度等于输出数组秩的元组，并返回与输入数组秩相等的对应输入坐标的元组。输出形状和输出类型可以选择性地提供。如果没有给出，它们等于输入的形状和类型。

  例如：

  >>> a = np.arange(12).reshape(4,3).astype(np.float64)
  >>> def shift_func(output_coordinates):
  ...     return (output_coordinates[0] - 0.5, output_coordinates[1] - 0.5)
  ...
  >>> from scipy.ndimage import geometric_transform
  >>> geometric_transform(a, shift_func)
  array([[ 0\.    ,  0\.    ,  0\.    ],
   [ 0\.    ,  1.3625,  2.7375],
   [ 0\.    ,  4.8125,  6.1875],
   [ 0\.    ,  8.2625,  9.6375]]) 
  

  可选地，可以定义并传递额外的参数给过滤器函数。extra_arguments 和 extra_keywords 参数可用于传递额外参数的元组和/或传递给每次调用中的导数的命名参数的字典。例如，我们可以像以下示例中一样传递偏移量。

  >>> def shift_func(output_coordinates, s0, s1):
  ...     return (output_coordinates[0] - s0, output_coordinates[1] - s1)
  ...
  >>> geometric_transform(a, shift_func, extra_arguments = (0.5, 0.5))
  array([[ 0\.    ,  0\.    ,  0\.    ],
   [ 0\.    ,  1.3625,  2.7375],
   [ 0\.    ,  4.8125,  6.1875],
   [ 0\.    ,  8.2625,  9.6375]]) 
  

  或

  >>> geometric_transform(a, shift_func, extra_keywords = {'s0': 0.5, 's1': 0.5})
  array([[ 0\.    ,  0\.    ,  0\.    ],
   [ 0\.    ,  1.3625,  2.7375],
   [ 0\.    ,  4.8125,  6.1875],
   [ 0\.    ,  8.2625,  9.6375]]) 
  

  注意

  映射函数也可以用 C 语言编写，并使用scipy.LowLevelCallable进行传递。更多信息请参阅在 C 中扩展 scipy.ndimage。

• 函数map_coordinates使用给定的坐标数组应用任意坐标变换。输出的形状从坐标数组的形状派生，通过删除第一个轴。参数coordinates用于找到输出中每个点在输入中对应的坐标。coordinates沿第一个轴的值是输入数组中找到输出值的坐标。（另见 numarray coordinates函数。）由于坐标可以是非整数坐标，因此在这些坐标处的输入值由请求的阶数的样条插值确定。

  这里是一个示例，插值一个 2D 数组在(0.5, 0.5)和(1, 2)处：

  >>> a = np.arange(12).reshape(4,3).astype(np.float64)
  >>> a
  array([[  0.,   1.,   2.],
   [  3.,   4.,   5.],
   [  6.,   7.,   8.],
   [  9.,  10.,  11.]])
  >>> from scipy.ndimage import map_coordinates
  >>> map_coordinates(a, [[0.5, 2], [0.5, 1]])
  array([ 1.3625,  7.]) 
  

• 函数affine_transform将仿射变换应用于输入数组。给定的变换矩阵和偏移量用于查找输出中每个点在输入中对应的坐标。在计算的坐标处的输入值由请求的阶数的样条插值确定。变换矩阵必须是 2-D，或者也可以作为 1-D 序列或数组给出。在后一种情况下，假定矩阵是对角的。然后应用更有效的插值算法，利用问题的可分离性。输出形状和输出类型可以选择性地提供。如果未给出，则它们与输入的形状和类型相等。

• 函数shift返回输入的移位版本，使用请求的order的样条插值。

• 函数zoom返回输入的重新缩放版本，使用请求的order的样条插值。

• 函数rotate返回在由参数axes给出的两个轴定义的平面中旋转的输入数组，使用请求的order的样条插值。角度必须以度为单位给出。如果reshape为真，则输出数组的大小将适应包含旋转输入的内容。## 形态学

二进制形态学

• generate_binary_structure 函数生成一个用于二值形态学操作的二值结构元素。必须提供结构的rank。返回的结构大小在每个方向上都等于三。每个元素的值等于 1，如果从元素到中心的欧氏距离的平方小于或等于connectivity。例如，生成 2D 4 连接和 8 连接结构如下：

  >>> from scipy.ndimage import generate_binary_structure
  >>> generate_binary_structure(2, 1)
  array([[False,  True, False],
   [ True,  True,  True],
   [False,  True, False]], dtype=bool)
  >>> generate_binary_structure(2, 2)
  array([[ True,  True,  True],
   [ True,  True,  True],
   [ True,  True,  True]], dtype=bool) 
  

这是generate_binary_structure在 3D 中的视觉呈现：

大多数二值形态学函数可以用基本操作腐蚀和膨胀来表示，这些操作可以在这里看到：

• binary_erosion 函数使用给定的结构元素对任意秩的数组进行二值腐蚀。原点参数控制结构元素的放置，如过滤器函数中所述。如果没有提供结构元素，则使用generate_binary_structure生成一个连接度等于 1 的元素。border_value参数给出边界外数组的值。腐蚀重复iterations次。如果iterations小于 1，则重复腐蚀直到结果不再改变。如果给出mask数组，则仅在每次迭代时修改对应掩码元素处具有真值的元素。

• binary_dilation 函数实现了带有给定结构元素的任意秩数组的二值膨胀。原点参数控制结构元素的放置位置，如滤波函数中所述。如果没有提供结构元素，则使用 generate_binary_structure 生成连接等于一的元素。border_value 参数指定边界外数组的值。膨胀重复 iterations 次。如果 iterations 小于一，则膨胀重复直到结果不再改变。如果给定 mask 数组，则只有在每次迭代时对应掩模元素处为 true 的元素才会被修改。

这里是使用 binary_dilation 的示例，通过多次使用数据数组作为掩模，从边界膨胀空数组来找到所有接触边界的元素：

>>> struct = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
>>> a = np.array([[1,0,0,0,0], [1,1,0,1,0], [0,0,1,1,0], [0,0,0,0,0]])
>>> a
array([[1, 0, 0, 0, 0],
 [1, 1, 0, 1, 0],
 [0, 0, 1, 1, 0],
 [0, 0, 0, 0, 0]])
>>> from scipy.ndimage import binary_dilation
>>> binary_dilation(np.zeros(a.shape), struct, -1, a, border_value=1)
array([[ True, False, False, False, False],
 [ True,  True, False, False, False],
 [False, False, False, False, False],
 [False, False, False, False, False]], dtype=bool) 

binary_erosion 和 binary_dilation 函数都有一个 iterations 参数，允许腐蚀或膨胀重复多次。用给定结构 n 次腐蚀或膨胀等同于使用自身膨胀 n-1 次的结构进行腐蚀或膨胀。提供了一个函数来计算多次膨胀自身的结构：

• iterate_structure 函数通过将输入结构自身膨胀 迭代 - 1 次来返回一个结构。

  例如：

  >>> struct = generate_binary_structure(2, 1)
  >>> struct
  array([[False,  True, False],
   [ True,  True,  True],
   [False,  True, False]], dtype=bool)
  >>> from scipy.ndimage import iterate_structure
  >>> iterate_structure(struct, 2)
  array([[False, False,  True, False, False],
   [False,  True,  True,  True, False],
   [ True,  True,  True,  True,  True],
   [False,  True,  True,  True, False],
   [False, False,  True, False, False]], dtype=bool)
  
  If the origin of the original structure is equal to 0, then it is
  also equal to 0 for the iterated structure. If not, the origin
  must also be adapted if the equivalent of the *iterations*
  erosions or dilations must be achieved with the iterated
  structure. The adapted origin is simply obtained by multiplying
  with the number of iterations. For convenience, the
  :func:`iterate_structure` also returns the adapted origin if the
  *origin* parameter is not ``None``:
  
  .. code:: python
  
   >>> iterate_structure(struct, 2, -1)
   (array([[False, False,  True, False, False],
   [False,  True,  True,  True, False],
   [ True,  True,  True,  True,  True],
   [False,  True,  True,  True, False],
   [False, False,  True, False, False]], dtype=bool), [-2, -2]) 
  

其他形态学操作可以通过腐蚀和膨胀来定义。以下函数为方便起见提供了一些这些操作：

• binary_opening 函数实现了使用给定结构元素进行任意秩数组的二进制开运算。二进制开运算相当于使用相同结构元素进行二进制腐蚀，然后进行二进制膨胀。原点参数控制结构元素的放置，如滤波函数中描述的那样。如果没有提供结构元素，则使用generate_binary_structure生成具有连接性等于一的元素。iterations 参数指定进行的腐蚀次数，然后是相同次数的膨胀。

• binary_closing 函数实现了使用给定结构元素进行任意秩数组的二进制闭运算。二进制闭运算相当于使用相同结构元素进行二进制膨胀，然后进行二进制腐蚀。原点参数控制结构元素的放置，如滤波函数中描述的那样。如果没有提供结构元素，则使用generate_binary_structure生成具有连接性等于一的元素。iterations 参数指定进行的膨胀次数，然后是相同次数的腐蚀。

• binary_fill_holes 函数用于关闭二进制图像中物体的孔洞，其中结构定义了孔洞的连接性。原点参数控制结构元素的放置，如滤波函数中描述的那样。如果没有提供结构元素，则使用generate_binary_structure生成具有连接性等于一的元素。

• binary_hit_or_miss 函数实现了对具有给定结构元素的任意秩数组的二进制击中或未击中变换。击中或未击中变换通过对输入与第一个结构元素进行腐蚀，对输入的逻辑 not 与第二个结构元素进行腐蚀，然后对这两个腐蚀结果进行逻辑 and 来计算。原点参数控制结构元素的放置位置，如 Filter functions 中所述。如果 origin2 等于 None，则设为 origin1 参数的值。如果未提供第一个结构元素，则使用 generate_binary_structure 生成一个连接度等于一的结构元素。如果未提供 structure2，则设置为 structure1 的逻辑 not。### 灰度形态学

灰度形态学操作是在具有任意值的数组上操作的二进制形态学操作的等价物。下面我们描述了侵蚀、膨胀、开运算和闭运算的灰度等价物。这些操作的实现方式与 Filter functions 中描述的滤波器类似，并且我们参考此部分来描述滤波器核和足迹的处理以及数组边界的处理。灰度形态学操作可以选择接受一个 structure 参数，该参数给出结构元素的值。如果未提供此参数，则假定结构元素是一个值为零的平坦结构元素。结构的形状可以通过 footprint 参数进行可选定义。如果未提供此参数，则假定结构是矩形的，大小等于 structure 数组的尺寸，或者由 size 参数定义（如果未提供 structure）。仅当 structure 和 footprint 都未给出时，size 参数才会被使用，此时结构元素被假定为矩形和平坦的，其尺寸由 size 给出。如果提供了 size 参数，则必须是一个尺寸序列或一个单独的数字，此时滤波器在每个轴上的尺寸是相同的。如果提供了 footprint 参数，则必须是一个数组，通过其非零元素定义核的形状。

类似于二进制腐蚀和膨胀，灰度腐蚀和膨胀也有相应的操作：

• grey_erosion 函数计算多维灰度侵蚀。

• grey_dilation 函数计算多维灰度膨胀。

灰度开运算和闭运算操作可以类似于它们的二值化对应物定义：

• grey_opening 函数实现了任意秩数组的灰度开运算。灰度开运算等效于灰度腐蚀后再进行灰度膨胀。

• grey_closing 函数实现了任意秩数组的灰度闭运算。灰度闭运算等效于灰度膨胀后再进行灰度腐蚀。

• morphological_gradient 函数实现了任意秩数组的灰度形态梯度。灰度形态梯度等于灰度膨胀与灰度腐蚀的差。

• morphological_laplace 函数实现了任意秩数组的灰度形态拉普拉斯。灰度形态拉普拉斯等于灰度膨胀与灰度腐蚀的和减去两倍的输入。

• white_tophat 函数实现了任意秩数组的白顶帽滤波器。白顶帽等于输入与灰度开运算的差。

• black_tophat 函数实现了任意秩数组的黑顶帽滤波器。黑顶帽等于灰度闭运算与输入的差。

距离变换用于计算对象中每个元素到背景的最小距离。以下函数实现了三种不同距离度量的距离变换：欧几里得距离、城市街区距离和棋盘距离。

• 函数distance_transform_cdt使用 chamfer 类型算法计算输入的距离变换，通过将每个对象元素（定义为大于零的值）替换为到背景的最短距离（所有非对象元素）。结构确定所做 chamfering 的类型。如果结构等于“cityblock”，则使用generate_binary_structure生成一个距离平方等于 1 的结构。如果结构等于“chessboard”，则使用generate_binary_structure生成一个距离平方等于数组秩的结构。这些选择对应于二维中城市街区距离和棋盘距离度量的常见解释。

  除了距离变换外，还可以计算特征变换。在这种情况下，返回结果的第一个轴上最接近背景元素的索引。可以使用return_distances和return_indices标志来指示是否返回距离变换、特征变换或两者。

  可以使用distances和indices参数来提供必须是正确大小和类型（numpy.int32）的可选输出数组。用于实现此功能的算法的基本原理在[2]中描述。

• 函数distance_transform_edt通过将每个对象元素（定义为大于零的值）替换为到背景的最短欧几里得距离，计算输入的精确欧几里得距离变换。

  除了距离变换外，还可以计算特征变换。在这种情况下，返回结果的第一个轴上最接近背景元素的索引。可以使用return_distances和return_indices标志来指示是否返回距离变换、特征变换或两者。

  可选地，可以通过sampling参数给出每个轴向的采样，它应该是与输入秩相等的长度序列，或者是一个单一的数，在这种情况下，假定沿所有轴向的采样是相等的。

  可以使用distances和indices参数来提供必须是正确大小和类型（numpy.float64和numpy.int32）的可选输出数组。用于实现此功能的算法的基本原理在[3]中描述。

• 函数distance_transform_bf使用蛮力算法计算输入的距离变换，通过将每个对象元素（由大于零的值定义）替换为到背景（所有非对象元素）的最短距离。度量必须是“euclidean”、“cityblock”或“chessboard”之一。

  除了距离变换之外，还可以计算特征变换。在这种情况下，返回结果的第一个轴上最接近的背景元素的索引。return_distances和return_indices标志可以用于指示是否必须返回距离变换、特征变换或两者。

  可选地，沿每个轴的采样可以由sampling参数给出，该参数应该是与输入秩相等的长度序列，或者是一个单独的数字，其中假定沿所有轴的采样相等。这个参数仅在欧几里得距离变换的情况下使用。

  distances和indices参数可用于提供必须具有正确大小和类型(numpy.float64和numpy.int32)的可选输出数组。

  注意

  此函数使用缓慢的蛮力算法，函数distance_transform_cdt可用于更有效地计算城市块和棋盘距离变换。函数distance_transform_edt可用于更有效地计算精确的欧几里得距离变换。

分割和标记

分割是从背景中分离感兴趣的对象的过程。可能最简单的方法是强度阈值分割，可以轻松使用numpy函数完成：

>>> a = np.array([[1,2,2,1,1,0],
...               [0,2,3,1,2,0],
...               [1,1,1,3,3,2],
...               [1,1,1,1,2,1]])
>>> np.where(a > 1, 1, 0)
array([[0, 1, 1, 0, 0, 0],
 [0, 1, 1, 0, 1, 0],
 [0, 0, 0, 1, 1, 1],
 [0, 0, 0, 0, 1, 0]]) 

结果是一个二进制图像，其中仍需要识别和标记各个对象。函数label生成一个数组，其中每个对象被分配一个唯一的编号：

• 函数label生成一个数组，其中输入中的对象用整数索引标记。它返回一个元组，包括对象标签数组和找到的对象数，除非给出了output参数，在这种情况下仅返回对象数。对象的连接性由一个结构元素定义。例如，在 2D 中使用 4 连接结构元素：

  >>> a = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
  >>> s = [[0, 1, 0], [1,1,1], [0,1,0]]
  >>> from scipy.ndimage import label
  >>> label(a, s)
  (array([[0, 1, 1, 0, 0, 0],
   [0, 1, 1, 0, 2, 0],
   [0, 0, 0, 2, 2, 2],
   [0, 0, 0, 0, 2, 0]]), 2) 
  

  这两个对象之间没有连接，因为我们无法将结构元素放置在既与第一个对象又与第二个对象重叠的方式。然而，8 连接结构元素会导致只有一个对象：

  >>> a = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
  >>> s = [[1,1,1], [1,1,1], [1,1,1]]
  >>> label(a, s)[0]
  array([[0, 1, 1, 0, 0, 0],
   [0, 1, 1, 0, 1, 0],
   [0, 0, 0, 1, 1, 1],
   [0, 0, 0, 0, 1, 0]]) 
  

  如果未提供结构元素，则通过调用generate_binary_structure（参见二元形态学）生成一个结构元素，使用连接性为 1（在 2D 中是第一个示例的 4 连接结构）。输入可以是任何类型，任何不等于零的值被视为对象的一部分。如果需要在移除不需要的对象后重新标记对象索引数组，则这是有用的。只需再次将标签函数应用于索引数组。例如：

  >>> l, n = label([1, 0, 1, 0, 1])
  >>> l
  array([1, 0, 2, 0, 3])
  >>> l = np.where(l != 2, l, 0)
  >>> l
  array([1, 0, 0, 0, 3])
  >>> label(l)[0]
  array([1, 0, 0, 0, 2]) 
  

  注意

  label使用的结构元素被假定为对称的。

有大量其他分割方法，例如，可以通过导数滤波器估计边界来获得对象的边界。其中一种方法是分水岭分割。函数watershed_ift生成一个数组，其中每个对象被分配一个唯一的标签，从定位对象边界的数组生成，例如通过梯度幅度滤波器。它使用一个包含对象初始标记的数组：

• 函数watershed_ift应用从标记算法开始的分水岭算法，使用图像森林变换，如[4]中所述。

• 该函数的输入是应用变换的数组，以及一个由唯一标签指定对象的标记数组，其中任何非零值均为标记。例如：

  >>> input = np.array([[0, 0, 0, 0, 0, 0, 0],
  ...                   [0, 1, 1, 1, 1, 1, 0],
  ...                   [0, 1, 0, 0, 0, 1, 0],
  ...                   [0, 1, 0, 0, 0, 1, 0],
  ...                   [0, 1, 0, 0, 0, 1, 0],
  ...                   [0, 1, 1, 1, 1, 1, 0],
  ...                   [0, 0, 0, 0, 0, 0, 0]], np.uint8)
  >>> markers = np.array([[1, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 2, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0]], np.int8)
  >>> from scipy.ndimage import watershed_ift
  >>> watershed_ift(input, markers)
  array([[1, 1, 1, 1, 1, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 2, 2, 2, 2, 2, 1],
   [1, 2, 2, 2, 2, 2, 1],
   [1, 2, 2, 2, 2, 2, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 1, 1, 1, 1, 1]], dtype=int8) 
  

  这里，使用了两个标记来指定一个对象（marker = 2）和背景（marker = 1）。这些标记的处理顺序是任意的：将背景的标记移动到数组的右下角会产生不同的结果：

  >>> markers = np.array([[0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 2, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 1]], np.int8)
  >>> watershed_ift(input, markers)
  array([[1, 1, 1, 1, 1, 1, 1],
   [1, 1, 1, 1, 1, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 1, 1, 1, 1, 1],
   [1, 1, 1, 1, 1, 1, 1]], dtype=int8) 
  

  结果是对象（marker = 2）较小，因为第二个标记被先处理了。如果第一个标记本应指定背景对象，则可能不是期望的效果。因此，watershed_ift明确地将具有负值的标记作为背景标记并在常规标记之后处理它们。例如，将第一个标记替换为负标记会产生与第一个示例类似的结果：

  >>> markers = np.array([[0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 2, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, -1]], np.int8)
  >>> watershed_ift(input, markers)
  array([[-1, -1, -1, -1, -1, -1, -1],
   [-1, -1,  2,  2,  2, -1, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1, -1,  2,  2,  2, -1, -1],
   [-1, -1, -1, -1, -1, -1, -1]], dtype=int8) 
  

  对象的连通性由结构元素定义。如果未提供结构元素，则通过调用 generate_binary_structure 函数生成一个结构元素（参见二进制形态学），使用连通性为一（在 2D 中为 4 连通结构）。例如，使用 8 连通结构来处理上一个示例将得到一个不同的对象：

  >>> watershed_ift(input, markers,
  ...               structure = [[1,1,1], [1,1,1], [1,1,1]])
  array([[-1, -1, -1, -1, -1, -1, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1, -1, -1, -1, -1, -1, -1]], dtype=int8) 
  

  注意

  watershed_ift 的实现将输入数据类型限制为 numpy.uint8 和 numpy.uint16。

对象测量

给定一个带标签对象的数组，可以测量各个对象的属性。find_objects 函数可用于生成一个切片列表，每个对象对应一个完全包含该对象的最小子数组：

• find_objects 函数查找标记数组中的所有对象，并返回一个切片列表，对应于包含对象的数组中最小的区域。

  例如：

  >>> a = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
  >>> l, n = label(a)
  >>> from scipy.ndimage import find_objects
  >>> f = find_objects(l)
  >>> a[f[0]]
  array([[1, 1],
   [1, 1]])
  >>> a[f[1]]
  array([[0, 1, 0],
   [1, 1, 1],
   [0, 1, 0]]) 
  

  函数 find_objects 返回所有对象的切片，除非 max_label 参数大于零，此时仅返回前 max_label 个对象。如果 label 数组中缺少索引，则返回 None 而不是切片。例如：

  >>> from scipy.ndimage import find_objects
  >>> find_objects([1, 0, 3, 4], max_label = 3)
  [(slice(0, 1, None),), None, (slice(2, 3, None),)] 
  

由 find_objects 生成的切片列表有助于找到数组中对象的位置和尺寸，也可用于对各个对象进行测量。例如，我们想要找到图像中一个对象强度的总和：

>>> image = np.arange(4 * 6).reshape(4, 6)
>>> mask = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
>>> labels = label(mask)[0]
>>> slices = find_objects(labels) 

然后我们可以计算第二个对象中元素的总和：

>>> np.where(labels[slices[1]] == 2, image[slices[1]], 0).sum()
80 

然而，这并不特别高效，对于其他类型的测量可能也更加复杂。因此，定义了一些测量函数，这些函数接受对象标签数组和要测量的对象的索引。例如，可以通过以下方式计算强度的总和：

>>> from scipy.ndimage import sum as ndi_sum
>>> ndi_sum(image, labels, 2)
80 

对于大数组和小对象，调用数组切片后再调用测量函数效率更高：

>>> ndi_sum(image[slices[1]], labels[slices[1]], 2)
80 

或者，我们可以通过单个函数调用对多个标签进行测量，返回一个结果列表。例如，在我们的示例中，要测量背景和第二个对象的值之和，我们提供一个标签列表：

>>> ndi_sum(image, labels, [0, 2])
array([178.0, 80.0]) 

下面描述的测量函数都支持index参数，以指示应测量哪些对象。index的默认值为None。这表示所有标签大于零的元素都应视为单个对象并进行测量。因此，在这种情况下，labels数组被视为由大于零的元素定义的掩码。如果index是一个数字或数字序列，则给出被测量对象的标签。如果index是一个序列，则返回结果的列表。如果返回多个结果的函数在index为单个数字时返回其结果作为元组，或在index为序列时返回其结果作为列表的元组。

• sum 函数计算由index给定标签的对象元素的总和，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• mean 函数计算由index给定标签的对象元素的平均值，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• variance 函数计算由index给定标签的对象元素的方差，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• standard_deviation 函数计算由index给定标签的对象元素的标准差，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• minimum 函数计算由index给定标签的对象元素的最小值，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• maximum 函数计算由index给定标签的对象元素的最大值，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• minimum_position 函数计算由index给定标签的对象元素的最小位置，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• maximum_position 函数计算由index给定标签的对象元素的最大位置，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• extrema 函数计算由index给定标签的对象元素的最小值、最大值及其位置，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。结果是一个元组，包含最小值、最大值、最小位置和最大位置。其结果与上述minimum、maximum、minimum_position 和 maximum_position 函数的结果组成的元组相同。

• center_of_mass 函数计算由index给定标签的对象的质心，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• histogram 函数使用 index 指定的标签的对象计算直方图，并使用 labels 数组作为对象标签。如果 index 为 None，则所有具有非零标签值的元素将被视为单个对象。如果 label 为 None，则计算时使用 input 的所有元素。直方图由其最小值 (min)、最大值 (max) 和箱数 (bins) 定义。它们作为 numpy.int32 类型的 1-D 数组返回。## 在 C 中扩展 scipy.ndimage

scipy.ndimage 中的几个函数接受回调参数。这可以是一个 Python 函数或一个包含指向 C 函数指针的 scipy.LowLevelCallable。使用 C 函数通常更有效率，因为它避免了在数组的许多元素上调用 Python 函数的开销。要使用 C 函数，您必须编写一个 C 扩展，其中包含回调函数和一个返回 scipy.LowLevelCallable 的 Python 函数，该函数包含指向回调函数的指针。

支持回调的一个函数示例是 geometric_transform，它接受定义从所有输出坐标到输入数组中对应坐标的映射的回调函数。考虑以下 Python 示例，它使用 geometric_transform 实现了一个移位函数。

from scipy import ndimage

def transform(output_coordinates, shift):
    input_coordinates = output_coordinates[0] - shift, output_coordinates[1] - shift
    return input_coordinates

im = np.arange(12).reshape(4, 3).astype(np.float64)
shift = 0.5
print(ndimage.geometric_transform(im, transform, extra_arguments=(shift,))) 

我们还可以使用以下 C 代码来实现回调函数：

/* example.c */

#include  <Python.h>
#include  <numpy/npy_common.h>

static  int
_transform(npy_intp  *output_coordinates,  double  *input_coordinates,
  int  output_rank,  int  input_rank,  void  *user_data)
{
  npy_intp  i;
  double  shift  =  *(double  *)user_data;

  for  (i  =  0;  i  <  input_rank;  i++)  {
  input_coordinates[i]  =  output_coordinates[i]  -  shift;
  }
  return  1;
}

static  char  *transform_signature  =  "int (npy_intp *, double *, int, int, void *)";

static  PyObject  *
py_get_transform(PyObject  *obj,  PyObject  *args)
{
  if  (!PyArg_ParseTuple(args,  ""))  return  NULL;
  return  PyCapsule_New(_transform,  transform_signature,  NULL);
}

static  PyMethodDef  ExampleMethods[]  =  {
  {"get_transform",  (PyCFunction)py_get_transform,  METH_VARARGS,  ""},
  {NULL,  NULL,  0,  NULL}
};

/* Initialize the module */
static  struct  PyModuleDef  example  =  {
  PyModuleDef_HEAD_INIT,
  "example",
  NULL,
  -1,
  ExampleMethods,
  NULL,
  NULL,
  NULL,
  NULL
};

PyMODINIT_FUNC
PyInit_example(void)
{
  return  PyModule_Create(&example);
} 

关于编写 Python 扩展模块的更多信息可以在这里找到。如果 C 代码在文件 example.c 中，则可以在将其添加到 meson.build（查看 meson.build 文件中的示例）之后进行编译并遵循其中的步骤。完成后，运行脚本：

import ctypes
import numpy as np
from scipy import ndimage, LowLevelCallable

from example import get_transform

shift = 0.5

user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)
callback = LowLevelCallable(get_transform(), ptr)
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

产生与原始的 Python 脚本相同的结果。

在 C 版本中，_transform 是回调函数，参数 output_coordinates 和 input_coordinates 的作用与 Python 版本相同，而 output_rank 和 input_rank 则提供了 len(output_coordinates) 和 len(input_coordinates) 的等效值。变量 shift 通过 user_data 而不是 extra_arguments 传递。最后，C 回调函数返回一个整数状态，成功时为 1，否则为 0。

函数 py_transform 将回调函数封装在 PyCapsule 中。主要步骤如下：

• 初始化一个 PyCapsule。第一个参数是指向回调函数的指针。

• 第二个参数是函数签名，必须与 ndimage 预期的完全匹配。

• 上述示例中，我们使用了 scipy.LowLevelCallable 来指定由 ctypes 生成的 user_data。

  另一种方法是在胶囊上下文中提供数据，可以通过 PyCapsule_SetContext 设置，且在 scipy.LowLevelCallable 中省略指定 user_data。然而，在此方法中，我们需要处理数据的分配/释放 —— 在销毁胶囊后释放数据可以通过在 PyCapsule_New 的第三个参数中指定非空回调函数来完成。

用于 ndimage 的 C 回调函数都遵循此模式。下一节列出了接受 C 回调函数并给出函数原型的 ndimage 函数。

另见

支持低级回调参数的函数有：

generic_filter，generic_filter1d，geometric_transform

下面，我们展示了使用 Numba、Cython、ctypes 或 cffi 写代码的替代方式，而不是在 C 中编写包装器代码。

Numba

Numba 提供了在 Python 中轻松编写低级函数的方式。我们可以使用 Numba 来编写上述代码：

# example.py
import numpy as np
import ctypes
from scipy import ndimage, LowLevelCallable
from numba import cfunc, types, carray

@cfunc(types.intc(types.CPointer(types.intp),
                  types.CPointer(types.double),
                  types.intc,
                  types.intc,
                  types.voidptr))
def transform(output_coordinates_ptr, input_coordinates_ptr,
              output_rank, input_rank, user_data):
    input_coordinates = carray(input_coordinates_ptr, (input_rank,))
    output_coordinates = carray(output_coordinates_ptr, (output_rank,))
    shift = carray(user_data, (1,), types.double)[0]

    for i in range(input_rank):
        input_coordinates[i] = output_coordinates[i] - shift

    return 1

shift = 0.5

# Then call the function
user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)
callback = LowLevelCallable(transform.ctypes, ptr)

im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

Cython

从功能上看，与上述代码相同的代码可以在 Cython 中写成，而且更少繁文缛节，如下所示：

# example.pyx

from numpy cimport npy_intp as intp

cdef api int transform(intp *output_coordinates, double *input_coordinates,
                       int output_rank, int input_rank, void *user_data):
    cdef intp i
    cdef double shift = (<double *>user_data)[0]

    for i in range(input_rank):
        input_coordinates[i] = output_coordinates[i] - shift
    return 1 

# script.py

import ctypes
import numpy as np
from scipy import ndimage, LowLevelCallable

import example

shift = 0.5

user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)
callback = LowLevelCallable.from_cython(example, "transform", ptr)
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

cffi

使用 cffi，您可以与驻留在共享库（DLL）中的 C 函数进行接口。首先，我们需要编写共享库，在 Linux/OSX 上如下：

/*
 example.c
 Needs to be compiled with "gcc -std=c99 -shared -fPIC -o example.so example.c"
 or similar
 */

#include  <stdint.h>

int
_transform(intptr_t  *output_coordinates,  double  *input_coordinates,
  int  output_rank,  int  input_rank,  void  *user_data)
{
  int  i;
  double  shift  =  *(double  *)user_data;

  for  (i  =  0;  i  <  input_rank;  i++)  {
  input_coordinates[i]  =  output_coordinates[i]  -  shift;
  }
  return  1;
} 

调用库的 Python 代码如下：

import os
import numpy as np
from scipy import ndimage, LowLevelCallable
import cffi

# Construct the FFI object, and copypaste the function declaration
ffi = cffi.FFI()
ffi.cdef("""
int _transform(intptr_t *output_coordinates, double *input_coordinates,
 int output_rank, int input_rank, void *user_data);
""")

# Open library
lib = ffi.dlopen(os.path.abspath("example.so"))

# Do the function call
user_data = ffi.new('double *', 0.5)
callback = LowLevelCallable(lib._transform, user_data)
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

您可以在 cffi 文档中找到更多信息。

ctypes

使用 ctypes，与上述 cffi 的情况一样，C 代码和 so/DLL 的编译是一样的。Python 代码则有所不同：

# script.py

import os
import ctypes
import numpy as np
from scipy import ndimage, LowLevelCallable

lib = ctypes.CDLL(os.path.abspath('example.so'))

shift = 0.5

user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)

# Ctypes has no built-in intptr type, so override the signature
# instead of trying to get it via ctypes
callback = LowLevelCallable(lib._transform, ptr,
    "int _transform(intptr_t *, double *, int, int, void *)")

# Perform the call
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

你可以在 ctypes 文档中找到更多信息。

参考资料

文件 IO（scipy.io）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/io.html

另请参见

NumPy IO routines

MATLAB 文件

loadmat(file_name[, mdict, appendmat])	Load MATLAB file.
savemat(file_name, mdict[, appendmat, ...])	保存字典的名称和数组到 MATLAB 风格的.mat 文件中。
whosmat(file_name[, appendmat])	列出 MATLAB 文件中的变量。

基本功能

我们将从导入scipy.io开始，并为方便起见称其为sio：

>>> import scipy.io as sio 

如果您正在使用 IPython，请尝试在sio上进行制表符完成。在众多选项中，您会找到：

sio.loadmat
sio.savemat
sio.whosmat 

这些是您在处理 MATLAB 文件时最可能使用的高级功能。您还会发现：

sio.matlab 

这是导入loadmat、savemat和whosmat的包。在sio.matlab中，您会找到mio模块。该模块包含loadmat和savemat使用的机制。偶尔您可能会发现自己重新使用此机制。

我该如何开始？

您可能有一个.mat文件，想要将其读入 SciPy。或者，您想要从 SciPy / NumPy 传递一些变量到 MATLAB。

为了避免使用 MATLAB 许可证，让我们从Octave开始。Octave 具有与 MATLAB 兼容的保存和加载功能。在命令行上启动 Octave（对我来说是octave）：

octave:1> a = 1:12
a =

   1   2   3   4   5   6   7   8   9  10  11  12

octave:2> a = reshape(a, [1 3 4])
a =

ans(:,:,1) =

   1   2   3

ans(:,:,2) =

   4   5   6

ans(:,:,3) =

   7   8   9

ans(:,:,4) =

   10   11   12

octave:3> save -6 octave_a.mat a % MATLAB 6 compatible
octave:4> ls octave_a.mat
octave_a.mat 

现在，到 Python：

>>> mat_contents = sio.loadmat('octave_a.mat')
>>> mat_contents
{'a': array([[[  1.,   4.,   7.,  10.],
 [  2.,   5.,   8.,  11.],
 [  3.,   6.,   9.,  12.]]]),
 '__version__': '1.0',
 '__header__': 'MATLAB 5.0 MAT-file, written by
 Octave 3.6.3, 2013-02-17 21:02:11 UTC',
 '__globals__': []}
>>> oct_a = mat_contents['a']
>>> oct_a
array([[[  1.,   4.,   7.,  10.],
 [  2.,   5.,   8.,  11.],
 [  3.,   6.,   9.,  12.]]])
>>> oct_a.shape
(1, 3, 4) 

现在让我们试着换个角度：

>>> import numpy as np
>>> vect = np.arange(10)
>>> vect.shape
(10,)
>>> sio.savemat('np_vector.mat', {'vect':vect}) 

然后回到 Octave：

octave:8> load np_vector.mat
octave:9> vect
vect =

  0  1  2  3  4  5  6  7  8  9

octave:10> size(vect)
ans =

    1   10 

如果要检查 MATLAB 文件的内容而不将数据读入内存，请使用whosmat命令：

>>> sio.whosmat('octave_a.mat')
[('a', (1, 3, 4), 'double')] 

whosmat返回一个元组列表，每个文件中的数组（或其他对象）都有一个。每个元组包含数组的名称、形状和数据类型。

MATLAB 结构

MATLAB 结构有点像 Python 字典，但字段名称必须是字符串。任何 MATLAB 对象都可以是字段的值。与 MATLAB 中的所有对象一样，结构实际上是结构数组，其中单个结构是形状为（1，1）的数组。

octave:11> my_struct = struct('field1', 1, 'field2', 2)
my_struct =
{
  field1 =  1
  field2 =  2
}

octave:12> save -6 octave_struct.mat my_struct 

我们可以在 Python 中加载它：

>>> mat_contents = sio.loadmat('octave_struct.mat')
>>> mat_contents
{'my_struct': array([[([[1.0]], [[2.0]])]],
 dtype=[('field1', 'O'), ('field2', 'O')]), '__version__': '1.0', '__header__': 'MATLAB 5.0 MAT-file, written by Octave 3.6.3, 2013-02-17 21:23:14 UTC', '__globals__': []}
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct.shape
(1, 1)
>>> val = oct_struct[0,0]
>>> val
([[1.0]], [[2.0]])
>>> val['field1']
array([[ 1.]])
>>> val['field2']
array([[ 2.]])
>>> val.dtype
dtype([('field1', 'O'), ('field2', 'O')]) 

在 SciPy 版本从 0.12.0 开始，MATLAB 结构返回为 NumPy 结构化数组，其字段命名为结构字段。您可以在上面的dtype输出中看到字段名称。还要注意：

>>> val = oct_struct[0,0] 

和：

octave:13> size(my_struct)
ans =

   1   1 

因此，在 MATLAB 中，结构数组必须至少是 2 维的，并且我们在读入 SciPy 时复制了这一点。如果您希望将所有长度为 1 的维度挤出，请尝试这样做：

>>> mat_contents = sio.loadmat('octave_struct.mat', squeeze_me=True)
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct.shape
() 

有时，将 MATLAB 结构加载为 Python 对象而不是 NumPy 结构化数组更方便 - 这可以使 Python 中的访问语法与 MATLAB 中的语法更加相似。为此，请使用struct_as_record=False参数设置为loadmat。

>>> mat_contents = sio.loadmat('octave_struct.mat', struct_as_record=False)
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct[0,0].field1
array([[ 1.]]) 

struct_as_record=False 与 squeeze_me 配合使用效果很好：

>>> mat_contents = sio.loadmat('octave_struct.mat', struct_as_record=False, squeeze_me=True)
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct.shape # but no - it's a scalar
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'mat_struct' object has no attribute 'shape'
>>> type(oct_struct)
<class 'scipy.io.matlab.mio5_params.mat_struct'>
>>> oct_struct.field1
1.0 

可以以多种方式保存结构数组。一种简单的方法是使用字典：

>>> a_dict = {'field1': 0.5, 'field2': 'a string'}
>>> sio.savemat('saved_struct.mat', {'a_dict': a_dict}) 

被加载为：

octave:21> load saved_struct
octave:22> a_dict
a_dict =

  scalar structure containing the fields:

    field2 = a string
    field1 =  0.50000 

您还可以像这样将结构体再次保存回 MATLAB（或者在我们的情况下是 Octave）：

>>> dt = [('f1', 'f8'), ('f2', 'S10')]
>>> arr = np.zeros((2,), dtype=dt)
>>> arr
array([(0.0, ''), (0.0, '')],
 dtype=[('f1', '<f8'), ('f2', 'S10')])
>>> arr[0]['f1'] = 0.5
>>> arr[0]['f2'] = 'python'
>>> arr[1]['f1'] = 99
>>> arr[1]['f2'] = 'not perl'
>>> sio.savemat('np_struct_arr.mat', {'arr': arr}) 

MATLAB 单元数组

MATLAB 中的单元数组与 Python 列表相似，数组中的元素可以包含任何类型的 MATLAB 对象。事实上，它们最类似于 NumPy 对象数组，这就是我们如何将它们加载到 NumPy 中的方式。

octave:14> my_cells = {1, [2, 3]}
my_cells =
{
  [1,1] =  1
  [1,2] =

     2   3

}

octave:15> save -6 octave_cells.mat my_cells 

回到 Python：

>>> mat_contents = sio.loadmat('octave_cells.mat')
>>> oct_cells = mat_contents['my_cells']
>>> print(oct_cells.dtype)
object
>>> val = oct_cells[0,0]
>>> val
array([[ 1.]])
>>> print(val.dtype)
float64 

保存到 MATLAB 单元数组只需创建一个 NumPy 对象数组：

>>> obj_arr = np.zeros((2,), dtype=np.object)
>>> obj_arr[0] = 1
>>> obj_arr[1] = 'a string'
>>> obj_arr
array([1, 'a string'], dtype=object)
>>> sio.savemat('np_cells.mat', {'obj_arr':obj_arr}) 

octave:16> load np_cells.mat
octave:17> obj_arr
obj_arr =
{
  [1,1] = 1
  [2,1] = a string
} 

IDL 文件

readsav(文件名[, idict, python_dict, ...])	读取 IDL 的 .sav 文件。

Matrix Market 文件

mminfo(源)	从类似于 Matrix Market 文件的 '源' 返回大小和存储参数。
mmread(源)	从类似于 Matrix Market 的 '源' 中读取内容到矩阵中。
mmwrite(目标, a[, 注释, 字段, ...])	将稀疏或密集数组 a 写入类似于 Matrix Market 的 '目标' 文件。

Wav 声音文件（scipy.io.wavfile）

read(文件名[, mmap])	打开 WAV 文件。
write(文件名, rate, 数据)	将 NumPy 数组写入 WAV 文件。

Arff 文件（scipy.io.arff）

loadarff(f)	读取 arff 文件。

Netcdf

netcdf_file(文件名[, 模式, mmap, 版本, ...])	用于 NetCDF 数据的文件对象。

允许读取 NetCDF 文件（使用 pupynere 包的版本）

可执行教程

插值过渡指南

原文：docs.scipy.org/doc/scipy-1.12.0/notebooks/interp_transition_guide.html

本笔记本包含三组演示：

• 用于遗留 bug 兼容scipy.interpolate.interp2d的 FITPACK 低级替代品；

• 建议用于新代码中的scipy.interpolate.interp2d替代品；

• 展示了基于 2D FITPACK 的线性插值失败模式及推荐的替代方案。

注意： 由于本笔记本展示了interp2d的用法（标记为已弃用），我们将简单起见静默处理弃用警告：

import warnings
warnings.filterwarnings('ignore') 

1. 如何过渡到不再使用interp2d

interp2d在 2D 常规网格上和插值 2D 分散数据时会静默切换。这种切换基于（拉直后的）x、y和z数组的长度。简而言之，对于常规网格，请使用scipy.interpolate.RectBivariateSpline；对于分散插值，请使用bisprep/bisplev组合。下面我们提供了逐点转换的文字示例，这应该完全保留interp2d的结果。

1.1 interp2d在常规网格上的应用

我们从（稍作修改的）文档字符串示例开始。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp2d, RectBivariateSpline

x = np.arange(-5.01, 5.01, 0.25)
y = np.arange(-5.01, 7.51, 0.25)
xx, yy = np.meshgrid(x, y)
z = np.sin(xx**2 + 2.*yy**2)
f = interp2d(x, y, z, kind='cubic') 

这是“常规网格”代码路径的示例，因为

z.size == len(x) * len(y) 

True 

还要注意x.size != y.size：

x.size, y.size 

(41, 51) 

现在，让我们构建一个方便的函数来构造插值器并绘制它。

def plot(f, xnew, ynew):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
    znew = f(xnew, ynew)

    ax1.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')

    im = ax2.imshow(znew)
    plt.colorbar(im, ax=ax2)

    plt.show()
    return znew 

绘图：

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 7.51, 1e-2)
znew_i = plot(f, xnew, ynew) 

替代方案：使用RectBivariateSpline，结果完全相同

注意转置：首先在构造函数中，其次在评估结果时需要转置。这是为了撤消interp2d的转置操作。

r = RectBivariateSpline(x, y, z.T)

rt = lambda xnew, ynew: r(xnew, ynew).T
znew_r = plot(rt, xnew, ynew) 

from numpy.testing import assert_allclose
assert_allclose(znew_i, znew_r, atol=1e-14) 

1.2. 使用点的完整坐标进行interp2d（分散插值）

在这里，我们展示了前一个练习中的网格平铺以说明功能。

xxr = xx.ravel()
yyr = yy.ravel()
zzr = z.ravel()

f = interp2d(xxr, yyr, zzr, kind='cubic') 

注意这是“非常规网格”代码路径，用于分散数据，其中len(x) == len(y) == len(z)。

len(xxr) == len(yyr) == len(zzr) 

True 

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 7.51, 1e-2)
znew_i = plot(f, xnew, ynew) 

替换：直接使用scipy.interpolate.bisplrep / scipy.interpolate.bisplev

from scipy.interpolate import bisplrep, bisplev
tck = bisplrep(xxr, yyr, zzr, kx=3, ky=3, s=0)
# convenience: make up a callable from bisplev
ff = lambda xnew, ynew: bisplev(xnew, ynew, tck).T   # Note the transpose, to mimic what interp2d does

znew_b = plot(ff, xnew, ynew) 

assert_allclose(znew_i, znew_b, atol=1e-15) 

2. 替代interp2d：正则网格

对于新代码，推荐的替代方案是RegularGridInterpolator。这是一个独立的实现，不基于 FITPACK。支持最近邻、线性插值和奇次张量积样条。

样条结节保证与数据点重合。

注意，这里：

1. 元组参数是(x, y)

2. z数组需要转置

3. 关键字名称是method，而不是kind

4. bounds_error参数默认为True。

from scipy.interpolate import RegularGridInterpolator as RGI

r = RGI((x, y), z.T, method='linear', bounds_error=False) 

评估：创建一个 2D 网格。使用indexing='ij'和sparse=True以节省一些内存：

xxnew, yynew = np.meshgrid(xnew, ynew, indexing='ij', sparse=True) 

评估时，请注意元组参数：

znew_reggrid = r((xxnew, yynew)) 

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))

# Again, note the transpose to undo the `interp2d` convention
znew_reggrid_t = znew_reggrid.T

ax1.plot(x, z[0, :], 'ro-', xnew, znew_reggrid_t[0, :], 'b-')

im = ax2.imshow(znew_reggrid_t)
plt.colorbar(im, ax=ax2) 

<matplotlib.colorbar.Colorbar at 0x7fa55ccd5f10> 

3. 散点 2D 线性插值：优先使用LinearNDInterpolator而不是SmoothBivariateSpline或bisplrep

对于 2D 散点线性插值，SmoothBivariateSpline和biplrep可能会发出警告，或者无法插值数据，或者产生带有结节远离数据点的样条。“相反，建议使用LinearNDInterpolator，它基于通过QHull对数据进行三角剖分。

# TestSmoothBivariateSpline::test_integral
from scipy.interpolate import SmoothBivariateSpline, LinearNDInterpolator

x = np.array([1,1,1,2,2,2,4,4,4])
y = np.array([1,2,3,1,2,3,1,2,3])
z = np.array([0,7,8,3,4,7,1,3,4]) 

现在，使用基于 Qhull 的数据三角剖分进行线性插值：

xy = np.c_[x, y]   # or just list(zip(x, y))
lut2 = LinearNDInterpolator(xy, z)

X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y) 

结果易于理解和解释：

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

ax.plot_wireframe(X, Y, lut2(X, Y))
ax.scatter(x, y, z,  'o', color='k', s=48) 

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa55cbd8250> 

注意，bisplrep做了一些不同的事情！它可能会将样条结节放在数据之外。

作为说明，考虑前面示例中的相同数据：

tck = bisplrep(x, y, z, kx=1, ky=1, s=0)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

xx = np.linspace(min(x), max(x))
yy = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xx, yy)
Z = bisplev(xx, yy, tck)
Z = Z.reshape(*X.shape).T

ax.plot_wireframe(X, Y, Z, rstride=2, cstride=2)
ax.scatter(x, y, z,  'o', color='k', s=48) 

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa55cc26310> 

此外，SmoothBivariateSpline无法插值数据。再次使用前面示例中的相同数据。

lut = SmoothBivariateSpline(x, y, z, kx=1, ky=1, s=0)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

xx = np.linspace(min(x), max(x))
yy = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xx, yy)

ax.plot_wireframe(X, Y, lut(xx, yy).T, rstride=4, cstride=4)
ax.scatter(x, y, z,  'o', color='k', s=48) 

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa55cc6ebe0> 

注意，SmoothBivariateSpline和bisplrep的结果都存在缺陷，不像LinearNDInterpolator那样。此处所示的问题是针对线性插值报告的，然而 FITPACK 的结节选择机制并不保证对高阶（如三次）样条曲面避免这些问题。

SciPy API

原文：docs.scipy.org/doc/scipy-1.12.0/reference/index.html

从 SciPy 导入

在 Python 中，库的公共 API 和私有实现细节之间的区分并不总是清晰的。与 Java 等其他语言不同，Python 中可以访问“私有”函数或对象。偶尔这可能很方便，但请注意，如果这样做，您的代码在未来版本中可能会无预警地中断。一些广泛认可的 Python 公共 API 规则包括：

• 方法/函数/类和模块属性名称以下划线开头的是私有的。

• 如果类名以下划线开头，则其所有成员都不是公共的，无论它们是否以下划线开头。

• 如果包中的模块名称以下划线开头，则其所有成员都不是公共的，无论它们是否以下划线开头。

• 如果模块或包定义了__all__，则这是官方定义的公共接口。

• 如果模块或包没有定义__all__，则所有不以下划线开头的名称都是公共的。

注意

阅读上述指南，可以得出结论，每个私有模块或对象都以下划线开头。但事实并非如此；下划线的存在确实标志着某些内容为私有，但缺少下划线并不意味着某些内容为公共的。

在 SciPy 中，有些模块的名称不以下划线开头，但应视为私有。为了澄清这些模块是哪些，我们在下面定义了 SciPy 的公共 API，并提供了一些关于如何从 SciPy 导入模块/函数/对象的建议。

从 SciPy 导入函数的指南

SciPy 子模块的命名空间中的所有内容都是公共的。通常在 Python 中建议使用命名空间。例如，函数 curve_fit（在 scipy/optimize/_minpack_py.py 中定义）应该这样导入：

import scipy
result = scipy.optimize.curve_fit(...) 

或者，您可以像这样使用子模块作为命名空间：

from scipy import optimize
result = optimize.curve_fit(...) 

注意

对于 scipy.io，推荐使用 import scipy，因为 io 也是 Python 标准库中的模块名称。

在某些情况下，公共 API 是更深层次的。例如，scipy.sparse.linalg 模块是公共的，它包含的函数在 scipy.sparse 命名空间中不可用。如果选择第二种形式，则代码更容易理解，例如，以下代码立即清楚地表明 lomax 是一个分布：

# first form
from scipy import stats
stats.lomax(...)

# second form
from scipy.stats import distributions
distributions.lomax(...) 

在这种情况下，如果文档中指明该子模块是公共的，则可以选择第二种形式。当然，您仍然可以使用：

import scipy
scipy.stats.lomax(...)
# or
scipy.stats.distributions.lomax(...) 

注意

SciPy 使用延迟加载机制，这意味着只有在首次尝试访问模块时才会将其加载到内存中。

注意

scipy 命名空间本身还包含从 numpy 导入的函数。这些函数仍然存在以保持向后兼容性，但应直接从 numpy 导入。

API 定义

下面列出的每个子模块都是公共的。这意味着这些子模块不太可能被重命名或以不兼容的方式进行更改，如果必须更改，将在 SciPy 的一个版本中引发弃用警告。

• scipy

• scipy.cluster

  • scipy.cluster.vq

  • scipy.cluster.hierarchy

• scipy.constants

• scipy.datasets

• scipy.fft

• scipy.fftpack

• scipy.integrate

• scipy.interpolate

• scipy.io

  • scipy.io.arff

  • scipy.io.matlab

  • scipy.io.wavfile

• scipy.linalg

  • scipy.linalg.blas

  • scipy.linalg.cython_blas

  • scipy.linalg.lapack

  • scipy.linalg.cython_lapack

  • scipy.linalg.interpolative

• scipy.misc

• scipy.ndimage

• scipy.odr

• scipy.optimize

  • scipy.optimize.cython_optimize

• scipy.signal

  • scipy.signal.windows

• scipy.sparse

  • scipy.sparse.linalg

  • scipy.sparse.csgraph

• scipy.spatial

  • scipy.spatial.distance

  • scipy.spatial.transform

• scipy.special

• scipy.stats

  • scipy.stats.contingency

  • scipy.stats.distributions

  • scipy.stats.mstats

  • scipy.stats.qmc

  • scipy.stats.sampling

SciPy 结构

所有 SciPy 模块应遵循以下约定。在此处，SciPy 模块 定义为位于 scipy/ 目录中的 Python 包，比如 yyy。

• 理想情况下，每个 SciPy 模块应尽可能自包含。即应最小化对其他包或模块的依赖。当然，假定对 NumPy 的依赖。

• 目录 yyy/ 包含：

  • 文件 meson.build 包含子模块的构建配置。

  • 目录 tests/ 包含文件 test_<name>.py，对应模块 yyy/<name>{.py,.so,/}。

• 私有模块应以下划线 _ 前缀，例如 yyy/_somemodule.py。

• 用户可见的函数应遵循 NumPy 文档风格。

• 模块的 __init__.py 应包含其主要参考文档，位于其 docstring 中。这与 Sphinx 文档在 doc/ 下的连接通过 Sphinx 的 automodule 指令。

  参考文档应首先使用 autosummary:: 指令列出模块内容的分类列表，随后解释了解模块使用的重要点。

  教程风格的文档应详细示例，需单独放置于 doc/source/tutorial/。

参见现有的 SciPy 子模块以获取指导。

SciPy 的主要命名空间

原文：docs.scipy.org/doc/scipy-1.12.0/reference/main_namespace.html

scipy 的主命名空间设计上仅包含很少的对象。只显示与测试、构建信息和版本控制相关的通用功能，以及一个类（LowLevelCallable(function[, user_data, ...])	低级回调函数。
show_config([mode])	显示构建和使用 SciPy 的库和系统信息
test	运行此命名空间的测试

唯一的公共属性是：

__version__	SciPy 版本字符串

子模块

cluster	聚类功能
constants	物理和数学常数及单位
datasets	载入 SciPy 数据集
fft	离散 Fourier 及相关变换
fftpack	离散 Fourier 变换（遗留）
integrate	数值积分和常微分方程组
interpolate	插值
io	科学数据格式读写
linalg	线性代数功能
misc	实用程序例程（已弃用）
ndimage	N 维图像处理和插值
odr	正交距离回归
optimize	数值优化
signal	信号处理
sparse	稀疏数组、线性代数和图算法
spatial	空间数据结构和算法
special	特殊函数
stats	统计函数

scipy.LowLevelCallable

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.LowLevelCallable.html#scipy.LowLevelCallable

class scipy.LowLevelCallable(function, user_data=None, signature=None)

低级回调函数。

SciPy 中的某些函数接受回调函数作为参数，这些函数可以是 Python 可调用对象或低级编译函数。使用编译的回调函数可以通过避免将数据包装在 Python 对象中来提高性能。

SciPy 中这种低级函数被包装在LowLevelCallable对象中，可以从 ctypes、cffi、Cython 获取的函数指针或包含在 Python PyCapsule对象中构造。

参见

接受低级可调用函数的函数：

scipy.integrate.quad, scipy.ndimage.generic_filter, scipy.ndimage.generic_filter1d, scipy.ndimage.geometric_transform

使用示例：

在 C 中扩展 scipy.ndimage, 使用低级回调函数加速积分

参数：

function

低级回调函数。

user_data

要传递到回调函数的用户数据。

signaturestr, 可选

函数的签名。如果省略，将从function中确定，如果可能的话。

注意

参数function可以是以下之一：

• 包含 C 函数签名的 PyCapsule

• ctypes 函数指针

• cffi 函数指针

低级回调函数的签名必须与其传递到的例程所期望的签名之一匹配。

如果从 PyCapsule 构造低级函数，则胶囊的名称必须是相应签名，格式为：

return_type (arg1_type, arg2_type, ...) 

例如：

"void (double)"
"double (double, int *, void *)" 

如果未显式提供user_data的值，则使用作为function传入的 PyCapsule 的上下文作为user_data。

属性：

function

给定回调函数。

user_data

给定的用户数据。

signature

函数的签名。

方法

from_cython(module, name[, user_data, signature])	从导出的 Cython 函数创建低级回调函数。

scipy.show_config

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.show_config.html#scipy.show_config

scipy.show_config(mode='stdout')

显示构建和使用 SciPy 的库和系统信息

参数：

mode{‘stdout’, ‘dicts’}，可选。

指示如何显示配置信息。‘stdout’ 输出到控制台，‘dicts’ 返回配置的字典。

返回：

out{dict, None}

如果 mode 是‘dicts’，则返回一个字典，否则返回 None

注意

1. 如果安装了pyyaml，‘stdout’ 模式将提供更可读的输出

scipy.test

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.test.html#scipy.test

scipy.test = <scipy._lib._testutils.PytestTester object>

运行此命名空间的测试

scipy.test() 运行 SciPy 的所有测试，使用默认设置。当从子模块使用时（例如 scipy.cluster.test()），仅运行该命名空间的测试。

参数：

label，可选

是否仅运行快速测试，还是包括标记为慢速的测试。默认为 'fast'。

verboseint，可选

测试输出详细程度，默认为 1。

extra_argvlist，可选

传递给 Pytest 的参数。

doctestsbool，可选

是否运行 doctests，默认为 False。

coveragebool，可选

是否启用代码覆盖率测量运行测试。默认为 False。

tests字符串列表，可选

要运行测试的模块名称列表。默认情况下，使用调用 test 函数的模块。

parallelint，可选

使用 pytest-xdist 并行运行测试，如果给定的数字大于 1。默认为 1。

聚类包（scipy.cluster）

原文：docs.scipy.org/doc/scipy-1.12.0/reference/cluster.html

在信息理论、目标检测、通信、压缩及其他领域中，聚类算法非常有用。vq 模块仅支持向量量化和 k-means 算法。

hierarchy 模块提供了用于层次聚类和凝聚聚类的函数。其功能包括从距离矩阵生成层次聚类，计算聚类的统计信息，切断链接以生成平坦聚类，并通过树状图可视化聚类。

常数（scipy.constants）

原文：docs.scipy.org/doc/scipy-1.12.0/reference/constants.html

物理和数学常数与单位。

数学常数

pi	圆周率 Pi
golden	黄金比例
golden_ratio	黄金比例

物理常数

c	真空中光速
speed_of_light	真空中光速
mu_0	磁常数 (\mu_0)
epsilon_0	电常数（真空介电常数），(\epsilon_0)
h	普朗克常数 (h)
Planck	普朗克常数 (h)
hbar	约化普朗克常数 (\hbar = h/(2\pi))
G	牛顿引力常数
gravitational_constant	牛顿引力常数
g	标准重力加速度
e	元电荷
elementary_charge	元电荷
R	摩尔气体常数
gas_constant	摩尔气体常数
alpha	精细结构常数
fine_structure	精细结构常数
N_A	阿伏伽德罗常数
Avogadro	阿伏伽德罗常数
k	玻尔兹曼常数
Boltzmann	玻尔兹曼常数
sigma	斯特藩-玻尔兹曼常数 (\sigma)
Stefan_Boltzmann	斯特藩-玻尔兹曼常数 (\sigma)
Wien	维恩位移定律常数
Rydberg	雷德伯常数
m_e	电子质量
electron_mass	电子质量
m_p	质子质量
proton_mass	质子质量
m_n	中子质量
neutron_mass	中子质量

常数数据库

除了上述变量外，scipy.constants 还包含 2018 年 CODATA 推荐值数据库 [CODATA2018] 中更多物理常数。

value(key)	物理常数索引键的值
unit(key)	物理常数单元，通过键索引
precision(key)	物理常数索引键的相对精度
find([sub, disp])	返回包含给定字符串的物理常数键列表
ConstantWarning	访问不再存在于当前 CODATA 数据集中的常数时的警告

scipy.constants.physical_constants

物理常数词典，格式为 physical_constants[name] = (value, unit, uncertainty)。

可用常数：

alpha particle mass	6.6446573357e-27 kg
alpha particle mass energy equivalent	5.9719201914e-10 J
alpha particle mass energy equivalent in MeV	3727.3794066 MeV
alpha particle mass in u	4.001506179127 u
alpha particle molar mass	0.0040015061777 kg mol^-1
alpha particle relative atomic mass	4.001506179127
alpha particle-electron mass ratio	7294.29954142
alpha particle-proton mass ratio	3.97259969009
Angstrom star	1.00001495e-10 m
atomic mass constant	1.6605390666e-27 kg
atomic mass constant energy equivalent	1.4924180856e-10 J
atomic mass constant energy equivalent in MeV	931.49410242 MeV
atomic mass unit-electron volt relationship	931494102.42 eV
atomic mass unit-hartree relationship	34231776.874 E_h
atomic mass unit-hertz relationship	2.25234271871e+23 Hz
atomic mass unit-inverse meter relationship	751300661040000.0 m^-1
atomic mass unit-joule relationship	1.4924180856e-10 J
atomic mass unit-kelvin relationship	10809540191600.0 K
atomic mass unit-kilogram relationship	1.6605390666e-27 kg
atomic unit of 1st hyperpolarizability	3.2063613061e-53 C³ m³ J^-2
atomic unit of 2nd hyperpolarizability	6.2353799905e-65 C⁴ m⁴ J^-3
atomic unit of action	1.054571817e-34 J s
atomic unit of charge	1.602176634e-19 C
atomic unit of charge density	1081202384570.0 C m^-3
atomic unit of current	0.00662361823751 A
atomic unit of electric dipole mom.	8.4783536255e-30 C m
atomic unit of electric field	514220674763.0 V m^-1
atomic unit of electric field gradient	9.7173624292e+21 V m^-2
atomic unit of electric polarizability	1.64877727436e-41 C² m² J^-1
atomic unit of electric potential	27.211386245988 V
atomic unit of electric quadrupole mom.	4.4865515246e-40 C m²
atomic unit of energy	4.3597447222071e-18 J
atomic unit of force	8.2387234983e-08 N
atomic unit of length	5.29177210903e-11 m
atomic unit of mag. dipole mom.	1.85480201566e-23 J T^-1
atomic unit of mag. flux density	235051.756758 T
atomic unit of magnetizability	7.8910366008e-29 J T^-2
atomic unit of mass	9.1093837015e-31 kg
atomic unit of momentum	1.9928519141e-24 kg m s^-1
atomic unit of permittivity	1.11265005545e-10 F m^-1
atomic unit of time	2.4188843265857e-17 s
atomic unit of velocity	2187691.26364 m s^-1
Avogadro constant	6.02214076e+23 mol^-1
Bohr magneton	9.2740100783e-24 J T^-1
Bohr magneton in eV/T	5.788381806e-05 eV T^-1
Bohr magneton in Hz/T	13996244936.1 Hz T^-1
Bohr magneton in inverse meter per tesla	46.686447783 m^-1 T^-1
Bohr magneton in K/T	0.67171381563 K T^-1
Bohr radius	5.29177210903e-11 m
Boltzmann constant	1.380649e-23 J K^-1
Boltzmann constant in eV/K	8.617333262e-05 eV K^-1
Boltzmann constant in Hz/K	20836619120.0 Hz K^-1
Boltzmann constant in inverse meter per kelvin	69.50348004 m^-1 K^-1
characteristic impedance of vacuum	376.730313668 ohm
classical electron radius	2.8179403262e-15 m
Compton wavelength	2.42631023867e-12 m
conductance quantum	7.748091729e-05 S
conventional value of ampere-90	1.00000008887 A
conventional value of coulomb-90	1.00000008887 C
conventional value of farad-90	0.9999999822 F
conventional value of henry-90	1.00000001779 H
conventional value of Josephson constant	483597900000000.0 Hz V^-1
conventional value of ohm-90	1.00000001779 ohm
conventional value of volt-90	1.00000010666 V
conventional value of von Klitzing constant	25812.807 ohm
conventional value of watt-90	1.00000019553 W
Cu x unit	1.00207697e-13 m
deuteron g factor	0.8574382338
deuteron mag. mom.	4.330735094e-27 J T^-1
deuteron mag. mom. to Bohr magneton ratio	0.000466975457
deuteron mag. mom. to nuclear magneton ratio	0.8574382338
deuteron mass	3.3435837724e-27 kg
deuteron mass energy equivalent	3.00506323102e-10 J
deuteron mass energy equivalent in MeV	1875.61294257 MeV
deuteron mass in u	2.013553212745 u
deuteron molar mass	0.00201355321205 kg mol^-1
deuteron relative atomic mass	2.013553212745
deuteron rms charge radius	2.12799e-15 m
deuteron-electron mag. mom. ratio	-0.0004664345551
deuteron-electron mass ratio	3670.48296788
deuteron-neutron mag. mom. ratio	-0.44820653
deuteron-proton mag. mom. ratio	0.30701220939
deuteron-proton mass ratio	1.99900750139
electron charge to mass quotient	-175882001076.0 C kg^-1
electron g factor	-2.00231930436256
electron gyromag. ratio	176085963023.0 s^-1 T^-1
electron gyromag. ratio in MHz/T	28024.9514242 MHz T^-1
electron mag. mom.	-9.2847647043e-24 J T^-1
electron mag. mom. anomaly	0.00115965218128
electron mag. mom. to Bohr magneton ratio	-1.00115965218128
electron mag. mom. to nuclear magneton ratio	-1838.28197188
electron mass	9.1093837015e-31 kg
electron mass energy equivalent	8.1871057769e-14 J
electron mass energy equivalent in MeV	0.51099895 MeV
electron mass in u	0.000548579909065 u
electron molar mass	5.4857990888e-07 kg mol^-1
electron relative atomic mass	0.000548579909065
electron to alpha particle mass ratio	0.0001370933554787
electron to shielded helion mag. mom. ratio	864.058257
electron to shielded proton mag. mom. ratio	-658.2275971
electron volt	1.602176634e-19 J
electron volt-atomic mass unit relationship	1.07354410233e-09 u
electron volt-hartree relationship	0.036749322175655 E_h
electron volt-hertz relationship	241798924200000.0 Hz
electron volt-inverse meter relationship	806554.3937 m^-1
electron volt-joule relationship	1.602176634e-19 J
electron volt-kelvin relationship	11604.51812 K
electron volt-kilogram relationship	1.782661921e-36 kg
electron-deuteron mag. mom. ratio	-2143.9234915
electron-deuteron mass ratio	0.0002724437107462
electron-helion mass ratio	0.0001819543074573
electron-muon mag. mom. ratio	206.7669883
electron-muon mass ratio	0.00483633169
electron-neutron mag. mom. ratio	960.9205
electron-neutron mass ratio	0.00054386734424
electron-proton mag. mom. ratio	-658.21068789
electron-proton mass ratio	0.000544617021487
electron-tau mass ratio	0.000287585
electron-triton mass ratio	0.0001819200062251
elementary charge	1.602176634e-19 C
elementary charge over h-bar	1519267447000000.0 A J^-1
Faraday constant	96485.33212 C mol^-1
Fermi coupling constant	1.1663787e-05 GeV^-2
fine-structure constant	0.0072973525693
first radiation constant	3.741771852e-16 W m²
first radiation constant for spectral radiance	1.191042972e-16 W m² sr^-1
Hartree energy	4.3597447222071e-18 J
Hartree energy in eV	27.211386245988 eV
hartree-atomic mass unit relationship	2.92126232205e-08 u
hartree-electron volt relationship	27.211386245988 eV
hartree-hertz relationship	6579683920502000.0 Hz
hartree-inverse meter relationship	21947463.13632 m^-1
hartree-joule relationship	4.3597447222071e-18 J
hartree-kelvin relationship	315775.02480407 K
hartree-kilogram relationship	4.8508702095432e-35 kg
helion g factor	-4.255250615
helion mag. mom.	-1.074617532e-26 J T^-1
helion mag. mom. to Bohr magneton ratio	-0.001158740958
helion mag. mom. to nuclear magneton ratio	-2.127625307
helion mass	5.0064127796e-27 kg
helion mass energy equivalent	4.4995394125e-10 J
helion mass energy equivalent in MeV	2808.39160743 MeV
helion mass in u	3.014932247175 u
helion molar mass	0.00301493224613 kg mol^-1
helion relative atomic mass	3.014932247175
helion shielding shift	5.996743e-05
helion-electron mass ratio	5495.88528007
helion-proton mass ratio	2.99315267167
hertz-atomic mass unit relationship	4.4398216652e-24 u
hertz-electron volt relationship	4.135667696e-15 eV
hertz-hartree relationship	1.519829846057e-16 E_h
hertz-inverse meter relationship	3.3356409519815204e-09 m^-1
hertz-joule relationship	6.62607015e-34 J
hertz-kelvin relationship	4.799243073e-11 K
hertz-kilogram relationship	7.372497323e-51 kg
hyperfine transition frequency of Cs-133	9192631770.0 Hz
inverse fine-structure constant	137.035999084
inverse meter-atomic mass unit relationship	1.3310250501e-15 u
inverse meter-electron volt relationship	1.239841984e-06 eV
inverse meter-hartree relationship	4.556335252912e-08 E_h
inverse meter-hertz relationship	299792458.0 Hz
每米-焦耳关系	1.986445857e-25 J
每米-开尔文关系	0.01438776877 K
每米-千克关系	2.210219094e-42 kg
电导量子的倒数	12906.40372 ohm
Josephson 常数	483597848400000.0 Hz V^-1
焦耳-原子质量单位关系	6700535256.5 u
焦耳-电子伏特关系	6.241509074e+18 eV
焦耳-哈特里关系	2.2937122783963e+17 E_h
焦耳-赫兹关系	1.509190179e+33 Hz
焦耳-每米关系	5.034116567e+24 m^-1
焦耳-开尔文关系	7.242970516e+22 K
焦耳-千克关系	1.1126500560536185e-17 kg
开尔文-原子质量单位关系	9.2510873014e-14 u
开尔文-电子伏特关系	8.617333262e-05 eV
开尔文-哈特里关系	3.1668115634556e-06 E_h
开尔文-赫兹关系	20836619120.0 Hz
开尔文-每米关系	69.50348004 m^-1
开尔文-焦耳关系	1.380649e-23 J
开尔文-千克关系	1.536179187e-40 kg
千克-原子质量单位关系	6.0221407621e+26 u
千克-电子伏特关系	5.609588603e+35 eV
千克-哈特里关系	2.0614857887409e+34 E_h
千克-赫兹关系	1.356392489e+50 Hz
千克-每米关系	4.524438335e+41 m^-1
千克-焦耳关系	8.987551787368176e+16 J
千克-开尔文关系	6.50965726e+39 K
硅晶格常数	5.431020511e-10 m
理想 Si 的晶格间距 (220)	1.920155716e-10 m
Loschmidt 常数 (273.15 K, 100 kPa)	2.651645804e+25 m^-3
Loschmidt 常数 (273.15 K, 101.325 kPa)	2.686780111e+25 m^-3
光效	683.0 lm W^-1
磁通量子	2.067833848e-15 Wb
Mo x 单位	1.00209952e-13 m
摩尔气体常数	8.314462618 J mol^-1 K^-1
摩尔质量常数	0.00099999999965 kg mol^-1
碳-12 的摩尔质量	0.0119999999958 kg mol^-1
摩尔普朗克常数	3.990312712e-10 J Hz^-1 mol^-1
理想气体的摩尔体积 (273.15 K, 100 kPa)	0.02271095464 m³ mol^-1
理想气体的摩尔体积 (273.15 K, 101.325 kPa)	0.02241396954 m³ mol^-1
硅的摩尔体积	1.205883199e-05 m³ mol^-1
μ子康普顿波长	1.17344411e-14 m
μ子 g 因子	-2.0023318418
μ子磁矩	-4.4904483e-26 J T^-1
μ子磁矩异常	0.00116592089
μ子磁矩对玻尔磁子比	-0.00484197047
μ子磁矩对核磁子比	-8.89059703
μ子质量	1.883531627e-28 kg
μ子质能等效	1.692833804e-11 J
μ子质能等效（以 MeV 计）	105.6583755 MeV
μ子质量（以 u 计）	0.1134289259 u
μ子的摩尔质量	0.0001134289259 kg mol^-1
μ子-电子质量比	206.768283
μ子-中子质量比	0.112454517
μ子-质子磁矩比	-3.183345142
μ子-质子质量比	0.1126095264
μ子-τ子质量比	0.0594635
电子伏秒中的自然单位行动	1.054571817e-34 J s
电子伏-秒中的自然单位行动	6.582119569e-16 eV s
自然能量单位	8.1871057769e-14 J
以 MeV 为单位的自然能量单位	0.51099895 MeV
自然长度单位	3.8615926796e-13 m
自然质量单位	9.1093837015e-31 kg
自然动量单位	2.730924488e-22 kg m s^-1
以 MeV/c 为单位的自然单位动量	0.5109989461 MeV/c
自然时间单位	1.28808866819e-21 s
以 m/s 为单位的自然速度单位	299792458.0 m s^-1
中子的康普顿波长	1.31959090581e-15 m
中子的 g 因子	-3.82608545
中子的旋磁比	183247171.0 s^-1 T^-1
中子的旋磁比（以 MHz/T 为单位）	29.1646931 MHz T^-1
中子的磁矩	-9.6623651e-27 J T^-1
中子磁矩对玻尔磁子比值	-0.00104187563
中子磁矩对核磁子比值	-1.91304273
中子的质量	1.67492749804e-27 kg
中子的质能等效	1.50534976287e-10 J
中子的质能等效（以 MeV 为单位）	939.56542052 MeV
中子的质量（以 u 为单位）	1.00866491595 u
中子的摩尔质量	0.0010086649156 kg mol^-1
中子的相对原子质量	1.00866491595
中子到屏蔽质子磁矩比	-0.68499694
中子-电子磁矩比	0.00104066882
中子-电子质量比	1838.68366173
中子-μ子质量比	8.89248406
中子-质子磁矩比	-0.68497934
中子-质子质量差的能量等效	2.30557435e-30 kg
中子-质子质量差的能量等效	2.07214689e-13 J
中子-质子质量差的能量等效（以 MeV 为单位）	1.29333236 MeV
中子-质子质量差的 u 值	0.00138844933 u
中子-质子质量比	1.00137841931
中子-τ子质量比	0.528779
万有引力常数	6.6743e-11 m³ kg^-1 s^-2
以 GeV/c² 为单位的万有引力常数与 h-bar c 之比	6.70883e-39 (GeV/c²)^-2
核磁子	5.0507837461e-27 J T^-1
每电子伏每特斯拉中的核磁子	3.15245125844e-08 eV T^-1
每米每特斯拉中的核磁子	0.0254262341353 m^-1 T^-1
核磁子在 K/T 中的值	0.00036582677756 K T^-1
每 MHz 每特斯拉中的核磁子	7.6225932291 MHz T^-1
普朗克常数	6.62607015e-34 J Hz^-1
以 eV/Hz 为单位的普朗克常数	4.135667696e-15 eV Hz^-1
普朗克长度	1.616255e-35 m
普朗克质量	2.176434e-08 kg
以 GeV 为单位的普朗克质量能量等效	1.22089e+19 GeV
普朗克温度	1.416784e+32 K
普朗克时间	5.391247e-44 s
质子电荷与质量的比值	95788331.56 C kg^-1
质子的康普顿波长	1.32140985539e-15 m
质子的 g 因子	5.5856946893
质子的旋磁比（以 MHz/T 为单位）	267522187.44 s^-1 T^-1
质子的旋磁比（以 MHz/T 为单位）	42.577478518 MHz T^-1
proton mag. mom.	1.41060679736e-26 J T^-1
proton mag. mom. to Bohr magneton ratio	0.0015210322023
proton mag. mom. to nuclear magneton ratio	2.79284734463
proton mag. shielding correction	2.5689e-05
proton mass	1.67262192369e-27 kg
proton mass energy equivalent	1.50327761598e-10 J
proton mass energy equivalent in MeV	938.27208816 MeV
proton mass in u	1.007276466621 u
proton molar mass	0.00100727646627 kg mol^-1
proton relative atomic mass	1.007276466621
proton rms charge radius	8.414e-16 m
proton-electron mass ratio	1836.15267343
proton-muon mass ratio	8.88024337
proton-neutron mag. mom. ratio	-1.45989805
proton-neutron mass ratio	0.99862347812
proton-tau mass ratio	0.528051
quantum of circulation	0.00036369475516 m² s^-1
quantum of circulation times 2	0.00072738951032 m² s^-1
reduced Compton wavelength	3.8615926796e-13 m
reduced muon Compton wavelength	1.867594306e-15 m
reduced neutron Compton wavelength	2.1001941552e-16 m
reduced Planck constant	1.054571817e-34 J s
reduced Planck constant in eV s	6.582119569e-16 eV s
reduced Planck constant times c in MeV fm	197.3269804 MeV fm
reduced proton Compton wavelength	2.10308910336e-16 m
reduced tau Compton wavelength	1.110538e-16 m
Rydberg constant	10973731.56816 m^-1
Rydberg constant times c in Hz	3289841960250800.0 Hz
Rydberg constant times hc in eV	13.605693122994 eV
Rydberg constant times hc in J	2.1798723611035e-18 J
Sackur-Tetrode constant (1 K, 100 kPa)	-1.15170753706
Sackur-Tetrode constant (1 K, 101.325 kPa)	-1.16487052358
second radiation constant	0.01438776877 m K
shielded helion gyromag. ratio	203789456.9 s^-1 T^-1
shielded helion gyromag. ratio in MHz/T	32.43409942 MHz T^-1
shielded helion mag. mom.	-1.07455309e-26 J T^-1
shielded helion mag. mom. to Bohr magneton ratio	-0.001158671471
shielded helion mag. mom. to nuclear magneton ratio	-2.127497719
shielded helion to proton mag. mom. ratio	-0.7617665618
shielded helion to shielded proton mag. mom. ratio	-0.7617861313
shielded proton gyromag. ratio	267515315.1 s^-1 T^-1
shielded proton gyromag. ratio in MHz/T	42.57638474 MHz T^-1
shielded proton mag. mom.	1.41057056e-26 J T^-1
shielded proton mag. mom. to Bohr magneton ratio	0.001520993128
shielded proton mag. mom. to nuclear magneton ratio	2.792775599
shielding difference of d and p in HD	2.02e-08
shielding difference of t and p in HT	2.414e-08
speed of light in vacuum	299792458.0 m s^-1
standard acceleration of gravity	9.80665 m s^-2
standard atmosphere	101325.0 Pa
standard-state pressure	100000.0 Pa
Stefan-Boltzmann constant	5.670374419e-08 W m^-2 K^-4
tau Compton wavelength	6.97771e-16 米
tau energy equivalent	1776.86 兆电子伏特
tau mass	3.16754e-27 kg
tau mass energy equivalent	2.84684e-10 焦耳
tau mass in u	1.90754 u
tau molar mass	0.00190754 千克 mol^-1
tau-electron mass ratio	3477.23
tau-muon mass ratio	16.817
tau-neutron mass ratio	1.89115
tau-proton mass ratio	1.89376
Thomson cross section	6.6524587321e-29 平方米
triton g factor	5.957924931
triton mag. mom.	1.5046095202e-26 焦耳特^-1
triton mag. mom. to Bohr magneton ratio	0.0016223936651
triton mag. mom. to nuclear magneton ratio	2.9789624656
triton mass	5.0073567446e-27 千克
triton mass energy equivalent	4.500387806e-10 焦耳
triton mass energy equivalent in MeV	2808.92113298 兆电子伏特
triton mass in u	3.01550071621 u
triton molar mass	0.00301550071517 千克 mol^-1
triton relative atomic mass	3.01550071621
triton to proton mag. mom. ratio	1.0666399191
triton-electron mass ratio	5496.92153573
triton-proton mass ratio	2.99371703414
unified atomic mass unit	1.6605390666e-27 千克
vacuum electric permittivity	8.8541878128e-12 法拉每米
vacuum mag. permeability	1.25663706212e-06 牛顿安培^-2
von Klitzing constant	25812.80745 欧姆
W to Z mass ratio	0.88153
weak mixing angle	0.2229
Wien frequency displacement law constant	58789257570.0 赫兹开尔文^-1
Wien wavelength displacement law constant	0.002897771955 米开尔文

单位

SI 前缀

quetta	(10^{30})
ronna	(10^{27})
yotta	(10^{24})
zetta	(10^{21})
exa	(10^{18})
peta	(10^{15})
tera	(10^{12})
giga	(10^{9})
mega	(10^{6})
kilo	(10^{3})
hecto	(10^{2})
deka	(10^{1})
deci	(10^{-1})
centi	(10^{-2})
milli	(10^{-3})
micro	(10^{-6})
nano	(10^{-9})
pico	(10^{-12})
femto	(10^{-15})
atto	(10^{-18})
zepto	(10^{-21})
yocto	(10^{-24})
ronto	(10^{-27})
quecto	(10^{-30})

二进制前缀

kibi	(2^{10})
mebi	(2^{20})
gibi	(2^{30})
tebi	(2^{40})
pebi	(2^{50})
exbi	(2^{60})
zebi	(2^{70})
yobi	(2^{80})

质量

gram	(10^{-3}) 千克
metric_ton	(10^{3}) 千克
grain	一粒在千克中的质量
lb	一磅（常衡制）在千克中的质量
pound	一磅（常衡制）在千克中的质量
blob	一英寸版本的蛞蝓在千克中的质量（自 1.0.0 版本添加）
slinch	一英寸版本的蛞蝓在千克中的质量（自 1.0.0 版本添加）
slug	一只蛞蝓在千克中的质量（自 1.0.0 版本添加）
oz	一盎司在千克中的质量
ounce	一盎司在千克中的质量
stone	一英石在千克中的质量
grain	一粒在千克中的质量
long_ton	一长吨在千克中的质量
short_ton	一短吨在千克中的质量
troy_ounce	一金衡盎司在千克中的质量
troy_pound	一金衡磅的千克数
carat	一克拉的千克数
m_u	原子质量常数（单位：千克）
u	原子质量常数（单位：千克）
atomic_mass	原子质量常数（单位：千克）

角度

degree	度的弧度数
arcmin	角分的弧度数
arcminute	角分的弧度数
arcsec	角秒的弧度数
arcsecond	角秒的弧度数

时间

minute	一分钟的秒数
hour	一小时的秒数
day	一天的秒数
week	一周的秒数
year	一年（365 天）的秒数
Julian_year	一儒略年（365.25 天）的秒数

长度

inch	一英寸的米数
foot	一英尺的米数
yard	一码的米数
mile	一英里的米数
mil	一毫的米数
pt	一点的米数
point	一点的米数
survey_foot	一测量英尺的米数
survey_mile	一测量英里的米数
nautical_mile	一海里的米数
fermi	一费米的米数
angstrom	一埃的米数
micron	一微米的米数
au	一天文单位的米数
astronomical_unit	一天文单位的米数
light_year	一光年的米数
parsec	一秒差距的米数

压力

atm	标准大气压的帕斯卡数
atmosphere	标准大气压的帕斯卡数
bar	一巴的帕斯卡数
torr	一托（毫米汞柱）的帕斯卡数
mmHg	一托（毫米汞柱）的帕斯卡数
psi	一磅力每平方英寸的帕斯卡数

面积

hectare	一公顷的平方米数
acre	一英亩的平方米数

体积

liter	一升的立方米数
litre	一升的立方米数
gallon	一加仑（美国）的立方米数
gallon_US	一加仑（美国）的立方米数
gallon_imp	一加仑（英国）的立方米数
fluid_ounce	一液体盎司（美国）的立方米数
fluid_ounce_US	一液体盎司（美国）的立方米数
fluid_ounce_imp	一液体盎司（英国）的立方米数
bbl	一桶的立方米数
barrel	一桶的立方米数

速度

kmh	千米每小时的米每秒数
mph	英里每小时的米每秒数
mach	一马赫（近似值，15 摄氏度，1 标准大气压）的米每秒数
speed_of_sound	一马赫（近似值，15 摄氏度，1 标准大气压）的米每秒数
knot	一节的米每秒数

温度

zero_Celsius	摄氏零度对应的开尔文数
degree_Fahrenheit	一华氏度（仅温差）的开尔文数
convert_temperature（val, old_scale, new_scale）	将温度从一个温度标度转换到另一个（包括摄氏度、开尔文、华氏度和兰氏度）的函数。

能量

eV	一电子伏特的焦耳数
electron_volt	一电子伏特的焦耳数
calorie	一卡（热化学）的焦耳数
calorie_th	一卡路里（热化学）等于多少焦耳
calorie_IT	一卡路里（国际蒸汽表卡路里，1956 年）等于多少焦耳
erg	一爱尔格等于多少焦耳
Btu	一英热单位（国际蒸汽表）等于多少焦耳
Btu_IT	一英热单位（国际蒸汽表）等于多少焦耳
Btu_th	一英热单位（热化学）等于多少焦耳
ton_TNT	一吨 TNT 等于多少焦耳

功率

hp	一马力等于多少瓦特
horsepower	一马力等于多少瓦特

力量

dyn	一达因等于多少牛顿
dyne	一达因等于多少牛顿
lbf	一磅力等于多少牛顿
pound_force	一磅力等于多少牛顿
kgf	一公斤力等于多少牛顿
kilogram_force	一公斤力等于多少牛顿

光学

lambda2nu(lambda_)	将波长转换为光学频率
nu2lambda(nu)	将光学频率转换为波长。

参考文献

[CODATA2018]

CODATA 2018 年推荐的基础物理常数。

physics.nist.gov/cuu/Constants/

scipy.constants.value

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.value.html#scipy.constants.value

scipy.constants.value(key)

由键索引的 physical_constants 中的值

参数：

keyPython 字符串

字典中的键 physical_constants

返回：

value浮点数

key 对应的 physical_constants 中的值

示例

>>> from scipy import constants
>>> constants.value('elementary charge')
1.602176634e-19 

scipy.constants.unit

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.unit.html#scipy.constants.unit

scipy.constants.unit(key)

字典中由关键字索引的单位

参数：

关键Python 字符串

字典中的关键字 physical_constants

返回：

单位Python 字符串

对应于 关键字 的单位 physical_constants

示例

>>> from scipy import constants
>>> constants.unit('proton mass')
'kg' 

scipy.constants.precision

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.precision.html#scipy.constants.precision

scipy.constants.precision(key)

根据 key 索引的 physical_constants 中的相对精度

参数:

keyPython 字符串

字典physical_constants中的键

返回:

prec浮点数

相对精度，在physical_constants中对应于 key

示例

>>> from scipy import constants
>>> constants.precision('proton mass')
5.1e-37 

scipy.constants.find

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.find.html#scipy.constants.find

scipy.constants.find(sub=None, disp=False)

返回包含给定字符串的physical_constant键列表。

参数：

sub字符串

要搜索键的子字符串。默认情况下，返回所有键。

disp布尔

如果为 True，则打印找到的键并返回 None。否则，返回不打印任何内容的键列表。

返回：

keys列表或 None

如果disp为 False，则返回键列表。否则，返回 None。

示例

>>> from scipy.constants import find, physical_constants 

哪些键在physical_constants字典中包含‘boltzmann’？

>>> find('boltzmann')
['Boltzmann constant',
 'Boltzmann constant in Hz/K',
 'Boltzmann constant in eV/K',
 'Boltzmann constant in inverse meter per kelvin',
 'Stefan-Boltzmann constant'] 

获取名为‘Boltzmann constant in Hz/K’的常数：

>>> physical_constants['Boltzmann constant in Hz/K']
(20836619120.0, 'Hz K^-1', 0.0) 

查找键中包含‘radius’的常数：

>>> find('radius')
['Bohr radius',
 'classical electron radius',
 'deuteron rms charge radius',
 'proton rms charge radius']
>>> physical_constants['classical electron radius']
(2.8179403262e-15, 'm', 1.3e-24) 

scipy.constants.ConstantWarning

docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.ConstantWarning.html#scipy.constants.ConstantWarning

exception scipy.constants.ConstantWarning

访问一个不再包含在当前 CODATA 数据集中的常数。

with_traceback()

Exception.with_traceback(tb) – 将 self.traceback 设置为 tb 并返回 self。

scipy.constants.convert_temperature

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.convert_temperature.html#scipy.constants.convert_temperature

scipy.constants.convert_temperature(val, old_scale, new_scale)

将温度刻度从摄氏度、开尔文、华氏度和兰金刻度之一转换为另一个刻度。

参数：

valarray_like

将要转换的温度值（或数组）以原始刻度表示的数值。

old_scalestr

指定原始刻度的字符串，温度值将从中进行转换。支持的刻度有摄氏度（‘Celsius’、‘celsius’、‘C’或‘c’）、开尔文（‘Kelvin’、‘kelvin’、‘K’、‘k’）、华氏度（‘Fahrenheit’、‘fahrenheit’、‘F’或‘f’）和兰金（‘Rankine’、‘rankine’、‘R’或‘r’）。

new_scalestr

指定将温度值转换为的新刻度的字符串。支持的刻度有摄氏度（‘Celsius’、‘celsius’、‘C’或‘c’）、开尔文（‘Kelvin’、‘kelvin’、‘K’、‘k’）、华氏度（‘Fahrenheit’、‘fahrenheit’、‘F’或‘f’）和兰金（‘Rankine’、‘rankine’、‘R’或‘r’）。

返回：

resfloat or array of floats

转换后的温度值以新刻度表达。

注意

从版本 0.18.0 开始新增。

示例

>>> from scipy.constants import convert_temperature
>>> import numpy as np
>>> convert_temperature(np.array([-40, 40]), 'Celsius', 'Kelvin')
array([ 233.15,  313.15]) 

scipy.constants.lambda2nu。

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.lambda2nu.html#scipy.constants.lambda2nu

scipy.constants.lambda2nu(lambda_)

将波长转换为光学频率。

参数：

lambda_数组形式

波长需要转换。

返回：

nu浮点数或浮点数数组

等效光学频率。

注意事项。

计算nu = c / lambda，其中 c = 299792458.0，即真空中的光速，单位为米/秒。

示例

>>> from scipy.constants import lambda2nu, speed_of_light
>>> import numpy as np
>>> lambda2nu(np.array((1, speed_of_light)))
array([  2.99792458e+08,   1.00000000e+00]) 

scipy.constants.nu2lambda

原文：docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.constants.nu2lambda.html#scipy.constants.nu2lambda

scipy.constants.nu2lambda(nu)

将光频率转换为波长。

参数：

nu类似数组

要转换的光频率。

返回值：

lambda浮点数或浮点数数组

等效波长（们）。

注意事项

计算 lambda = c / nu，其中 c = 299792458.0，即真空中的光速（米/秒）。

示例

>>> from scipy.constants import nu2lambda, speed_of_light
>>> import numpy as np
>>> nu2lambda(np.array((1, speed_of_light)))
array([  2.99792458e+08,   1.00000000e+00])