
本期目录 Oct.18, 2019 |
一、简介 |
二、安装 |
三、常用子模块 |
四、应用 4.1简介 4.2统计假设与检验 stats包 4.3信号特征 4.4寻优 4.5求解 4.6曲线拟合 curve-fit 4.7插值 4.8模式聚类 |
01
简介
Scipy是一个高级的科学计算库,它和Numpy联系很密切,Scipy一般都是操控Numpy数组来进行科学计算,Scipy让Python成为了半个MATLAB。Scipy包含的功能有最优化、线性代数、积分、插值、拟合、特殊函数、快速傅里叶变换、信号处理和图像处理、常微分方程求解和其他科学与工程中常用的计算,而这些功能都是我们在之后进行数据分析需要的。
02
安装
Scipy依赖于Numpy库,因此安装Scipy时应先安装Numpy库,Scipy安装与其他库一样,可通过pip install Scipy安装,也可以自行下载源代码,然后用pip install 路径+文件名全称(包括.后缀文件名)进行安装,源码下载链接:https://pypi.python.org/pypi/scipy/1.0.0,选择对应版本下载即可。
03
常用子模块

接下来我们简单的介绍一下这些模块的应用!
04
应用
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 pltx=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()

输出结果
END

作者有话说
Scipy的简单介绍及应用到这里就结束啦,能看到最后的大家都棒棒哒!总的来说,Scipy库的强大不是这一篇推文可以概括完的,但还是希望这篇推文能给正在学习python的你带来帮助!本期内容结束,感想阅读和关注~
本期作者:张怡
















