Python机器学习及分析工具：Scipy篇

2018-08-28 本文已影响70人殉道者之花火

Scipy是一个用于数学、科学、工程领域的常用软件包，可以处理插值、积分、优化、图像处理、常微分方程数值解的求解、信号处理等问题。它用于有效计算Numpy矩阵，使Numpy和Scipy协同工作，高效解决问题。
Scipy是由针对特定任务的子模块组成：

模块名	应用领域
scipy.cluster	向量计算/Kmeans
scipy.constants	物理和数学常量
scipy.fftpack	傅立叶变换
scipy.integrate	积分程序
scipy.interpolate	插值
scipy.io	数据输入输出
scipy.linalg	线性代数程序
scipy.ndimage	n维图像包
scipy.odr	正交距离回归
scipy.optimize	优化
scipy.signal	信号处理
scipy.sparse	稀疏矩阵
scipy.spatial	空间数据结构和算法
scipy.special	一些特殊的数学函数
scipy.stats	统计

scipy.io

载入和保存matlab文件

from scipy import io as spio
from numpy as np
x = np.ones((3,3))
spio.savemat('f.mat',{'a':a})
data = spio.loadmat('f.mat',struct_as_record=True)
data['a']

读取图片

from scipy import misc
misc.imread('picture')

scipy.special

special库中的特殊函数都是超越函数，所谓超越函数是指变量之间的关系不能用有限次加、减、乘、除、乘方、开方运算表示的函数。如初等函数中的三角函数、反三角函数与对数函数、指数函数都是初等超越函数，一般来说非初等函数都是超越函数。

初等函数：指由基本初等函数经过有限次四则运算与复合运算所得到的函数

下面我们对scipy.special中的部分常用函数进行说明:

$\gamma$ 函数
$\gamma$ 函数，也称伽玛函数，是由欧拉积分定义的函数，是阶乘函数在实数域与复数域上的扩展，当 $Re(z)>0$ 时, $\gamma$ 函数可以定义为： $\gamma(z)=\int_{0}^{+\infty}\frac{t^{z-1}}{e^t}\text{d}t$
由解析延拓原理可以将这个定义拓展到整个复数域上，非正整数除外。
在scipy.special中使用scipy.special.gamma()实现 $\gamma$ 函数的计算，如果对于精度有更高要求，可以使用采用对数坐标的scipy.special.gammaln()函数进行计算。
贝塞尔函数
在介绍贝塞尔函数之前，先对贝塞尔方程进行说明，一般来说，我们将形如
$x^2\frac{d^2y}{dx^2}+x\frac{dy}{dx}+(x^2-\alpha^2)y=0$
的二阶线性常微分方程称为贝塞尔方程,而贝塞尔方程的标准解函数 $y(x)$ 就是贝塞尔函数。其中，参数 $\alpha$ 被称为对应贝塞尔函数的阶。贝塞尔方程是一个二阶线性齐次常微分方程，必然存在两个线性无关的解，然而通常情况下它的解无法用初等函数表示，但是当 $\alpha=n$ 时，注意到 $x=0$ 是贝尔塞方程的正则奇点，则由常微分方程的广义幂级数解法可以得出贝塞尔方程的两个广义幂级数解：

$y_{1}=J_{n}(x)=\sum_{k=0}^{\infty}\frac{(-1)^k}{\gamma(n+k+ 1)\gamma(k+1)}(\frac{x}{2})^{2k+n}$
$y_{2}=J_{-n}(x)=\sum_{k=0}^{\infty}\frac{(-1)^k}{\gamma(-n+k+ 1)\gamma(k+1)}(\frac{x}{2})^{2k-n}$

其中 $y_{1}$ 被称为第一类贝塞尔函数， $y_{2}$ 被称为第二类贝塞尔函数（诺依曼函数）。（求解方法参考常微分方程教程 7.4 广义幂级数解法）。
贝塞尔方程是在柱坐标或球坐标下使用分离变量法求解拉普拉斯方程和亥姆霍兹方程时得到的,因此贝塞尔函数在波动问题以及各种涉及有势场的问题中占有非常重要的地位，其应用有：

在圆柱形波导中的电磁波传播问题
圆柱体中的热传导问题
圆形或环形薄膜的震动膜态分析问题
信号处理中的调频合成（FMsynthesis）
波动声学

在scipy.special中使用scipy.special.jn()计算 $n$ 阶贝塞尔函数

椭圆函数
椭圆函数是定义在有限复平面上亚纯的双周期函数。所谓的双周期是指具有两个基本周期的单复变函数，即存在 $\omega_{1},\omega_{2}$ 两个非0复数，对任意整数 $m,n$ 有：
$f(z+n\omega_{1}+m\omega_{2})=f(z)$ 于是 $\{n\omega_{1}+m\omega_{2}: n,m \in Z\}$ 构成 $f(z)$ 的全部周期。
在scipy.special中使用scipy.special.ellipj()函数计算椭圆函数。
Erf(高斯曲线的面积)
高斯曲线是指高斯分布也就是我们常说的正态分布
$X～N(\mu,\sigma^2)$
其概率密度函数为：
$f(x)=\frac{1} {\sigma\sqrt{2\pi}}e_{-\frac{(x-\mu)^2}{2\sigma^2}}$
其概率密度函数曲线就是高斯曲线也叫钟形曲线，如下：

不同参数下的高斯曲线
当 $\mu=0,\sigma=1$ 时，高斯分布被称为标准正态分布。
scipy.special使用scipy.special.erf()计算高斯曲线的面积。

scipy.linalg

scipy.linalg.det():计算方阵的行列式
scipy.linalg.inv():计算方阵的逆
scipy.linalg.svd():奇异值分解

scipy.fftpack

快速傅立叶变换（FFT）,是快速计算序列的离散傅立叶变换（DFT）或其逆变换的方法。FFT会通过把DFT矩阵分解为稀疏因子之积来快速计算此类变换。

傅立叶变换将函数的时域与频域相关联

scipy.fftpack使用：

scipy.fftpack.fftfreq():生成样本序列
scipy.fftpack.fft():计算快速傅立叶变换

scipy.optimize

scipy.optimize模块提供了函数最值、曲线拟合和求根的算法。

函数最值
以寻找函数 $f(x)=x^2+10sin(x)$ 的最小值为例进行说明：
首先绘制目标函数的图形：

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

#定义目标函数
def f(x):
    return x**2+10*np.sin(x)

#绘制目标函数的图形
plt.figure(figsize=(10,5))
x = np.arange(-10,10,0.1)
plt.xlabel('x')
plt.ylabel('y')
plt.title('optimize')
plt.plot(x,f(x),'r-',label='$f(x)=x^2+10sin(x)$')
#图像中的最低点函数值
a = f(-1.3)
plt.annotate('min',xy=(-1.3,a),xytext=(3,40),arrowprops=dict(facecolor='black',shrink=0.05))
plt.legend()
plt.show()

图形输出如下：

先确定最小值所在区间

显然这是一个非凸优化问题，对于这类函数得最小值问题一般是从给定的初始值开始进行一个梯度下降，在optimize中一般使用bfgs算法。

optimize.fmin_bfgs(f,0)

运行结果

结果显示在经过五次迭代之后找到了一个局部最低点-7.945823，显然这并不是函数的全局最小值，只是该函数的一个局部最小值，这也是拟牛顿算法（BFGS）的局限性，如果一个函数有多个局部最小值，拟牛顿算法可能找到这些局部最小值而不是全局最小值，这取决与初始点的选取。在我们不知道全局最低点，并且使用一些临近点作为初始点，那将需要花费大量的时间来获得全局最优。此时可以采用暴力搜寻算法，它会评估范围网格内的每一个点。对于本例，如下：

grid = (-10, 10, 0.1)
xmin_global = optimize.brute(f, (grid,))
print(xmin_global)

搜寻结果如下：

全局最小值
但是当函数的定义域大到一定程度时，scipy.optimize.brute() 变得非常慢。scipy.optimize.anneal() 提供了一个解决思路，使用模拟退火算法。

可以使用scipy.optimize.fminbound(function,a,b)得到指定范围([a,b])内的局部最低点。

函数零点

scipy.optimize.fsolve(f,x):函数可以求解f=0的零点，x是根据函数图形特征预先估计的一个零点。

曲线拟合

scipy.optimize.curve_fit():非线性最小二乘拟合

from scipy import optimize

xdata = np.linspace(-10, 10, num=20)
ydata = f(xdata) + np.random.randn(xdata.size)
def f2(x, a, b):
    return a*x**2 + b*np.sin(x)

guess = [2, 2]
params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess)
print(params)

scipy.optimize.leatsq():最小二乘法拟合

'''
使用最小二乘法拟合直线
'''
import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt

#训练数据
Xi = np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
Yi = np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])

#定义拟合函数形式
def func(p,x):
    k,b = p
    return k*x+b

#定义误差函数
def error(p,x,y,s):
    print(s)
    return func(p,x)-y

#随机给出参数的初始值
p = [10,2]

#使用leastsq()函数进行参数估计
s = '参数估计次数'
Para = leastsq(error,p,args=(Xi,Yi,s))
k,b = Para[0]
print('k=',k,'\n','b=',b)

#图形可视化
plt.figure(figsize = (8,6))
#绘制训练数据的散点图
plt.scatter(Xi,Yi,color='r',label='Sample Point',linewidths = 3)
plt.xlabel('x')
plt.ylabel('y')
x = np.linspace(0,10,1000)
y = k*x+b
plt.plot(x,y,color= 'orange',label = 'Fitting Line',linewidth = 2)
plt.legend()
plt.show()

拟合效果如下：

最小二乘法拟合直线

'''
使用最小二乘法拟合正弦函数
'''
import numpy as np
from scipy.optimize import leastsq
import  matplotlib.pyplot as plt 

#定义拟合函数图形
def func(x,p):
    A,k,theta = p
    return A*np.sin(2*np.pi*k*x+theta)

#定义误差函数
def error(p,x,y):
    return y-func(x,p)

#生成训练数据
#随机给出参数的初始值
p0 = [10,0.34,np.pi/6]
A,k,theta = p0
x = np.linspace(0,2*np.pi,1000)
#随机指定参数

y0 = func(x,[A,k,theta])
#randn(m)从标准正态分布中返回m个值，在本例作为噪声
y1 = y0 + 2*np.random.randn(len(x))

#进行参数估计
Para = leastsq(error,p0,args=(x,y1))
A,k,theta = Para[0]
print('A=',A,'k=',k,'theta=',theta)


'''
图形可视化
'''
plt.figure(figsize=(20,8))
ax1 = plt.subplot(2,1,1)
ax2 = plt.subplot(2,1,2)

#在ax1区域绘图
plt.sca(ax1)
#绘制散点图
plt.scatter(x,y1,color='red',label='Sample Point',linewidth = 3)
plt.xlabel('x')
plt.xlabel('y')
y = func(x,p0)
plt.plot(x,y0,color='black',label='sine',linewidth=2)

#在ax2区域绘图
plt.sca(ax2)
e = y-y1
plt.plot(x,e,color='orange',label='error',linewidth=1)

#显示图例和图形
plt.legend()
plt.show()

最小二乘法拟合曲线
关于Scipy更多内容后续会慢慢进行更新，如有需要请参考：
Scipy:高级科学计算
Scipy官方文档