

高级的科学计算库——Scipy

数据皮皮侠

2019-10-17

导读：Scipy是一个高级的科学计算库，它和Numpy联系很密切，Scipy一般都是操控Numpy数组来进行科学计算，Scipy让Python成为了半个MATLAB。

本期目录

Oct.18, 2019

一、简介

二、安装

三、常用子模块

四、应用

4.1简介

4.2统计假设与检验 stats包

4.3信号特征

4.4寻优

4.5求解

4.6曲线拟合 curve-fit

4.7插值

4.8模式聚类

简介

Scipy是一个高级的科学计算库，它和Numpy联系很密切，Scipy一般都是操控Numpy数组来进行科学计算，Scipy让Python成为了半个MATLAB。Scipy包含的功能有最优化、线性代数、积分、插值、拟合、特殊函数、快速傅里叶变换、信号处理和图像处理、常微分方程求解和其他科学与工程中常用的计算，而这些功能都是我们在之后进行数据分析需要的。

安装

Scipy依赖于Numpy库，因此安装Scipy时应先安装Numpy库，Scipy安装与其他库一样，可通过pip install Scipy安装，也可以自行下载源代码，然后用pip install 路径+文件名全称（包括.后缀文件名）进行安装，源码下载链接：https://pypi.python.org/pypi/scipy/1.0.0，选择对应版本下载即可。

常用子模块

接下来我们简单的介绍一下这些模块的应用！

应用

4.1 scipy.io文件的输入和输出

scipy.io模块提供了多种功能来解决不同格式的文件输入和输出，包括Matlab，Wave，Arff，Matrix Market等等，最常见的是Matlab格式的。

（1）导入相关依赖包

from scipy import io as spioimport numpy as np

（2）创建一个数组并保存到文件中

a = np.ones((3,3))spio.savemat('file.mat',{'a':a})data = spio.loadmat('file.mat', struct_as_record=True)data['a']

载入txt文件：

numpy.loadtxt()/numpy.savetxt()

智能导入文本/csv文件：

numpy.genfromtxt()/numpy.recfromcsv()

高速，有效率但numpy特有的二进制格式：

numpy.save()/numpy.load()

4.2 统计假设与检验 stats包

stats提供了产生连续性分布的函数：均匀分布(uniform)、正态分布（norm）、贝塔分布（beta）；

产生离散分布的函数：伯努利分布（bernoulli）、几何分布（geom）泊松分布 poisson。

使用时，调用分布的rvs方法即可。

例如：

import scipy.stats as statsx=stats.uniform.rvs(size=20)#产生20个在[0，1]均匀分布的随机数

y=stats.beta.rvs(size=20,a=3,b=4)#产生20个服从参数a=3，b=4的贝塔分布随机数

z=stats.norm.rvs(size=20,loc=0,scale=1)#产生了20个服从[0，1]正态分布的随机数

x=stats.poisson.rvs(0.6,loc=0,size=20)#产生poisson分布

如果假设给定的样本服从某种分布，那么如何验证呢？

import numpy as npimport scipy.stats as statsnormDist=stats.norm(loc=2.5,scale=0.5)z=normDist.rvs(size=400)mean=np.mean(z)med=np.median(z)dev=np.std(z)print('mean=',mean,' med=',med,' dev=',dev)

设z是实验获得的数据，如何验证它是否是正态分布的？

statVal, pVal = stats.kstest(z,'norm',(mean,dev))print('pVal=',pVal)

计算得到

mean= 2.4788577873926516

med= 2.473711775463392

dev= 0.5343868193748538

pVal= 0.8908805240503633

结论：我们接受假设，即z数据是服从正态分布的

4.3 信号特征

低频的原始信号，加高频的噪声，如何去掉噪声？

快速傅里叶变换 FFT应用范围非常广，在物理学、化学、电子通讯、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等。原理是将时域中的测量信号转换到频域，然后再将得到的频域信号合成为时域中的信号；时域信号转换为频域信号时，将信号分解成幅值谱，显示与频率对应的幅值大小。一个信号是由多种频率的信号合成的，将这些信号分离，就能去掉信号中的噪声。

Scipy类库中提供了快速傅里叶变换包fftpack。

接下来我们通过一个FFT应用案例—产生带噪声的信号来进一步学习这个包吧~

import numpy as npimport matplotlib.pyplot as pltfrom scipy import fftpack as ffttimeStep = 0.02  # 定义采样点间隔timeVec = np.arange(0, 20, timeStep)  #  生成采样点# 模拟带高频噪声的信号sigsig = np.sin( np.pi / 5 * timeVec) + 0.1 * np.random.randn(timeVec.size)  # 表达式中后面一项为噪声plt.plot(timeVec, sig)plt.show()

输出结果

FFT信号变换 sig已知

n=sig.sizesig_fft = fft.fft(sig) # 正变换后的结果保存在 sig_fft中sampleFreq = fft.fftfreq(n, d=timeStep) # 获得每个采样点频率幅值pidxs = np.where(sampleFreq > 0) # 结果是对称的，仅需要使用谱的正值部分来找出信号频率freqs = sampleFreq[pidxs] # 取得所有正值部分power = np.abs(sig_fft)[pidxs]  # 找到对应的频率振幅值freq = freqs[power.argmax()]   # 假设信噪比足够强，则最大幅值对应的频率，就是信号的频率sig_fft[np.abs(sampleFreq) > freq] = 0 # 舍弃所有非信号频率main_sig = fft.ifft(sig_fft) # 用傅立叶反变换，重构去除噪声的信号plt.plot(timeVec, main_sig, linewidth=3)

输出结果

4.4 寻优

f(x)=x2-4*x+8 在x=2的位置，函数有最小值4。

这个题可以通过scipy.optimize包提供的求极值的函数fmin来计算：

from scipy.optimize import fminimport numpy as npdef f(x):    return x**2-4*x+8print (fmin(f, 0))

结果如下：

Optimization terminated successfully.

Current function value: 4.000000

Iterations: 27

Function evaluations: 54

这个函数不仅可以进行单维寻优，也可以进行多维寻优哦~例如，求解z=x2+y2+8的最小值

from scipy.optimize import fmindef myfunc(p): # 注意定义    x,y=p    return x**2+y**2+8p=(1,1)print (fmin(myfunc, p ))

结果如下：

Optimization terminated successfully.

Current function value: 8.000000

Iterations: 38

Function evaluations: 69

[-2.10235293e-05 2.54845649e-05]

4.5 求解

求解一元高次方程的根

例如：

from scipy import optimizeimport numpy as npdef f(x):    return x**2 + 20 * np.sin(x)ans = optimize.fsolve(f, -4)print(ans)print(f(ans))

结果为

[-2.75294663]

[1.687539e-14]

求解非线性方程组

scipy. optimize的fsolve函数也可以方便用于求解非线性方程组。

求解原则：将方程组的右边全部规划为0，将方程组定义为如下格式的python函数：

def f(x):      x1,x2,…, xn=x  return [f1(x1, x2,…, xn), f2(x1,x2,…, xn),….]

例如：

from scipy.optimize import fsolvefrom math import *def f(x):    x0, x1,x2 = x     return [2*x1+3, 4*x0*x0 + sin(x1*x2), x1*x2/2 - 3]ans = fsolve(f, [1.0,1.0,1.0])print (ans)print (f(ans))

结果为

[-0.26429884 -1.5 -4. ]

[0.0, 1.1482453876610066e-10, 6.4002136923591024e-12]

4.6 曲线拟合 curve-fit

给定的y=ax+b函数上的一系列采样点，并在这些采样点上增加一些噪声，然后利用scipy optimize包中提供的curve_fit方法，求解系数a和b

from scipy import optimizeimport matplotlib.pyplot as pltimport numpy as npdef f(x,a,b):    return a*x + bx = np.linspace(-10, 10, 30)y = f(x,2,1)ynew= y+ 3*np.random.normal(size=x.size)  # 产生带噪声的数据点popt, pcov = optimize.curve_fit(f,x,ynew)print(popt)print (pcov)plt.plot(x,y,color='r',label='orig')plt.plot(x,popt[0]*x+popt[1],color='b',label='fitting')plt.legend(loc='upper left')plt.scatter(x,ynew)plt.show()

结果为

popt=[2.06246704 1.65798483]

pcov=[[6.31189436e-03 1.45358283e-10]

[1.45358283e-10 2.24906575e-01]]

输出结果

popt列表包含每个参数的拟合值，此例求得的a为2.06246704，b为1.65798483。pcov列表的对角元素是每个参数的方差。通过方差，可以评判拟合的质量，方差越小，拟合越可靠。

4.7 插值

根据现有的试验点值，去预测中间的点值

采用线性、两次样条、三次样条插值

我们同样通过一个案例来学习

首先在[0,10π]区间里等间距产生了20个采样点，并计算其sin值，模拟试验采集得到的20个点，采用线性、两次样条、三次样条插值函数插值拟合原函数sin(x)

import numpy as npfrom scipy.interpolate import interp1dimport matplotlib.pyplot as plt
x=np.linspace(0,10*np.pi,20)  #产生20个点y=np.sin(x)
#  x,y 现在是假设的采样点
fl=interp1d(x,y,kind='linear')  # 线性插值函数fq=interp1d(x,y,kind='quadratic') # 二次样条插值fc=interp1d(x,y,kind='cubic') # 三次样条插值xint=np.linspace(x.min(),x.max(),1000) # 产生插值点yintl=fl(xint) # 由线性插值得到的函数值yintq=fq(xint) # 由二次样条插值函数计算得到的函数值yintc=fc(xint) # 由三次样条插值函数计算得到的函数值plt.plot(xint,yintl,color='r', linestyle='--',label='linear')plt.plot(xint,yintq,color='b' ,label='quadr')plt.plot(xint,yintc,color='b' ,linestyle='-.',label='cubic')plt.legend()plt.show()

输出结果

4.8 模型聚类

Scipy聚类分析中主要提供了矢量量化分析和系统聚类法

import numpy as npfrom scipy.cluster import vqimport matplotlib.pyplot as pltclass1=np.random.randn(30,2)+10  # 产生第一个正态分布类，基础抬高10class2=np.random.randn(40,2)-10  # 产生第二个正态分布类，基础降低10class3=np.random.randn(50,2)     # 产生第三个正态分布类data=np.vstack([class1,class2,class3]) #将数据叠合到一起，形成一个矩阵centroid, var=vq.kmeans(data,3)   # 用k均值聚类法聚类，指定按3个类别聚类，获取类中心和方差key,distance=vq.vq(data,centroid)  # 根据聚类中心，将不同的样本分类vqclass1=data[key==0]  # 取出类别0vqclass2=data[key==1]vqclass3=data[key==2]
plt.scatter(vqclass1[:,0],vqclass1[:,1],marker='o', color="r",label='c1')  # 为类0 制图plt.scatter(vqclass2[:,0],vqclass2[:,1],marker='1', color="g" ,label='c2')plt.scatter(vqclass3[:,0],vqclass3[:,1],marker='2', color="b",label='c3')plt.legend(loc='upper left')plt.show()