SciPy 快速指南

---

---

SciPy - 简介

SciPy，读作“Sigh Pi”，是一个基于Python的开源科学计算库，根据BSD许可证发布，用于执行数学、科学和工程计算。

SciPy 库依赖于 NumPy，NumPy 提供了方便快捷的 N 维数组操作。SciPy 库构建在 NumPy 数组之上，并提供了许多用户友好且高效的数值计算方法，例如数值积分和优化的例程。它们都可以在所有流行的操作系统上运行，安装快速且免费。NumPy 和 SciPy 易于使用，但功能强大，足以被世界上一些领先的科学家和工程师所依赖。

SciPy 子包

SciPy 被组织成涵盖不同科学计算领域的子包。这些子包在下面的表格中进行了总结：

数据结构

SciPy 使用的基本数据结构是由 NumPy 模块提供的多维数组。NumPy 提供了一些用于线性代数、傅里叶变换和随机数生成的函数，但其通用性不如 SciPy 中的等效函数。

SciPy - 环境设置

标准 Python 发行版不包含任何 SciPy 模块。一种轻量级的替代方法是使用流行的 Python 包安装程序安装 SciPy，

pip install pandas

如果我们安装了**Anaconda Python 包**，Pandas 将默认安装。以下是不同操作系统中安装它们的软件包和链接。

Windows

**Anaconda**（来自 https://www.continuum.io）是 SciPy 栈的免费 Python 发行版。它也适用于 Linux 和 Mac。

**Canopy** (https://www.enthought.com/products/canopy/) 提供免费版和商业版，其中包含适用于 Windows、Linux 和 Mac 的完整 SciPy 栈。

**Python(x,y)** - 这是一个带有 SciPy 栈和 Spyder IDE 的免费 Python 发行版，适用于 Windows 操作系统。（可从 https://python-xy.github.io/ 下载）

Linux

使用各个 Linux 发行版的包管理器来安装 SciPy 栈中的一个或多个包。

Ubuntu

我们可以使用以下路径在 Ubuntu 中安装 Python。

sudo apt-get install python-numpy python-scipy 
python-matplotlibipythonipython-notebook python-pandas python-sympy python-nose

Fedora

我们可以使用以下路径在 Fedora 中安装 Python。

sudo yum install numpyscipy python-matplotlibipython python-pandas 
sympy python-nose atlas-devel

SciPy - 基本功能

默认情况下，所有 NumPy 函数都可通过 SciPy 命名空间访问。导入 SciPy 后，无需显式导入 NumPy 函数。NumPy 的主要对象是同质多维数组。它是一个元素（通常是数字）表，所有元素类型相同，由正整数元组索引。在 NumPy 中，维度称为轴。轴的数量称为秩。

现在，让我们回顾一下 NumPy 中向量和矩阵的基本功能。由于 SciPy 建立在 NumPy 数组之上，因此了解 NumPy 基础知识是必要的，因为线性代数的大部分内容都只处理矩阵。

NumPy 向量

向量可以通过多种方式创建。下面描述其中一些方法。

将 Python 类数组对象转换为 NumPy 数组

让我们考虑以下示例。

import numpy as np
list = [1,2,3,4]
arr = np.array(list)
print arr

上述程序的输出如下所示。

[1 2 3 4]

NumPy 本身的数组创建

NumPy 具有从头开始创建数组的内置函数。下面解释其中一些函数。

使用 zeros()

zeros(shape) 函数将创建一个用 0 值填充的数组，其形状由参数指定。默认 dtype 为 float64。让我们考虑以下示例。

import numpy as np
print np.zeros((2, 3))

上述程序的输出如下所示。

array([[ 0., 0., 0.],
[ 0., 0., 0.]])

使用 ones()

ones(shape) 函数将创建一个用 1 值填充的数组。在其他方面与 zeros 相同。让我们考虑以下示例。

import numpy as np
print np.ones((2, 3))

上述程序的输出如下所示。

array([[ 1., 1., 1.],
[ 1., 1., 1.]])

使用 arange()

arange() 函数将创建具有规律递增值的数组。让我们考虑以下示例。

import numpy as np
print np.arange(7)

上述程序将生成以下输出。

array([0, 1, 2, 3, 4, 5, 6])

定义值的 数据类型

让我们考虑以下示例。

import numpy as np
arr = np.arange(2, 10, dtype = np.float)
print arr
print "Array Data Type :",arr.dtype

上述程序将生成以下输出。

[ 2. 3. 4. 5. 6. 7. 8. 9.]
Array Data Type : float64

使用 linspace()

linspace() 函数将创建具有指定数量元素的数组，这些元素将均匀地分布在指定的起始值和结束值之间。让我们考虑以下示例。

import numpy as np
print np.linspace(1., 4., 6)

上述程序将生成以下输出。

array([ 1. , 1.6, 2.2, 2.8, 3.4, 4. ])

矩阵

矩阵是一个特殊的二维数组，它在运算过程中保持其二维特性。它具有一些特殊的运算符，例如 *（矩阵乘法）和 **（矩阵幂）。让我们考虑以下示例。

import numpy as np
print np.matrix('1 2; 3 4')

上述程序将生成以下输出。

matrix([[1, 2],
[3, 4]])

矩阵的共轭转置

此功能返回 `self` 的（复数）共轭转置。让我们考虑以下示例。

import numpy as np
mat = np.matrix('1 2; 3 4')
print mat.H

上述程序将生成以下输出。

matrix([[1, 3],
        [2, 4]])

矩阵的转置

此功能返回 `self` 的转置。让我们考虑以下示例。

import numpy as np
mat = np.matrix('1 2; 3 4')
mat.T

上述程序将生成以下输出。

matrix([[1, 3],
        [2, 4]])

当我们转置矩阵时，我们创建一个新的矩阵，其行是原始矩阵的列。另一方面，共轭转置为每个矩阵元素交换行和列索引。矩阵的逆矩阵是一个矩阵，如果它与原始矩阵相乘，则结果为单位矩阵。

SciPy - 聚类分析

**K 均值聚类**是一种在未标记数据集中查找聚类和聚类中心的方法。直观地说，我们可能会将一个聚类理解为——由一组数据点组成，其点与点之间的距离与到聚类外部点的距离相比很小。给定一组初始的 K 个中心，K 均值算法迭代以下两个步骤：

• 对于每个中心，识别比任何其他中心都更靠近它的训练点子集（其聚类）。

• 计算每个聚类中数据点每个特征的均值，该均值向量成为该聚类的新的中心。

这两个步骤重复进行，直到中心不再移动或分配不再改变。然后，可以将新点**x**分配到最接近原型的聚类。SciPy 库通过 cluster 包提供了一个很好的 K 均值算法实现。让我们了解如何使用它。

SciPy 中的 K 均值实现

我们将了解如何在 SciPy 中实现 K 均值。

导入 K 均值

我们将看到每个导入函数的实现和用法。

from SciPy.cluster.vq import kmeans,vq,whiten

数据生成

我们必须模拟一些数据来探索聚类。

from numpy import vstack,array
from numpy.random import rand

# data generation with three features
data = vstack((rand(100,3) + array([.5,.5,.5]),rand(100,3)))

现在，我们必须检查数据。上述程序将生成以下输出。

array([[ 1.48598868e+00, 8.17445796e-01, 1.00834051e+00],
       [ 8.45299768e-01, 1.35450732e+00, 8.66323621e-01],
       [ 1.27725864e+00, 1.00622682e+00, 8.43735610e-01],
       …………….

按特征基础对一组观察值进行归一化。在运行 K 均值之前，最好使用白化重新缩放观察集的每个特征维度。每个特征除以其在所有观察值中的标准差，使其方差为单位。

对数据进行白化处理

我们必须使用以下代码对数据进行白化处理。

# whitening of data
data = whiten(data)

使用三个聚类计算 K 均值

现在，让我们使用以下代码使用三个聚类计算 K 均值。

# computing K-Means with K = 3 (2 clusters)
centroids,_ = kmeans(data,3)

上述代码对形成 K 个聚类的一组观察向量执行 K 均值。K 均值算法调整质心，直到无法取得足够的进展，即自上次迭代以来失真变化小于某个阈值。在这里，我们可以通过使用下面给出的代码打印 centroids 变量来观察聚类的质心。

print(centroids)

上述代码将生成以下输出。

print(centroids)[ [ 2.26034702  1.43924335  1.3697022 ]
                  [ 2.63788572  2.81446462  2.85163854]
                  [ 0.73507256  1.30801855  1.44477558] ]

使用下面给出的代码将每个值分配给一个聚类。

# assign each sample to a cluster
clx,_ = vq(data,centroids)

**vq** 函数将 `M` x `N` 的 **obs** 数组中的每个观察向量与质心进行比较，并将观察值分配给最接近的聚类。它返回每个观察值的聚类和失真。我们也可以检查失真。让我们使用以下代码检查每个观察值的聚类。

# check clusters of observation
print clx

上述代码将生成以下输出。

array([1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 2, 0, 2, 0, 1, 1, 1,
0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0,
0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 1,  0, 0, 0, 0, 1, 0, 0, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 0, 0,
2, 2, 2, 1, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2], dtype=int32)

上述数组的不同值 0、1、2 表示聚类。

SciPy - 常量

SciPy 常量包提供了广泛的常数，这些常数用于一般的科学领域。

SciPy 常量包

**scipy.constants 包**提供了各种常数。我们必须导入所需的常数并根据需要使用它们。让我们看看这些常量变量是如何导入和使用的。

首先，让我们通过考虑以下示例来比较 'pi' 值。

#Import pi constant from both the packages
from scipy.constants import pi
from math import pi

print("sciPy - pi = %.16f"%scipy.constants.pi)
print("math - pi = %.16f"%math.pi)

上述程序将生成以下输出。

sciPy - pi = 3.1415926535897931
math - pi = 3.1415926535897931

可用常量列表

下表简要描述了各种常数。

数学常数

物理常数

下表列出了最常用的物理常数。

单位

下表列出了 SI 单位。

这些单位范围从尧、泽、艾、拍、太……千、百……纳、皮……到仄。

其他重要常数

下表列出了 SciPy 中使用的其他重要常数。

记住所有这些有点困难。获取哪个键对应哪个函数的简单方法是使用 **scipy.constants.find()** 方法。让我们考虑以下示例。

import scipy.constants
res = scipy.constants.physical_constants["alpha particle mass"]
print res

上述程序将生成以下输出。

[
   'alpha particle mass',
   'alpha particle mass energy equivalent',
   'alpha particle mass energy equivalent in MeV',
   'alpha particle mass in u',
   'electron to alpha particle mass ratio'
]

此方法返回键列表，如果关键字不匹配则返回空。

SciPy - FFTpack (快速傅里叶变换)

傅里叶变换用于计算时域信号，以检查其在频域中的行为。傅里叶变换广泛应用于信号和噪声处理、图像处理、音频信号处理等领域。SciPy 提供了 fftpack 模块，允许用户计算快速傅里叶变换。

下面是一个正弦函数的例子，我们将使用 fftpack 模块计算其傅里叶变换。

快速傅里叶变换

让我们详细了解什么是快速傅里叶变换。

一维离散傅里叶变换

长度为 N 的序列 x[n] 的长度为 N 的 FFT y[k] 由 fft() 计算，逆变换由 ifft() 计算。让我们考虑以下示例

#Importing the fft and inverse fft functions from fftpackage
from scipy.fftpack import fft

#create an array with random n numbers
x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

#Applying the fft function
y = fft(x)
print y

上述程序将生成以下输出。

[ 4.50000000+0.j           2.08155948-1.65109876j   -1.83155948+1.60822041j
 -1.83155948-1.60822041j   2.08155948+1.65109876j ]

让我们看另一个例子

#FFT is already in the workspace, using the same workspace to for inverse transform

yinv = ifft(y)

print yinv

上述程序将生成以下输出。

[ 1.0+0.j   2.0+0.j   1.0+0.j   -1.0+0.j   1.5+0.j ]

scipy.fftpack 模块允许计算快速傅里叶变换。作为一个例子，一个（含噪声的）输入信号可能如下所示：

import numpy as np
time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)
sig = np.sin(2 * np.pi / period * time_vec) + 0.5 *np.random.randn(time_vec.size)
print sig.size

我们正在创建一个时间步长为 0.02 秒的信号。最后一条语句打印信号 sig 的大小。输出如下：

1000

我们不知道信号频率；我们只知道信号 sig 的采样时间步长。该信号应该来自一个实函数，因此傅里叶变换将是对称的。scipy.fftpack.fftfreq() 函数将生成采样频率，而 scipy.fftpack.fft() 将计算快速傅里叶变换。

让我们通过一个例子来理解这一点。

from scipy import fftpack
sample_freq = fftpack.fftfreq(sig.size, d = time_step)
sig_fft = fftpack.fft(sig)
print sig_fft

上述程序将生成以下输出。

array([ 
   25.45122234 +0.00000000e+00j,   6.29800973 +2.20269471e+00j,
   11.52137858 -2.00515732e+01j,   1.08111300 +1.35488579e+01j,
   …….])

离散余弦变换

离散余弦变换 (DCT) 将有限的数据点序列表示为不同频率振荡的余弦函数之和。SciPy 使用函数 dct 提供 DCT，并使用函数 idct 提供相应的 IDCT。让我们考虑以下示例。

from scipy.fftpack import dct
print dct(np.array([4., 3., 5., 10., 5., 3.]))

上述程序将生成以下输出。

array([ 60.,  -3.48476592,  -13.85640646,  11.3137085,  6.,  -6.31319305])

逆离散余弦变换根据其离散余弦变换 (DCT) 系数重建序列。idct 函数是 dct 函数的逆函数。让我们通过以下示例来理解这一点。

from scipy.fftpack import dct
print idct(np.array([4., 3., 5., 10., 5., 3.]))

上述程序将生成以下输出。

array([ 39.15085889, -20.14213562, -6.45392043, 7.13341236,
8.14213562, -3.83035081])

SciPy - Integrate (积分)

当一个函数无法解析地积分，或者解析地积分非常困难时，通常会转向数值积分方法。SciPy 有许多用于执行数值积分的例程。它们大部分都位于相同的 scipy.integrate 库中。下表列出了一些常用的函数。

单重积分

Quad 函数是 SciPy 积分函数的核心。数值积分有时称为求积，因此得名。它通常是执行函数 *f(x)* 在从 a 到 b 的给定固定范围内的单重积分的默认选择。

$$\int_{a}^{b} f(x)dx$$

quad 的一般形式是 scipy.integrate.quad(f, a, b)，其中“f”是要积分的函数的名称。“a”和“b”分别是下限和上限。让我们看一个高斯函数在 0 和 1 范围内的积分示例。

我们首先需要定义函数 → $f(x) = e^{-x^2}$，这可以使用 lambda 表达式完成，然后在该函数上调用 quad 方法。

import scipy.integrate
from numpy import exp
f= lambda x:exp(-x**2)
i = scipy.integrate.quad(f, 0, 1)
print i

上述程序将生成以下输出。

(0.7468241328124271, 8.291413475940725e-15)

quad 函数返回两个值，其中第一个数字是积分值，第二个值是积分值绝对误差的估计值。

注意 - 由于 quad 需要函数作为第一个参数，因此我们不能直接将 exp 作为参数传递。Quad 函数接受正无穷大和负无穷大作为极限。Quad 函数可以积分单个变量的标准预定义 NumPy 函数，例如 exp、sin 和 cos。

多重积分

二重积分和三重积分的机制已被封装到 dblquad、tplquad 和 nquad 函数中。这些函数分别整合四个或六个参数。所有内积分的极限都需要定义为函数。

二重积分

dblquad 的一般形式为 scipy.integrate.dblquad(func, a, b, gfun, hfun)。其中，func 是要积分的函数的名称，“a”和“b”分别是 x 变量的下限和上限，而 gfun 和 hfun 是定义 y 变量的下限和上限的函数的名称。

例如，让我们执行二重积分方法。

$$\int_{0}^{1/2} dy \int_{0}^{\sqrt{1-4y^2}} 16xy \:dx$$

我们使用 lambda 表达式定义函数 f、g 和 h。请注意，即使 g 和 h 是常数，就像在许多情况下一样，它们也必须定义为函数，就像我们这里为下限所做的那样。

import scipy.integrate
from numpy import exp
from math import sqrt
f = lambda x, y : 16*x*y
g = lambda x : 0
h = lambda y : sqrt(1-4*y**2)
i = scipy.integrate.dblquad(f, 0, 0.5, g, h)
print i

上述程序将生成以下输出。

(0.5, 1.7092350012594845e-14)

除了上面描述的例程之外，scipy.integrate 还有一些其他的积分例程，包括 nquad，它执行 n 重多重积分，以及实现各种积分算法的其他例程。但是，quad 和 dblquad 将满足我们对数值积分的大部分需求。

SciPy - Interpolate (插值)

在本章中，我们将讨论插值如何在 SciPy 中发挥作用。

什么是插值？

插值是在线或曲线上的两点之间查找值的過程。为了帮助我们记住它的含义，我们应该将单词的第一部分“inter”理解为“进入”，这提醒我们要查看我们最初拥有的数据的“内部”。插值这个工具不仅在统计学中很有用，而且在科学、商业或需要预测落入两个现有数据点之间的值时也很有用。

让我们创建一些数据，看看如何使用 scipy.interpolate 包进行此插值。

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
x = np.linspace(0, 4, 12)
y = np.cos(x**2/3+4)
print x,y

上述程序将生成以下输出。

(
   array([0.,  0.36363636,  0.72727273,  1.09090909,  1.45454545, 1.81818182, 
          2.18181818,  2.54545455,  2.90909091,  3.27272727,  3.63636364,  4.]),
            
   array([-0.65364362,  -0.61966189,  -0.51077021,  -0.31047698,  -0.00715476,
           0.37976236,   0.76715099,   0.99239518,   0.85886263,   0.27994201,
          -0.52586509,  -0.99582185])
)

现在，我们有两个数组。假设这两个数组是空间中点的两个维度，让我们使用以下程序进行绘图，看看它们是什么样的。

plt.plot(x, y,’o’)
plt.show()

上述程序将生成以下输出。

一维插值

scipy.interpolate 中的 interp1d 类是一种方便的方法，可以根据固定的数据点创建函数，可以使用线性插值在给定数据定义的域内的任何位置进行评估。

使用以上数据，让我们创建一个插值函数并绘制一个新的插值图。

f1 = interp1d(x, y,kind = 'linear')

f2 = interp1d(x, y, kind = 'cubic')

使用 interp1d 函数，我们创建了两个函数 f1 和 f2。这些函数对于给定的输入 x 返回 y。第三个变量 kind 表示插值技术的类型。“线性”、“最近邻”、“零”、“Slinear”、“二次”、“三次”是几种插值技术。

现在，让我们创建一个新的更长输入，以查看插值的明显差异。我们将对新数据使用旧数据的相同函数。

xnew = np.linspace(0, 4,30)

plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')

plt.legend(['data', 'linear', 'cubic','nearest'], loc = 'best')

plt.show()

上述程序将生成以下输出。

样条曲线

为了绘制穿过数据点的平滑曲线，绘图员曾经使用过称为机械样条曲线的薄而灵活的木条、硬橡胶、金属或塑料条。要使用机械样条曲线，需要将大头针放置在设计曲线上的明智选择的点上，然后弯曲样条曲线，使其接触到这些大头针中的每一个。

显然，通过这种构造，样条曲线在这些大头针处对曲线进行插值。它可用于在其他图纸中复制曲线。放置大头针的点称为节点。我们可以通过调整节点的位置来改变样条曲线定义的曲线的形状。

单变量样条曲线

一维平滑样条曲线拟合给定的数据集。scipy.interpolate 中的 UnivariateSpline 类是一种方便的方法，可以根据固定的数据点创建函数类 – scipy.interpolate.UnivariateSpline(x, y, w = None, bbox = [None, None], k = 3, s = None, ext = 0, check_finite = False)。

参数 - 下面是单变量样条曲线的参数。

• 这将 k 次样条曲线 y = spl(x) 拟合到提供的 x、y 数据。

• ‘w’ - 指定样条拟合的权重。必须为正。如果没有（默认值），则所有权重都相等。

• ‘s’ - 通过指定平滑条件来指定节点数。

• ‘k’ - 平滑样条的次数。必须 <= 5。默认为 k = 3，三次样条曲线。

• Ext - 控制不在节点序列定义的区间内的元素的外推模式。

• 如果 ext = 0 或 'extrapolate'，则返回外推值。

• 如果 ext = 1 或 'zero'，则返回 0

• 如果 ext = 2 或 'raise'，则引发 ValueError

• 如果 ext = 3 或 'const'，则返回边界值。

• check_finite – 是否检查输入数组仅包含有限数字。

让我们考虑以下示例。

import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
x = np.linspace(-3, 3, 50)
y = np.exp(-x**2) + 0.1 * np.random.randn(50)
plt.plot(x, y, 'ro', ms = 5)
plt.show()

使用平滑参数的默认值。

spl = UnivariateSpline(x, y)
xs = np.linspace(-3, 3, 1000)
plt.plot(xs, spl(xs), 'g', lw = 3)
plt.show()

手动更改平滑量。

spl.set_smoothing_factor(0.5)
plt.plot(xs, spl(xs), 'b', lw = 3)
plt.show()

SciPy - 输入和输出

Scipy.io（输入和输出）包提供了广泛的功能，可以处理不同格式的文件。其中一些格式是：

• MATLAB
• IDL
• 矩阵市场
• Wave
• Arff
• Netcdf 等。

让我们详细讨论最常用的文件格式：

MATLAB

以下是用于加载和保存 .mat 文件的函数。

让我们考虑以下示例。

import scipy.io as sio
import numpy as np

#Save a mat file
vect = np.arange(10)
sio.savemat('array.mat', {'vect':vect})

#Now Load the File
mat_file_content = sio.loadmat(‘array.mat’)
Print mat_file_content

上述程序将生成以下输出。

{
   'vect': array([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]]), '__version__': '1.0', 
   '__header__': 'MATLAB 5.0 MAT-file Platform: posix, Created on: Sat Sep 30 
   09:49:32 2017', '__globals__': []
}

我们可以看到包含元信息的数组。如果我们想在不将数据读入内存的情况下检查 MATLAB 文件的内容，请使用如下所示的 whosmat 命令。

import scipy.io as sio
mat_file_content = sio.whosmat(‘array.mat’)
print mat_file_content

上述程序将生成以下输出。

[('vect', (1, 10), 'int64')]

SciPy - Linalg (线性代数)

SciPy 使用优化的 ATLAS LAPACK 和 BLAS 库构建。它具有非常快速的线性代数能力。所有这些线性代数例程都期望一个可以转换为二维数组的对象。这些例程的输出也是二维数组。

SciPy.linalg 与 NumPy.linalg 的比较

scipy.linalg 包含 numpy.linalg 中的所有函数。此外，scipy.linalg 还有一些 numpy.linalg 中没有的其他高级函数。与 numpy.linalg 相比，使用 scipy.linalg 的另一个优点是它始终与 BLAS/LAPACK 支持一起编译，而对于 NumPy，这是可选的。因此，SciPy 版本的速度可能更快，具体取决于 NumPy 的安装方式。

线性方程

scipy.linalg.solve 功能求解线性方程 a * x + b * y = Z，以求解未知的 x、y 值。

例如，假设需要求解以下联立方程。

x + 3y + 5z = 10

2x + 5y + z = 8

2x + 3y + 8z = 3

为了求解上述方程组中的 x、y、z 值，我们可以使用矩阵求逆法找到解向量，如下所示。

$$

但是，最好使用linalg.solve命令，因为它速度更快，数值稳定性也更好。

solve 函数接受两个输入 ‘a’ 和 ‘b’，其中 ‘a’ 表示系数，‘b’ 表示右侧的值，并返回解数组。

让我们考虑以下示例。

#importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy arrays
a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])
b = np.array([2, 4, -1])

#Passing the values to the solve function
x = linalg.solve(a, b)

#printing the result array
print x

上述程序将生成以下输出。

array([ 2., -2., 9.])

求解行列式

方阵 A 的行列式通常表示为 |A|，它是线性代数中经常使用的量。在 SciPy 中，可以使用det()函数计算它。它接收一个矩阵作为输入，并返回一个标量值。

让我们考虑以下示例。

#importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy array
A = np.array([[1,2],[3,4]])

#Passing the values to the det function
x = linalg.det(A)

#printing the result
print x

上述程序将生成以下输出。

-2.0

特征值和特征向量

特征值-特征向量问题是最常用的线性代数运算之一。我们可以通过考虑以下关系式来找到方阵 (A) 的特征值 (λ) 和相应的特征向量 (v)：

Av = λv

scipy.linalg.eig计算普通或广义特征值问题的特征值。此函数返回特征值和特征向量。

让我们考虑以下示例。

#importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy array
A = np.array([[1,2],[3,4]])

#Passing the values to the eig function
l, v = linalg.eig(A)

#printing the result for eigen values
print l

#printing the result for eigen vectors
print v

上述程序将生成以下输出。

array([-0.37228132+0.j, 5.37228132+0.j]) #--Eigen Values
array([[-0.82456484, -0.41597356], #--Eigen Vectors
       [ 0.56576746, -0.90937671]])

奇异值分解

奇异值分解 (SVD) 可以被认为是特征值问题到非方阵矩阵的扩展。

scipy.linalg.svd将矩阵 ‘a’ 分解为两个酉矩阵 ‘U’ 和 ‘Vh’ 以及一个包含奇异值（实数，非负）的一维数组 ‘s’，使得 a == U*S*Vh，其中 ‘S’ 是一个形状合适的矩阵，其主对角线为 ‘s’，其余元素为零。

让我们考虑以下示例。

#importing the scipy and numpy packages
from scipy import linalg
import numpy as np

#Declaring the numpy array
a = np.random.randn(3, 2) + 1.j*np.random.randn(3, 2)

#Passing the values to the eig function
U, s, Vh = linalg.svd(a)

# printing the result
print U, Vh, s

上述程序将生成以下输出。

(
   array([
      [ 0.54828424-0.23329795j, -0.38465728+0.01566714j,
      -0.18764355+0.67936712j],
      [-0.27123194-0.5327436j , -0.57080163-0.00266155j,
      -0.39868941-0.39729416j],
      [ 0.34443818+0.4110186j , -0.47972716+0.54390586j,
      0.25028608-0.35186815j]
   ]),

   array([ 3.25745379, 1.16150607]),

   array([
      [-0.35312444+0.j , 0.32400401+0.87768134j],
      [-0.93557636+0.j , -0.12229224-0.33127251j]
   ])
)

SciPy - Ndimage (多维图像)

SciPy 的 ndimage 子模块专门用于图像处理。这里，ndimage 指的是 n 维图像。

图像处理中最常见的一些任务如下：

• 输入/输出，显示图像
• 基本操作——裁剪、翻转、旋转等。
• 图像滤波——去噪、锐化等。
• 图像分割——标记对应于不同对象的像素
• 分类
• 特征提取
• 配准

让我们讨论如何使用 SciPy 实现其中一些任务。

打开和写入图像文件

SciPy 中的misc 包包含一些图像。我们使用这些图像来学习图像处理操作。让我们考虑以下示例。

from scipy import misc
f = misc.face()
misc.imsave('face.png', f) # uses the Image module (PIL)

import matplotlib.pyplot as plt
plt.imshow(f)
plt.show()

上述程序将生成以下输出。

任何原始格式的图像都是矩阵格式中数字表示的颜色的组合。机器仅根据这些数字来理解和处理图像。RGB 是一种流行的表示方式。

让我们看看上面图像的统计信息。

from scipy import misc
face = misc.face(gray = False)
print face.mean(), face.max(), face.min()

上述程序将生成以下输出。

110.16274388631184, 255, 0

现在，我们知道图像是由数字构成的，因此数字值的任何变化都会改变原始图像。让我们对图像执行一些几何变换。基本几何操作是裁剪。

from scipy import misc
face = misc.face(gray = True)
lx, ly = face.shape
# Cropping
crop_face = face[lx / 4: - lx / 4, ly / 4: - ly / 4]
import matplotlib.pyplot as plt
plt.imshow(crop_face)
plt.show()

上述程序将生成以下输出。

我们还可以执行一些基本操作，例如将图像倒置，如下所示。

# up <-> down flip
from scipy import misc
face = misc.face()
flip_ud_face = np.flipud(face)

import matplotlib.pyplot as plt
plt.imshow(flip_ud_face)
plt.show()

上述程序将生成以下输出。

除此之外，我们还有rotate()函数，它可以将图像旋转指定的角度。

# rotation
from scipy import misc,ndimage
face = misc.face()
rotate_face = ndimage.rotate(face, 45)

import matplotlib.pyplot as plt
plt.imshow(rotate_face)
plt.show()

上述程序将生成以下输出。

滤波器

让我们讨论滤波器如何帮助图像处理。

什么是图像处理中的滤波？

滤波是一种修改或增强图像的技术。例如，您可以滤波图像以强调某些特征或去除其他特征。使用滤波实现的图像处理操作包括平滑、锐化和边缘增强。

滤波是一种邻域运算，其中输出图像中任何给定像素的值是通过对相应输入像素邻域中的像素值应用某种算法来确定的。现在让我们使用 SciPy ndimage 执行一些操作。

模糊

模糊广泛用于减少图像中的噪声。我们可以执行滤波操作并查看图像的变化。让我们考虑以下示例。

from scipy import misc
face = misc.face()
blurred_face = ndimage.gaussian_filter(face, sigma=3)
import matplotlib.pyplot as plt
plt.imshow(blurred_face)
plt.show()

上述程序将生成以下输出。

sigma 值表示五级模糊程度。我们可以通过调整 sigma 值来查看图像质量的变化。有关模糊的更多详细信息，请点击 → DIP（数字图像处理）教程。

边缘检测

让我们讨论边缘检测如何帮助图像处理。

什么是边缘检测？

边缘检测是一种图像处理技术，用于查找图像中对象的边界。它通过检测亮度的不连续性来工作。边缘检测用于图像分割和数据提取，应用领域包括图像处理、计算机视觉和机器视觉。

最常用的边缘检测算法包括：

• Sobel
• Canny
• Prewitt
• Roberts
• 模糊逻辑方法

让我们考虑以下示例。

import scipy.ndimage as nd
import numpy as np

im = np.zeros((256, 256))
im[64:-64, 64:-64] = 1
im[90:-90,90:-90] = 2
im = ndimage.gaussian_filter(im, 8)

import matplotlib.pyplot as plt
plt.imshow(im)
plt.show()

上述程序将生成以下输出。

图像看起来像一个彩色的正方形块。现在，我们将检测这些彩色块的边缘。这里，ndimage 提供了一个名为Sobel的函数来执行此操作。而 NumPy 提供Hypot函数将两个结果矩阵组合成一个。

让我们考虑以下示例。

import scipy.ndimage as nd
import matplotlib.pyplot as plt

im = np.zeros((256, 256))
im[64:-64, 64:-64] = 1
im[90:-90,90:-90] = 2
im = ndimage.gaussian_filter(im, 8)

sx = ndimage.sobel(im, axis = 0, mode = 'constant')
sy = ndimage.sobel(im, axis = 1, mode = 'constant')
sob = np.hypot(sx, sy)

plt.imshow(sob)
plt.show()

上述程序将生成以下输出。

SciPy - Optimize (优化)

scipy.optimize 包提供了几个常用的优化算法。该模块包含以下方面：

• 使用各种算法（例如 BFGS、Nelder-Mead 单纯形、牛顿共轭梯度、COBYLA 或 SLSQP）对多元标量函数进行无约束和约束最小化 (minimize())

• 全局（蛮力）优化例程（例如，anneal()、basinhopping()）

• 最小二乘法最小化 (leastsq()) 和曲线拟合 (curve_fit()) 算法

• 标量单变量函数最小化器 (minimize_scalar()) 和求根器 (newton())

• 使用各种算法（例如混合 Powell、Levenberg-Marquardt 或大规模方法，如 Newton-Krylov）的多元方程组求解器 (root())

多元标量函数的无约束和约束最小化

minimize() 函数为scipy.optimize中多元标量函数的无约束和约束最小化算法提供了一个公共接口。为了演示最小化函数，考虑最小化 NN 变量的 Rosenbrock 函数的问题：

$$f(x) = \sum_{i = 1}^{N-1} \:100(x_i - x_{i-1}^{2})$$

此函数的最小值为 0，当 xi = 1 时达到。

Nelder-Mead 单纯形算法

在以下示例中，minimize() 例程与Nelder-Mead 单纯形算法 (method = 'Nelder-Mead') 一起使用（通过 method 参数选择）。让我们考虑以下示例。

import numpy as np
from scipy.optimize import minimize

def rosen(x):

x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead')

print(res.x)

上述程序将生成以下输出。

[7.93700741e+54  -5.41692163e+53  6.28769150e+53  1.38050484e+55  -4.14751333e+54]

单纯形算法可能是最小化相当良好的函数的最简单方法。它只需要函数评估，并且是简单最小化问题的不错选择。但是，因为它不使用任何梯度评估，所以找到最小值可能需要更长的时间。

另一个只需要函数调用就能找到最小值的优化算法是Powell 方法，可以通过在 minimize() 函数中设置 method = 'powell' 来使用。

最小二乘法

求解具有变量边界的非线性最小二乘问题。给定残差 f(x)（n 个实变量的 m 维实函数）和损失函数 rho(s)（标量函数），least_squares 找到代价函数 F(x) 的局部最小值。让我们考虑以下示例。

在这个例子中，我们找到了没有自变量边界的 Rosenbrock 函数的最小值。

#Rosenbrock Function
def fun_rosenbrock(x):
   return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])
   
from scipy.optimize import least_squares
input = np.array([2, 2])
res = least_squares(fun_rosenbrock, input)

print res

请注意，我们只提供了残差向量。该算法将代价函数构造为残差平方和，这给出了 Rosenbrock 函数。精确的最小值在 x = [1.0,1.0] 处。

上述程序将生成以下输出。

active_mask: array([ 0., 0.])
      cost: 9.8669242910846867e-30
      fun: array([ 4.44089210e-15, 1.11022302e-16])
      grad: array([ -8.89288649e-14, 4.44089210e-14])
      jac: array([[-20.00000015,10.],[ -1.,0.]])
   message: '`gtol` termination condition is satisfied.'
      nfev: 3
      njev: 3
   optimality: 8.8928864934219529e-14
      status: 1
      success: True
         x: array([ 1., 1.])

求根

让我们了解求根如何在 SciPy 中发挥作用。

标量函数

如果只有一个变量方程，则可以尝试四种不同的求根算法。这些算法中的每一个都需要一个区间端点，在这个区间中预期存在根（因为函数符号发生变化）。通常，brentq 是最佳选择，但在某些情况下或出于学术目的，其他方法可能有用。

不动点求解

与寻找函数的零点密切相关的问题是寻找函数的不动点的问题。函数的不动点是函数求值返回该点的点：g(x) = x。显然gg的不动点是f(x) = g(x)-x的根。等效地，ff的根是g(x) = f(x)+x的不动点。例程fixed_point提供了一种简单的迭代方法，使用Aitkens序列加速来估计gg的不动点，如果给定一个起始点。

方程组

可以使用root()函数找到一组非线性方程的根。可以使用多种方法，其中hybr（默认）和 lm 分别使用来自 MINPACK 的Powell 混合方法和Levenberg-Marquardt 方法。

以下示例考虑单变量超越方程。

x² + 2cos(x) = 0

其根可以如下找到：

import numpy as np
from scipy.optimize import root
def func(x):
   return x*2 + 2 * np.cos(x)
sol = root(func, 0.3)
print sol

上述程序将生成以下输出。

fjac: array([[-1.]])
fun: array([ 2.22044605e-16])
message: 'The solution converged.'
   nfev: 10
   qtf: array([ -2.77644574e-12])
      r: array([-3.34722409])
   status: 1
   success: True
      x: array([-0.73908513])

SciPy - Stats (统计)

所有统计函数都位于子包scipy.stats中，可以使用info(stats)函数获得这些函数的相当完整的列表。也可以从 stats 子包的docstring中获得可用随机变量的列表。该模块包含大量概率分布以及不断增长的统计函数库。

每个单变量分布都有自己的子类，如下表所述：

正态连续随机变量

一个概率分布，其中随机变量 X 可以取任何值是连续随机变量。location (loc) 关键字指定均值。scale (scale) 关键字指定标准差。

作为rv_continuous类的实例，norm对象继承了它的通用方法集合，并用特定于此特定分布的细节来完成它们。

为了计算多个点的 CDF，我们可以传递一个列表或一个 NumPy 数组。让我们考虑以下示例。

from scipy.stats import norm
import numpy as np
print norm.cdf(np.array([1,-1., 0, 1, 3, 4, -2, 6]))

上述程序将生成以下输出。

array([ 0.84134475, 0.15865525, 0.5 , 0.84134475, 0.9986501 ,
0.99996833, 0.02275013, 1. ])

为了找到分布的中位数，我们可以使用百分点函数 (PPF)，它是 CDF 的反函数。让我们通过以下示例来理解。

from scipy.stats import norm
print norm.ppf(0.5)

上述程序将生成以下输出。

0.0

为了生成一系列随机变量，我们应该使用 size 关键字参数，这在以下示例中显示。

from scipy.stats import norm
print norm.rvs(size = 5)

上述程序将生成以下输出。

array([ 0.20929928, -1.91049255, 0.41264672, -0.7135557 , -0.03833048])

以上输出不可重现。要生成相同的随机数，请使用 seed 函数。

均匀分布

可以使用 uniform 函数生成均匀分布。让我们考虑以下示例。

from scipy.stats import uniform
print uniform.cdf([0, 1, 2, 3, 4, 5], loc = 1, scale = 4)

上述程序将生成以下输出。

array([ 0. , 0. , 0.25, 0.5 , 0.75, 1. ])

构建离散分布

让我们生成一个随机样本，并将观察到的频率与概率进行比较。

二项分布

作为rv_discrete 类的实例，binom 对象继承了它的通用方法集合，并用特定于此特定分布的细节来完成它们。让我们考虑以下示例。

from scipy.stats import uniform
print uniform.cdf([0, 1, 2, 3, 4, 5], loc = 1, scale = 4)

上述程序将生成以下输出。

array([ 0. , 0. , 0.25, 0.5 , 0.75, 1. ])

描述性统计

基本统计数据，例如最小值、最大值、平均值和方差，以NumPy数组作为输入并返回相应的结果。scipy.stats包中提供的一些基本统计函数在下面的表格中进行了描述。

其中许多函数在scipy.stats.mstats中都有类似的版本，这些版本适用于掩码数组。让我们通过下面的示例来了解这一点。

from scipy import stats
import numpy as np
x = np.array([1,2,3,4,5,6,7,8,9])
print x.max(),x.min(),x.mean(),x.var()

上述程序将生成以下输出。

(9, 1, 5.0, 6.666666666666667)

T检验

让我们了解T检验在SciPy中的用途。

ttest_1samp

计算一组分数的均值的T检验。这是一个针对零假设的双尾检验，该零假设是独立观察样本‘a’的期望值（均值）等于给定的总体均值popmean。让我们考虑以下示例。

from scipy import stats
rvs = stats.norm.rvs(loc = 5, scale = 10, size = (50,2))
print stats.ttest_1samp(rvs,5.0)

上述程序将生成以下输出。

Ttest_1sampResult(statistic = array([-1.40184894, 2.70158009]),
pvalue = array([ 0.16726344, 0.00945234]))

比较两个样本

在以下示例中，有两个样本，它们可以来自相同的分布或不同的分布，我们想要检验这些样本是否具有相同的统计特性。

ttest_ind - 计算两个独立样本分数均值的T检验。这是一个针对零假设的双尾检验，该零假设是两个独立样本具有相同的平均值（期望值）。默认情况下，此检验假设总体具有相同的方差。

如果我们观察到来自相同或不同总体的两个独立样本，我们可以使用此检验。让我们考虑以下示例。

from scipy import stats
rvs1 = stats.norm.rvs(loc = 5,scale = 10,size = 500)
rvs2 = stats.norm.rvs(loc = 5,scale = 10,size = 500)
print stats.ttest_ind(rvs1,rvs2)

上述程序将生成以下输出。

Ttest_indResult(statistic = -0.67406312233650278, pvalue = 0.50042727502272966)

您可以使用相同长度但均值不同的新数组来测试相同内容。在loc中使用不同的值并进行相同的测试。

SciPy - CSGraph (压缩稀疏图)

CSGraph代表压缩稀疏图，它专注于基于稀疏矩阵表示的快速图算法。

图表示

首先，让我们了解什么是稀疏图以及它如何帮助图表示。

什么是稀疏图？

图只是一组节点，它们之间有链接。图几乎可以表示任何事物——社交网络连接（每个节点都是一个人，并连接到熟人）；图像（每个节点都是一个像素，并连接到相邻像素）；高维分布中的点（每个节点都连接到其最近邻）；以及你能想象到的几乎任何其他事物。

表示图数据的一种非常有效的方法是使用稀疏矩阵：我们称之为G。矩阵G的大小为N x N，G[i, j]给出节点'i'和节点'j'之间连接的值。稀疏图主要包含零——也就是说，大多数节点只有很少的连接。在大多数感兴趣的情况下，此属性都是正确的。

创建稀疏图子模块的动机是scikit-learn中使用的几种算法，其中包括以下算法：

• Isomap - 一种流形学习算法，需要在图中找到最短路径。

• 层次聚类 - 基于最小生成树的聚类算法。

• 谱分解 - 基于稀疏图拉普拉斯算子的投影算法。

作为一个具体的例子，假设我们想表示以下无向图：

该图具有三个节点，其中节点0和1由权重为2的边连接，节点0和2由权重为1的边连接。我们可以构建密集的、掩码的和稀疏的表示，如下例所示，记住无向图由对称矩阵表示。

G_dense = np.array([ [0, 2, 1],
                     [2, 0, 0],
                     [1, 0, 0] ])
                     
G_masked = np.ma.masked_values(G_dense, 0)
from scipy.sparse import csr_matrix

G_sparse = csr_matrix(G_dense)
print G_sparse.data

上述程序将生成以下输出。

array([2, 1, 2, 1])

这与之前的图相同，只是节点0和2由权重为零的边连接。在这种情况下，上面的密集表示会导致歧义——如果零是一个有意义的值，那么如何表示非边？在这种情况下，必须使用掩码表示或稀疏表示来消除歧义。

让我们考虑以下示例。

from scipy.sparse.csgraph import csgraph_from_dense
G2_data = np.array
([
   [np.inf, 2, 0 ],
   [2, np.inf, np.inf],
   [0, np.inf, np.inf]
])
G2_sparse = csgraph_from_dense(G2_data, null_value=np.inf)
print G2_sparse.data

上述程序将生成以下输出。

array([ 2., 0., 2., 0.])

使用稀疏图的单词接龙

单词接龙是由刘易斯·卡罗尔发明的一个游戏，在这个游戏中，单词通过每次更改一个字母来链接。例如：

APE → APT → AIT → BIT → BIG → BAG → MAG → MAN

在这里，我们用七步从“APE”到“MAN”，每次更改一个字母。问题是——我们能否使用相同的规则找到这两个单词之间更短的路径？这个问题很自然地表达为一个稀疏图问题。节点将对应于单个单词，我们将创建连接在最多相差一个字母的单词对之间的连接。

获取单词列表

首先，我们必须获得有效单词的列表。我正在运行Mac，Mac在以下代码块中给出的位置有一个词典。如果您使用的是不同的架构，您可能需要搜索一下才能找到您的系统词典。

wordlist = open('/usr/share/dict/words').read().split()
print len(wordlist)

上述程序将生成以下输出。

235886

现在我们想查看长度为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 = map(str.lower, word_list)
print len(word_list)

上述程序将生成以下输出。

1135

现在，我们有1135个有效的三个字母单词的列表（确切的数字可能会根据使用的特定列表而变化）。这些单词中的每一个都将成为我们图中的一个节点，我们将创建连接与仅相差一个字母的单词对相关的节点的边。

import numpy as np
word_list = np.asarray(word_list)

word_list.dtype
word_list.sort()

word_bytes = np.ndarray((word_list.size, word_list.itemsize),
   dtype = 'int8',
   buffer = word_list.data)
print word_bytes.shape

上述程序将生成以下输出。

(1135, 3)

我们将使用每个点之间的汉明距离来确定哪些单词对是连接的。汉明距离衡量两个向量之间不同条目的分数：任何两个汉明距离等于1/N的单词，其中N是单词中字母的数量，在单词接龙中是连接的。

from scipy.spatial.distance import pdist, squareform
from scipy.sparse import csr_matrix
hamming_dist = pdist(word_bytes, metric = 'hamming')
graph = csr_matrix(squareform(hamming_dist < 1.5 / word_list.itemsize))

在比较距离时，我们不使用等式，因为这对于浮点值来说可能不稳定。只要单词列表中没有两个条目相同，不等式就会产生所需的结果。现在，我们的图已经设置好了，我们将使用最短路径搜索来查找图中任意两个单词之间的路径。

i1 = word_list.searchsorted('ape')
i2 = word_list.searchsorted('man')
print word_list[i1],word_list[i2]

上述程序将生成以下输出。

ape, man

我们需要检查这些是否匹配，因为如果这些单词不在列表中，输出将出现错误。现在，我们只需要找到图中这两个索引之间的最短路径。我们将使用迪杰斯特拉算法，因为它允许我们仅查找一个节点的路径。

from scipy.sparse.csgraph import dijkstra
distances, predecessors = dijkstra(graph, indices = i1, return_predecessors = True)
print distances[i2]

上述程序将生成以下输出。

5.0

因此，我们看到“ape”和“man”之间的最短路径仅包含五个步骤。我们可以使用算法返回的前驱来重建此路径。

path = []
i = i2

while i != i1:
   path.append(word_list[i])
   i = predecessors[i]
   
path.append(word_list[i1])
print path[::-1]i2]

上述程序将生成以下输出。

['ape', 'ope', 'opt', 'oat', 'mat', 'man']

SciPy - Spatial (空间数据)

scipy.spatial包可以通过利用Qhull库来计算一组点的三角剖分、沃罗诺伊图和凸包。此外，它还包含用于最近邻点查询的KDTree实现以及用于各种度量中距离计算的实用程序。

德劳内三角剖分

让我们了解什么是德劳内三角剖分以及它们如何在SciPy中使用。

什么是德劳内三角剖分？

在数学和计算几何中，给定平面中离散点集P的德劳内三角剖分是三角剖分DT(P)，这样P中没有点位于DT(P)中任何三角形的圆周内。

我们可以通过SciPy计算相同的内容。让我们考虑以下示例。

from scipy.spatial import Delaunay
points = np.array([[0, 4], [2, 1.1], [1, 3], [1, 2]])
tri = Delaunay(points)
import matplotlib.pyplot as plt
plt.triplot(points[:,0], points[:,1], tri.simplices.copy())
plt.plot(points[:,0], points[:,1], 'o')
plt.show()

上述程序将生成以下输出。

共面点

让我们了解什么是共面点以及它们如何在SciPy中使用。

什么是共面点？

共面点是位于同一平面上的三个或更多个点。回想一下，平面是一个平面，它向各个方向无限延伸。它通常在数学教科书中显示为一个四面体。

让我们看看如何使用SciPy找到它。让我们考虑以下示例。

from scipy.spatial import Delaunay
points = np.array([[0, 0], [0, 1], [1, 0], [1, 1], [1, 1]])
tri = Delaunay(points)
print tri.coplanar

上述程序将生成以下输出。

array([[4, 0, 3]], dtype = int32)

这意味着点4位于三角形0和顶点3附近，但不包含在三角剖分中。

凸包

让我们了解什么是凸包以及它们如何在SciPy中使用。

什么是凸包？

在数学中，欧几里得平面或欧几里得空间（或更一般地，实数上的仿射空间）中的一组点X的凸包或凸包络是最小的凸集，它包含X。

让我们考虑以下示例以详细了解它。

from scipy.spatial import ConvexHull
points = np.random.rand(10, 2) # 30 random points in 2-D
hull = ConvexHull(points)
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 - ODR (正交距离回归)

ODR代表正交距离回归，用于回归研究。基本线性回归通常用于通过在图上绘制最佳拟合线来估计两个变量y和x之间的关系。

为此使用的一种数学方法称为最小二乘法，其目的是最小化每个点的平方误差之和。这里关键的问题是如何计算每个点的误差（也称为残差）？

在标准线性回归中，目标是根据X值预测Y值——因此，明智的做法是计算Y值的误差（如下图中的灰色线所示）。但是，有时更明智的做法是同时考虑X和Y的误差（如下图中的红色虚线所示）。

例如：当您知道您的X测量值不确定时，或者当您不想关注一个变量的误差超过另一个变量时。

正交距离回归（ODR）是一种可以做到这一点的方法（在这种情况下，正交意味着垂直——因此它计算垂直于线的误差，而不仅仅是“垂直”）。

scipy.odr对单变量回归的实现

以下示例演示了scipy.odr对单变量回归的实现。

import numpy as np
import matplotlib.pyplot as plt
from scipy.odr import *
import random

# Initiate some data, giving some randomness using random.random().
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([i**2 + random.random() for i in x])

# Define a function (quadratic in our case) to fit the data with.
def linear_func(p, x):
   m, c = p
   return m*x + c

# Create a model for fitting.
linear_model = Model(linear_func)

# Create a RealData object using our initiated data from above.
data = RealData(x, y)

# Set up ODR with the model and data.
odr = ODR(data, linear_model, beta0=[0., 1.])

# Run the regression.
out = odr.run()

# Use the in-built pprint method to give us results.
out.pprint()

上述程序将生成以下输出。

Beta: [ 5.51846098 -4.25744878]
Beta Std Error: [ 0.7786442 2.33126407]

Beta Covariance: [
   [ 1.93150969 -4.82877433]
   [ -4.82877433 17.31417201
]]

Residual Variance: 0.313892697582
Inverse Condition #: 0.146618499389
Reason(s) for Halting:
   Sum of squares convergence

SciPy - Special Package (特殊函数包)

special包中提供的函数是通用函数，它们遵循广播和自动数组循环。

让我们看看一些最常用的特殊函数：

• 三次根函数
• 指数函数
• 相对误差指数函数
• 对数和指数函数
• 朗伯函数
• 排列和组合函数
• 伽马函数

现在让我们简要了解一下这些函数。

三次根函数

此三次方根函数的语法为 – scipy.special.cbrt(x)。这将获取x 的逐元素立方根。

让我们考虑以下示例。

from scipy.special import cbrt
res = cbrt([10, 9, 0.1254, 234])
print res

上述程序将生成以下输出。

[ 2.15443469 2.08008382 0.50053277 6.16224015]

指数函数

指数函数的语法为 – scipy.special.exp10(x)。这将计算 10**x 的逐元素值。

让我们考虑以下示例。

from scipy.special import exp10
res = exp10([2, 9])
print res

上述程序将生成以下输出。

[1.00000000e+02  1.00000000e+09]

相对误差指数函数

此函数的语法为 – scipy.special.exprel(x)。它生成相对误差指数 (exp(x) - 1)/x。

当x 接近零时，exp(x) 接近 1，因此 exp(x) - 1 的数值计算可能会遭受灾难性的精度损失。然后实现 exprel(x) 来避免当x 接近零时发生的精度损失。

让我们考虑以下示例。

from scipy.special import exprel
res = exprel([-0.25, -0.1, 0, 0.1, 0.25])
print res

上述程序将生成以下输出。

[0.88479687 0.95162582 1.   1.05170918 1.13610167]

对数和指数函数

此函数的语法为 – scipy.special.logsumexp(x)。它有助于计算输入元素指数和的对数。

让我们考虑以下示例。

from scipy.special import logsumexp
import numpy as np
a = np.arange(10)
res = logsumexp(a)
print res

上述程序将生成以下输出。

9.45862974443

朗伯函数

此函数的语法为 – scipy.special.lambertw(x)。它也称为兰伯特 W 函数。兰伯特 W 函数 W(z) 定义为 w * exp(w) 的反函数。换句话说，W(z) 的值使得对于任何复数 z，z = W(z) * exp(W(z))。

兰伯特 W 函数是一个多值函数，具有无限多个分支。每个分支给出了方程 z = w exp(w) 的一个单独解。这里，分支由整数 k 索引。

让我们考虑以下示例。这里，兰伯特 W 函数是 w exp(w) 的反函数。

from scipy.special import lambertw
w = lambertw(1)
print w
print w * np.exp(w)

上述程序将生成以下输出。

(0.56714329041+0j)
(1+0j)

排列与组合

让我们分别讨论排列和组合以便更清晰地理解它们。

组合 − 组合函数的语法为 – scipy.special.comb(N,k)。让我们考虑以下示例 −

from scipy.special import comb
res = comb(10, 3, exact = False,repetition=True)
print res

上述程序将生成以下输出。

220.0

注意 − 仅当 exact = False 时才接受数组参数。如果 k > N，N < 0 或 k < 0，则返回 0。

排列 − 组合函数的语法为 – scipy.special.perm(N,k)。从 N 个事物中取 k 个事物的排列，即 N 的 k 排列。这也称为“部分排列”。

让我们考虑以下示例。

from scipy.special import perm
res = perm(10, 3, exact = True)
print res

上述程序将生成以下输出。

720

伽马函数

伽马函数通常被称为广义阶乘，因为 z*gamma(z) = gamma(z+1)，对于自然数 'n'，gamma(n+1) = n!。

组合函数的语法为 – scipy.special.gamma(x)。从 N 个事物中取 k 个事物的排列，即 N 的 k 排列。这也称为“部分排列”。

组合函数的语法为 – scipy.special.gamma(x)。从 N 个事物中取 k 个事物的排列，即 N 的 k 排列。这也称为“部分排列”。

from scipy.special import gamma
res = gamma([0, 0.5, 1, 5])
print res

上述程序将生成以下输出。

[inf  1.77245385  1.  24.]