SciPy - 快速指南


SciPy - 简介

SciPy,发音为 Sigh Pi,是一个科学 Python 开源软件,在 BSD 许可库下分发,用于执行数学、科学和工程计算。

SciPy 库依赖于 NumPy,它提供方便快捷的 N 维数组操作。SciPy 库专为与 NumPy 数组配合使用而构建,并提供许多用户友好且高效的数值实践,例如数值积分和优化例程。它们一起运行在所有流行的操作系统上,安装快速并且免费。NumPy 和 SciPy 易于使用,但功能强大,足以受到一些世界领先的科学家和工程师的依赖。

SciPy 子包

SciPy 被组织成涵盖不同科学计算领域的子包。下表总结了这些 -

scipy.cluster 矢量量化/Kmeans
scipy.常量 物理和数学常数
scipy.fftpack 傅里叶变换
scipy.integrate 集成例程
scipy.插值 插值法
scipy.io 数据输入输出
scipy.linalg 线性代数例程
scipy.ndimage n维图像包
scipy.odr 正交距离回归
scipy.优化 优化
scipy.信号 信号处理
scipy.稀疏 稀疏矩阵
scipy.spatial 空间数据结构和算法
scipy.special 任何特殊的数学函数
scipy.stats 统计数据

数据结构

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) - 它是一个免费的 Python 发行版,带有适用于 Windows 操作系统的 SciPy 堆栈和 Spyder IDE。(可从https://python-xy.github.io/下载)

Linux

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

乌班图

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

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

软呢帽

我们可以使用以下路径在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(shape) 函数将创建一个填充有 0 值且具有指定形状的数组。默认数据类型是 float64。让我们考虑下面的例子。

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

上述程序的输出如下。

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

使用一个()

Ones(shape) 函数将创建一个包含 1 个值的数组。它在所有其他方面与零相同。让我们考虑下面的例子。

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-Means 算法的良好实现。让我们了解如何使用它。

SciPy 中的 K 均值实现

我们将了解如何在 SciPy 中实现 K-Means。

导入 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-Means 之前,通过白化重新调整观察集的每个特征维度是有益的。每个特征除以所有观测值的标准差,得到单位方差。

数据白化

我们必须使用以下代码来美白数据。

# whitening of data
data = whiten(data)

使用三个集群计算 K 均值

现在让我们使用以下代码计算三个集群的 K 均值。

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

上面的代码对形成 K 个簇的一组观察向量执行 K-Means。K-Means 算法调整质心,直到无法取得足够的进展,即失真的变化,因为最后一次迭代小于某个阈值。在这里,我们可以通过使用下面给出的代码打印 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”ד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

可用常量列表

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

数学常数

先生。没有。 持续的 描述
1 圆周率 圆周率
2 金的 黄金比例

物理常数

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

先生。没有。 常数及说明
1

C

真空中的光速

2

光速

真空中的光速

3

H

普朗克常数

4

普朗克

普朗克常数 h

5

G

牛顿引力常数

6

e

基本电荷

7

气体摩尔常数

8

阿伏加德罗

阿伏伽德罗常数

9

k

玻尔兹曼常数

10

电子质量(OR) m_e

电子质量

11

质子质量 (OR) m_p

质子质量

12

中子质量(OR)m_n

中子质量

单位

下表列出了 SI 单位。

先生。没有。 单元 价值
1 0.001
2 1E-06
3 公斤 1000

这些单位范围从yotta、zetta、exa、peta、tera……kilo、hector、……nano、pico……到zepto。

其他重要常数

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

先生。没有。 单元 价值
1 公克 0.001公斤
2 Atomics质量 Atomics质量常数
3 程度 以弧度为单位的度数
4 分钟 一分钟几秒
5 一天一秒
6 英寸 一英寸等于米
7 微米 一微米换算为米
8 光年 一光年(米)
9 自动取款机 标准大气压,单位:帕斯卡
10 英亩 一英亩等于平方米
11 一升换算成立方米
12 加仑 一加仑换算成立方米
13 公里小时 公里每小时(米每秒)
14 度_华氏度 一华氏度开尔文
15 电子伏特 一电子伏特(焦耳)
16 生命值 一马力,单位为瓦
17 号 动态 一达因(牛顿)
18 λ2nu 将波长转换为光频率

记住所有这些有点困难。获取哪个键对应哪个函数的简单方法是使用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

对时域信号进行傅里叶变换计算,以检查其在频域中的Behave。傅里叶变换在信号和噪声处理、图像处理、音频信号处理等学科中得到应用。SciPy 提供了 fftpack 模块,可让用户计算快速傅里叶变换。

以下是正弦函数的示例,它将用于使用 fftpack 模块计算傅立叶变换。

快速傅立叶变换

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

一维离散傅立叶变换

通过fft()计算长度为N的序列x[n]的长度为N的FFT y[k],并使用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 - 集成

当函数无法进行解析积分或很难进行解析积分时,通常会转向数值积分方法。SciPy 有许多用于执行数值积分的例程。它们中的大多数都可以在同一个scipy.integrate库中找到。下表列出了一些常用的函数。

先生号 功能说明
1

四边形

单一集成

2

双四元组

双整合

3

tpl四边形

三重整合

4

四方

n重多重积分

5

固定四边形

高斯求积,n 阶

6

正交

高斯求积至公差

7

龙贝格

龙伯格整合

8

特拉兹

梯形法则

9

库姆特拉普兹

累积计算积分的梯形规则

10

辛普斯

辛普森法则

11

抢劫

龙伯格整合

12

多整型

解析多项式积分 (NumPy)

13

聚一维

polyint 的辅助函数 (NumPy)

单积分

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

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

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

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

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、tplquadnquad中。这些函数分别集成四个或六个参数。所有内积分的极限都需要定义为函数。

二重积分

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 还有许多其他积分例程,包括执行 n 重多重积分的 nquad,以及实现各种积分算法的其他例程。然而,quad 和 dblquad 将满足我们对数值积分的大部分需求。

SciPy - 插值

在本章中,我们将讨论插值在 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 表示插值技术的类型。“线性”、“最近”、“零”、“线性”、“二次”、“三次”是插值的几种技术。

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

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 = 无,ext = 0,check_finite = False)。

参数- 以下是单变量样条的参数。

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

  • 'w' - 指定样条拟合的权重。必须是积极的。如果没有(默认),则权重全部相等。

  • 's' - 通过指定平滑条件来指定结数。

  • 'k' - 平滑样条的次数。必须 <= 5。默认值为 k = 3,三次样条。

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

    • 如果 ext = 0 或“外推”,则返回外推值。

    • 如果 ext = 1 或“零”,则返回 0

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

    • 如果 ext = 3 of '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
  • 室内设计语言
  • 矩阵市场
  • 海浪
  • 阿尔夫
  • NETCDF 等

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

MATLAB

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

先生。没有。 功能说明
1

装载垫

加载 MATLAB 文件

2

保存马特

保存 MATLAB 文件

3

何斯马特

列出 MATLAB 文件内的变量

让我们考虑下面的例子。

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 LAPACKBLAS库构建。它具有非常快的线性代数能力。所有这些线性代数例程都需要一个可以转换为二维数组的对象。这些例程的输出也是一个二维数组。

SciPy.linalg 与 NumPy.linalg

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

线性方程组

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

作为示例,假设需要求解以下联立方程。

x + 3y + 5z = 10

2x + 5y + z = 8

2x + 3y + 8z = 3

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

$$\begin{bmatrix} x\\ y\\ z \end{bmatrix} = \begin{bmatrix} 1 & 3 & 5\\ 2 & 5 & 1\\ 2 & 3 & 8 \end{bmatrix}^ {-1} \begin{bmatrix} 10\\ 8\\ 3 \end{bmatrix} = \frac{1}{25} \begin{bmatrix} -232\\ 129\\ 19 \end{bmatrix} = \开始{bmatrix} -9.28\\ 5.16\\ 0.76 \end{bmatrix}.$$

但是,最好使用linalg.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()

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

图像模糊操作

西格玛值表示模糊程度,等级为 5。我们可以通过调整sigma值来看到图像质量的变化。有关模糊的更多详细信息,请单击 → DIP(数字图像处理)教程。

边缘检测

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

什么是边缘检测?

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

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

  • 索贝尔
  • 精明的
  • 普鲁伊特
  • 罗伯茨
  • 模糊逻辑方法

让我们考虑下面的例子。

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()

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

边缘检测-2

SciPy - 优化

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})$$

当 xi = 1 时,该函数的最小值为 0。

Nelder–Mead 单纯形算法

在以下示例中,minimum() 例程与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 分别使用Powell 的混合方法和MINPACK 的Levenberg-Marquardt 方法。

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

x 2 + 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 - 统计

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

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

先生。没有。 类别和描述
1

rv_连续

用于子类化的通用连续随机变量类

2

rv_离散

用于子类化的通用离散随机变量类

3

rv_直方图

生成由直方图给出的分布

正态连续随机变量

随机变量X可以取任意值的概率分布是连续随机变量。location (loc) 关键字指定平均值。尺度(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])

上述输出不可重现。要生成相同的随机数,请使用种子函数。

均匀分布

可以使用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. ])

描述性统计

Min、Max、Mean 和 Variance 等基本统计数据将 NumPy 数组作为输入并返回各自的结果。下表描述了scipy.stats 包中可用的一些基本统计函数。

先生。没有。 功能说明
1

描述()

计算传递数组的几个描述性统计数据

2

平均()

计算沿指定轴的几何平均值

3

hmean()

计算沿指定轴的调和平均值

4

峰度()

计算峰度

5

模式()

返回模态值

6

倾斜()

测试数据的偏度

7

f_oneway()

执行单向方差分析

8

IQR()

计算沿指定轴的数据的四分位数范围

9

zscore()

计算样本中每个值相对于样本均值和标准差的 z 分数

10

扫描电镜()

计算输入数组中值的平均值(或测量标准误差)的标准误差

其中几个函数在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 imp