Python SciPy 教程：是什么、库、函数和示例

Python 中的 SciPy

Python 中的 SciPy 是一个开源库，用于解决数学、科学、工程和技术问题。它允许用户通过一系列高级 Python 命令来操作和可视化数据。SciPy 是在 Python NumPy 扩展的基础上构建的。SciPy 也发音为“Sigh Pi”。

SciPy 的子包

• 文件输入/输出 – scipy.io
• 特殊函数 – scipy.special
• 线性代数运算 – scipy.linalg
• 插值 – scipy.interpolate
• 优化和拟合 – scipy.optimize
• 统计和随机数 – scipy.stats
• 数值积分 – scipy.integrate
• 快速傅里叶变换 – scipy.fftpack
• 信号处理 – scipy.signal
• 图像处理 – scipy.ndimage

为什么使用 SciPy

• SciPy 包含各种子包，有助于解决科学计算中最常见的问题。
• Python 中的 SciPy 包是除 C/C++ 的 GNU Scientific Library 或 Matlab 之外最常用的科学库。
• 易于使用和理解，并且计算能力强。
• 它可以操作 NumPy 库的数组。

NumPy vs SciPy

NumPy

• NumPy 使用 C 编写，用于数学或数值计算。
• 它比其他 Python 库更快
• NumPy 是数据科学执行基本计算的最有用的库。
• NumPy 仅包含数组数据类型，它执行最基本的操作，如排序、整形、索引等。

SciPy

• SciPy 构建于 NumPy 之上
• Python 中的 SciPy 模块是功能齐全的线性代数版本，而 NumPy 只包含一些功能。
• 许多新的数据科学功能在 SciPy 中可用，而不是 NumPy。

SciPy – 安装和环境设置

您还可以通过 pip 在 Windows 中安装 SciPy

Python3 -m pip install --user numpy scipy 

在 Linux 上安装 Scipy

sudo apt-get install  python-scipy python-numpy

在 Mac 上安装 SciPy

sudo port install py35-scipy py35-numpy

在我们开始学习 SciPy Python 之前，您需要了解 NumPy 的基本功能和不同类型的数组。

导入 SciPy 模块和 Numpy 的标准方法

from scipy import special   #same for other modules
import numpy as np

文件输入/输出包

Scipy 的 I/O 包提供了广泛的函数来处理不同的文件格式，如 Matlab、Arff、Wave、Matrix Market、IDL、NetCDF、TXT、CSV 和二进制格式。

我们来看一个在 MatLab 中经常使用的 Python SciPy 文件格式示例。

 import numpy as np
 from scipy import io as sio
 array = np.ones((4, 4))
 sio.savemat('example.mat', {'ar': array}) 
 data = sio.loadmat(‘example.mat', struct_as_record=True)
 data['ar']

输出

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

代码解释

• 第 1 和第 2 行： 使用 I/O 包和 Numpy 导入必要的 SciPy 库。
• 第 3 行： 创建一个 4x4 的二维全 1 数组。
• 第 4 行： 将数组保存在 example.mat 文件中。
• 第 5 行： 从 example.mat 文件获取数据。
• 第 6 行： 打印输出。

特殊函数包

• scipy.special 包包含大量的数学物理函数。
• SciPy 的特殊函数包括立方根、指数、对数指数和、Lambert、排列组合、Gamma、Bessel、超几何、Kelvin、beta、抛物柱面、相对误差指数等。
• 要获得这些函数的一行描述，请在 Python 控制台中输入：

help(scipy.special)	
Output:
NAME
    scipy.special

DESCRIPTION
    ========================================
    Special functions (:mod:`scipy.special`)
    ========================================
     
    .. module:: scipy.special
     
    Nearly all of the functions below are universal functions and follow
    broadcasting and automatic array-looping rules. Exceptions are noted.

立方根函数

立方根函数计算值的立方根。

语法

scipy.special.cbrt(x)

示例

from scipy.special import cbrt
#Find cubic root of 27 & 64 using cbrt() function
cb = cbrt([27, 64])
#print value of cb
print(cb)

输出： array([3., 4.])

指数函数

指数函数计算每个元素的 10**x。

示例

from scipy.special import exp10
#define exp10 function and pass value in its
exp = exp10([1,10])
print(exp)

输出： [1.e+01 1.e+10]

排列与组合

SciPy 还提供了计算排列和组合的功能。

组合 – scipy.special.comb(N,k)

示例

from scipy.special import comb
#find combinations of 5, 2 values using comb(N, k)
com = comb(5, 2, exact = False, repetition=True)
print(com)

输出 15.0

排列 –

scipy.special.perm(N,k)

示例

from scipy.special import perm
#find permutation of 5, 2 using perm (N, k) function
per = perm(5, 2, exact = True)
print(per)

输出 20

对数指数和函数

Log Sum Exponential 计算输入元素指数和的对数。

语法

scipy.special.logsumexp(x) 

Bessel 函数

N 阶整数阶计算函数

语法

scipy.special.jn()

SciPy 线性代数

• SciPy 的线性代数是 BLAS 和 ATLAS LAPACK 库的实现。
• 线性代数的性能比 BLAS 和 LAPACK 快得多。
• 线性代数例程接受二维数组对象，输出也是二维数组。

现在让我们用 scipy.linalg 进行一些测试，

计算二维矩阵的行列式，

from scipy import linalg
import numpy as np
#define square matrix
two_d_array = np.array([ [4,5], [3,2] ])
#pass values to det() function
linalg.det( two_d_array )

输出 -7.0

逆矩阵 –

scipy.linalg.inv()

SciPy 的逆矩阵计算任何方阵的逆。

让我们看看，

from scipy import linalg
import numpy as np
# define square matrix
two_d_array = np.array([ [4,5], [3,2] ])
#pass value to function inv()
linalg.inv( two_d_array )

输出

array( [[-0.28571429,  0.71428571],
       [ 0.42857143, -0.57142857]] )

特征值和特征向量

scipy.linalg.eig()

• 线性代数中最常见的问题是特征值和特征向量，可以使用 eig() 函数轻松解决。
• 现在我们来找出二维方阵 (X) 的特征值及其对应的特征向量。

示例

from scipy import linalg
import numpy as np
#define two dimensional array
arr = np.array([[5,4],[6,3]])
#pass value into function
eg_val, eg_vect = linalg.eig(arr)
#get eigenvalues
print(eg_val)
#get eigenvectors
print(eg_vect)

输出

[ 9.+0.j -1.+0.j] #eigenvalues
 [ [ 0.70710678 -0.5547002 ] #eigenvectors
   [ 0.70710678  0.83205029] ]

离散傅里叶变换 – scipy.fftpack

• DFT 是一种数学技术，用于将空间数据转换为频率数据。
• FFT（快速傅里叶变换）是计算 DFT 的算法。
• FFT 应用于多维数组。
• 频率定义了一个特定时间段内的信号或波长数量。

示例： 以 Matplotlib 库为例，展示一个波形。我们以 sin(20 × 2πt) 为例，这是一个简单的周期函数。

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np 

#Frequency in terms of Hertz
fre  = 5 
#Sample rate
fre_samp = 50
t = np.linspace(0, 2, 2 * fre_samp, endpoint = False )
a = np.sin(fre  * 2 * np.pi * t)
figure, axis = plt.subplots()
axis.plot(t, a)
axis.set_xlabel ('Time (s)')
axis.set_ylabel ('Signal amplitude')
plt.show()

输出

您可以看到，频率为 5 Hz，其信号在 1/5 秒内重复——这称为特定时间周期。

现在让我们利用 DFT 应用来使用这个正弦波。

from scipy import fftpack

A = fftpack.fft(a)
frequency = fftpack.fftfreq(len(a)) * fre_samp
figure, axis = plt.subplots()

axis.stem(frequency, np.abs(A))
axis.set_xlabel('Frequency in Hz')
axis.set_ylabel('Frequency Spectrum Magnitude')
axis.set_xlim(-fre_samp / 2, fre_samp/ 2)
axis.set_ylim(-5, 110)
plt.show()

输出

• 您可以清楚地看到输出是一个一维数组。
• 输入中的复数值在两个点之外均为零。
• 在 DFT 示例中，我们可视化了信号的幅度。

SciPy 中的优化与拟合 – scipy.optimize

• 优化提供了用于曲线拟合、多维或标量以及根拟合最小化的有用算法。
• 让我们以标量函数为例，来寻找最小标量函数。

%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
import numpy as np

def function(a):
       return   a*2 + 20 * np.sin(a)
plt.plot(a, function(a))
plt.show()
#use BFGS algorithm for optimization
optimize.fmin_bfgs(function, 0) 

输出

优化成功终止。

当前函数值：-23.241676

迭代次数：4

函数评估次数：18

梯度评估次数：6

array([-1.67096375])

• 在此示例中，优化是通过从初始点使用梯度下降算法完成的。
• 但是可能存在的问题是局部最小值而不是全局最小值。如果我们找不到全局最小值的邻域，那么我们需要应用全局优化，并使用 basinhopping() 作为全局最小化函数，它结合了局部优化器。

optimize.basinhopping(function, 0)

输出

fun: -23.241676238045315
 lowest_optimization_result:
      fun: -23.241676238045315
 hess_inv: array([[0.05023331]])
      jac: array([4.76837158e-07])
  message: 'Optimization terminated successfully.'
     nfev: 15
      nit: 3
     njev: 5
   status: 0
  success: True
        x: array([-1.67096375])
                    message: ['requested number of basinhopping iterations completed successfully']
      minimization_failures: 0
                       nfev: 1530
                        nit: 100
                       njev: 510
               x: array([-1.67096375])

Nelder-Mead 算法

• Nelder-Mead 算法通过 method 参数进行选择。
• 对于行为良好的函数，它提供了最直接的最小化方法。
• Nelder-Mead 算法不用于梯度评估，因为它可能需要更长的时间来找到解决方案。

import numpy as np
from scipy.optimize import minimize
#define function f(x)
def f(x):   
    return .4*(1 - x[0])**2
  
optimize.minimize(f, [2, -1], method="Nelder-Mead")

输出

final_simplex: (array([[ 1.        , -1.27109375],
       [ 1.        , -1.27118835],
       [ 1.        , -1.27113762]]), array([0., 0., 0.]))
           fun: 0.0
       message: 'Optimization terminated successfully.'
          nfev: 147
           nit: 69
        status: 0
       success: True
             x: array([ 1.        , -1.27109375])

SciPy 图像处理 – scipy.ndimage

• scipy.ndimage 是 SciPy 的一个子模块，主要用于执行图像相关操作。
• ndimage 意为“n”维图像。
• SciPy 图像处理提供几何变换（旋转、裁剪、翻转）、图像滤波（锐化和去噪）、显示图像、图像分割、分类和特征提取。
• SciPy 中的MISC 包包含预建图像，可用于执行图像操作任务。

示例： 让我们以图像的几何变换为例。

from scipy import misc
from matplotlib import pyplot as plt
import numpy as np
#get face image of panda from misc package
panda = misc.face()
#plot or show image of face
plt.imshow( panda )
plt.show()

输出

现在我们垂直翻转当前图像。

#Flip Down using scipy misc.face image  
flip_down = np.flipud(misc.face())
plt.imshow(flip_down)
plt.show()

输出

示例： 使用 Scipy 旋转图像，

from scipy import ndimage, misc
from matplotlib import pyplot as plt
panda = misc.face()
#rotatation function of scipy for image – image rotated 135 degree
panda_rotate = ndimage.rotate(panda, 135)
plt.imshow(panda_rotate)
plt.show()

输出

SciPy 集成 – 数值积分

• 当我们在无法进行解析积分的函数上进行积分时，我们需要转向数值积分。
• SciPy 提供了使用数值积分功能来积分函数。
• scipy.integrate 库提供单重积分、重积分、三重积分、多重积分、高斯积分、Romberg、梯形和辛普森规则。

示例： 现在以单重积分为例。

这里 a 是上限，b 是下限。

from scipy import integrate
# take f(x) function as f
f = lambda x : x**2
#single integration with a = 0 & b = 1  
integration = integrate.quad(f, 0 , 1)
print(integration)

输出

(0.33333333333333337, 3.700743415417189e-15)

这里函数返回两个值，第一个值是积分值，第二个值是积分的估计误差。

示例：现在以 SciPy 的双重积分为例。我们计算以下方程的双重积分，

from scipy import integrate
import numpy as np
#import square root function from math lib
from math import sqrt
# set  fuction f(x)
f = lambda x, y : 64 *x*y
# lower limit of second integral
p = lambda x : 0
# upper limit of first integral
q = lambda y : sqrt(1 - 2*y**2)
# perform double integration
integration = integrate.dblquad(f , 0 , 2/4,  p, q)
print(integration)

输出

(3.0, 9.657432734515774e-14)

您可以看到上面的输出与之前的输出相同。

摘要

• SciPy（发音为“Sigh Pi”）是一个开源的 Python 库，用于数学、科学计算、工程和技术计算。
• SciPy 包含各种子包，有助于解决科学计算中最常见的问题。
• SciPy 构建于 NumPy 之上

包名	描述
scipy.io	文件输入/输出
scipy.special	特殊函数
scipy.linalg	线性代数运算
scipy.interpolate	插值
scipy.optimize	优化和拟合
scipy.stats	统计和随机数
scipy.integrate	数值积分
scipy.fftpack	快速傅里叶变换
scipy.signal	信号处理
scipy.ndimage	图像处理 –