数据分析与可视化(工具篇)——Scipy 使用指北

scipy包含各种专用于科学计算中常见问题的工具箱。其不同的子模块对应不同的应用,如插值、积分、优化、图像处理、统计、特殊函数等。

scipy可以与其他标准科学计算库进行比较,例如 GSL(用于 C 和 C++ 的 GNU 科学库)或 Matlab 的工具箱。scipy是 Python 中科学计算的核心包;它旨在有效地在numpy 数组上运行,以便 numpy 和 scipy 共同使用。

 

目录

1. 文件输入/输出:scipy.io

2. 特殊函数:scipy.special

3. 线性代数运算:scipy.linalg

4. 插值:scipy.interpolate

5. 优化和拟合:scipy.optimize

6. 统计和随机数:scipy.stats

6.1 分布:直方图和概率密度函数

6.2 平均值、中位数和百分位数

6.3 统计检验

7. 数值积分:scipy.integrate

7.1 函数积分

7.2 积分微分方程

8. 快速傅里叶变换:scipy.fftpack

9. 信号处理:scipy.signal

10. 图像处理:scipy.ndimage

10.1 图像的几何变换

10.2 图像过滤

10.3 数学形态学


scipy 非常丰富,这里我们只介绍一些重点部分,帮助我们了解如何将scipy用于科学计算。

scipy 由特定于任务的子模块组成:

scipy.cluster矢量量化/Kmeans
scipy.constants物理和数学常数
scipy.fftpack傅里叶变换
scipy.integrate积分
scipy.interpolate插值
scipy.io数据输入输出
scipy.linalg线性代数
scipy.ndimagen维图像包
scipy.odrOrthogonal distance regression
scipy.optimize优化
scipy.signal信号处理
scipy.sparse稀疏矩阵
scipy.spatial空间数据结构和算法
scipy.special任何特殊的数学函数
scipy.stats统计数据

它们都依赖于numpy,但大多是相互独立的。导入 Numpy 和这些 Scipy 模块的标准方法是: 

>>> import numpy as np
>>> from scipy import stats

主scipy命名空间大部分包含了真正的numpy函数(scipy.cos , np.cos)。这些都是由于历史原因造成的;没有理由在代码中使用import scipy。

1. 文件输入/输出:scipy.io

Matlab 文件:加载和保存:

>>> from scipy import io as spio
>>> a = np.ones((3, 3))
>>> spio.savemat('file.mat', {'a': a}) 
>>> data = spio.loadmat('file.mat')
>>> data['a']
array([[1.,  1.,  1.],
       [1.,  1.,  1.],
       [1.,  1.,  1.]])

 Python/Matlab 不匹配,例如

>>> a = np.ones(3)
>>> a
array([1.,  1.,  1.])
>>> spio.savemat('file.mat', {'a': a})
>>> spio.loadmat('file.mat')['a']
array([[1.,  1.,  1.]])

 图像文件:读取图像:

>>> import imageio
>>> imageio.imread('fname.png')    
Array(...)
>>> # Matplotlib also has a similar function
>>> import matplotlib.pyplot as plt
>>> plt.imread('fname.png')    
array(...)
  • 加载文本文件:numpy.loadtxt()/numpy.savetxt()
  • 巧妙加载文本/csv文件: numpy.genfromtxt()/numpy.recfromcsv()
  • 快速高效,但特定于numpy的二进制格式: numpy.save()/numpy.load()
  • scikit-image 中更高级的图像输入/输出: skimage.io

2. 特殊函数:scipy.special

特殊函数是超越函数。scipy.special模块的文档写得很好,所以我们不会在这里列出所有的函数。常用的有:

贝塞尔函数,如scipy.special.jn()(nth integer order Bessel function)
椭圆函数(scipy.special.ellipj()对于雅可比椭圆函数,...)
Gamma function: scipy.special.gamma(),还要注意 scipy.special.gammaln()这将使 Gamma 的对数具有更高的数值精度。
Erf,高斯曲线下的面积: scipy.special.erf()

3. 线性代数运算:scipy.linalg

scipy.linalg模块提供标准的线性代数运算。

from scipy import linalg
>>> arr = np.array([[1, 2],
...                 [3, 4]])
>>> linalg.det(arr)
-2.0
>>> arr = np.array([[3, 2],
...                 [6, 4]])
>>> linalg.det(arr) 
0.0
>>> linalg.det(np.ones((3, 4)))
Traceback (most recent call last):
...
ValueError: expected square matrix
  • scipy.linalg.inv()函数计算方阵的逆:
>>> arr = np.array([[1, 2],
...                 [3, 4]])
>>> iarr = linalg.inv(arr)
>>> iarr
array([[-2. ,  1. ],
       [ 1.5, -0.5]])
>>> np.allclose(np.dot(arr, iarr), np.eye(2))
True
  • 可以使用更高级的操作,例如奇异值分解(SVD):
arr = np.array([[3, 2], ... [6, 4]])
>>> linalg.inv(arr) Traceback (most recent call last): ... ...
LinAlgError: singular matrix
  • 得到的数组谱为:
>>> spec array([14.88982544, 0.45294236, 0.29654967])
  • 原始矩阵可以通过svdwith的输出的矩阵乘法重新组合

np.dot:

>>> sarr = np.diag(spec)
>>> svd_mat = uarr.dot(sarr).dot(vharr) >>> np.allclose(svd_mat, arr) True
  • SVD 常用于统计和信号处理。许多其他标准分解(QR、LU、Cholesky、Schur)以及线性系统的求解器都包含在scipy.linalg.

4. 插值:scipy.interpolate

scipy.interpolate对于从实验数据拟合函数非常有用,因此可以对不存在测度的点进行评估。该模块基于FITPACK Fortran。

通过想象接近正弦函数的实验数据:

>>> measured_time = np.linspace(0, 1, 10)
>>> noise = (np.random.random(10)*2 - 1) * 1e-1
>>> measures = np.sin(2 * np.pi * measured_time) + noise

scipy.interpolate.interp1d 可以创建一个线性插值函数:

>>> from scipy.interpolate import interp1d
>>> linear_interp = interp1d(measured_time, measures)

可以在感兴趣的点估算结果:

>>> interpolation_time = np.linspace(0, 1, 50)
>>> linear_results = linear_interp(interpolation_time)

三次插值也可以通过提供kind可选关键字参数来选择:

>>> cubic_interp = interp1d(measured_time, measures, kind='cubic')
>>> cubic_results = cubic_interp(interpolation_time)

scipy.interpolate.interp2d类似于 scipy.interpolate.interp1d,但适用于二维数组。请注意,对于interp系列,插值点必须保持在给定数据点的范围内。

5. 优化和拟合:scipy.optimize

优化是找到最小化或相等的数值解的问题。

scipy.optimize模块提供函数最小化(标量或多维)、曲线拟合和求根的算法。

>>> from scipy import optimize

6. 统计和随机数:scipy.stats

该模块scipy.stats包含随机过程的统计工具和概率描述。各种随机过程的随机数生成器可以在 numpy.random 中找到。

6.1 分布:直方图和概率密度函数

对给定随机过程的观察,它们的直方图是随机过程 PDF(概率密度函数)的估计量:

>>> samples = np.random.normal(size=1000)
>>> bins = np.arange(-4, 5)
>>> bins
array([-4, -3, -2, -1,  0,  1,  2,  3,  4])
>>> histogram = np.histogram(samples, bins=bins, density=True)[0]
>>> bins = 0.5*(bins[1:] + bins[:-1])
>>> bins
array([-3.5, -2.5, -1.5, -0.5,  0.5,  1.5,  2.5,  3.5])
>>> from scipy import stats
>>> pdf = stats.norm.pdf(bins)  # norm is a distribution object

>>> plt.plot(bins, histogram) 
[<matplotlib.lines.Line2D object at ...>]
>>> plt.plot(bins, pdf) 
[<matplotlib.lines.Line2D object at ...>]

分布对象

scipy.stats.norm是一个分布对象:每个分布都表示为一个对象。norm 是正态分布,包含有 PDF、CDF 等等。

如果我们知道随机过程属于给定的随机过程族,例如正态过程,我们可以对观测值进行最大似然拟合以估计基础分布的参数。在这里,我们将正态过程拟合到观察到的数据:

loc, std = stats.norm.fit(samples)
print(loc,std)
-0.004822316383745015 1.0050446891199836

6.2 平均值、中位数和百分位数

均值

>>> np.mean(samples)     
-0.004822316383745015

中位数

>>> np.median(samples)     
 0.03361659877914626

与平均值不同,中位数对分布的尾部不敏感。

中位数也是百分位数 50,因为 50% 的观察值低于它:

>>> stats.scoreatpercentile(samples, 50)     
 0.03361659877914626

同样,我们可以计算百分位数 90:

>>> stats.scoreatpercentile(samples, 90)     
1.2565894469998924

百分位数是累积分布函数估计量:。

6.3 统计检验

统计检验是一个决策指标。例如,如果我们有两组观察值,我们假设它们是从高斯过程生成的,我们可以使用 T 检验来确定两组观察值的均值是否显着不同:

a = np.random.normal(0, 1, size=100)
b = np.random.normal(1, 1, size=10)
stats.ttest_ind(a, b)
Out[14]: Ttest_indResult(statistic=-3.6062755503726986, pvalue=0.0004718656448315962)

输出结果由以下部分组成:

  • T统计值: 它是一个数字,其符号与两个随机过程的差值成正比,其大小与这个差值的显著性有关。
  • p值:两个过程相同的概率。如果它接近于1,那么这两个过程几乎肯定是相同的。它越接近于零,过程就越有可能有不同的方法。

7. 数值积分:scipy.integrate

7.1 函数积分

最通用的是scipy.integrate.quad(). 计算

>>> from scipy.integrate import quad
>>> res, err = quad(np.sin, 0, np.pi/2)
>>> np.allclose(res, 1)   # Res是结果,应该接近1
True
>>> np.allclose(err, 1 - res)  # Err是误差
True

7.2 积分微分方程

scipy.integrate还具有用于常微分方程 (ODE) 的例程。特别是scipy.integrate.odeint()求解以下形式的 ODE:

dy/dt = rhs(y1, y2, .., t0,...)

作为介绍,让我们计算在t=0-4之间的常微分方程 dy/dt = -2y,其初始条件y(t=0) = 1。首先需要定义计算位置导数的函数:

>>> def calc_derivative(ypos, time):
...     return -2 * ypos

8. 快速傅里叶变换:scipy.fftpack

scipy.fftpack模块计算快速傅里叶变换 (FFT) 并提供处理它们的实用程序。主要功能有:

  • scipy.fftpack.fft() 计算 FFT
  • scipy.fftpack.fftfreq() 生成采样频率
  • scipy.fftpack.ifft() 计算从频率空间到信号空间的逆 FFT

例如,一个(噪声)输入信号 ( sig) 及其 FFT:

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

9. 信号处理:scipy.signal

scipy.signal 用于典型的信号处理:一维、规则采样信号。

10. 图像处理:scipy.ndimage

scipy.ndimage 提供将 n 维数组作为图像进行操作。

10.1 图像的几何变换

改变方向,分辨率,..

>>> from scipy import misc  # Load an image
>>> face = misc.face(gray=True)

>>> from scipy import ndimage # Shift, roate and zoom it
>>> shifted_face = ndimage.shift(face, (50, 50))
>>> shifted_face2 = ndimage.shift(face, (50, 50), mode='nearest')
>>> rotated_face = ndimage.rotate(face, 30)
>>> cropped_face = face[50:-50, 50:-50]
>>> zoomed_face = ndimage.zoom(face, 2)
>>> zoomed_face.shape
(1536, 2048)
>>> plt.subplot(151)    
<matplotlib.axes._subplots.AxesSubplot object at 0x...>

>>> plt.imshow(shifted_face, cmap=plt.cm.gray)    
<matplotlib.image.AxesImage object at 0x...>

>>> plt.axis('off')
(-0.5, 1023.5, 767.5, -0.5)

>>> # etc.

10.2 图像过滤

生成一张嘈杂的脸:

>>> from scipy import misc
>>> face = misc.face(gray=True)
>>> face = face[:512, -512:]  # crop out square on right
>>> import numpy as np
>>> noisy_face = np.copy(face).astype(np.float)
>>> noisy_face += face.std() * 0.5 * np.random.standard_normal(face.shape)

在其上应用各种过滤器:

>>> blurred_face = ndimage.gaussian_filter(noisy_face, sigma=3)
>>> median_face = ndimage.median_filter(noisy_face, size=5)
>>> from scipy import signal
>>> wiener_face = signal.wiener(noisy_face, (5, 5))

在其它过滤器scipy.ndimage.filtersscipy.signal 可应用于图像。

10.3 数学形态学

数学形态学源于集合论。它表征和转换几何结构。尤其是二进制(黑白)图像,可以使用该理论进行转换:要转换的集合是相邻非零值像素的集合。该理论还扩展到灰度图像。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值