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

简介

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

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函数可以定义为：

  \displaystyle \gamma(z)=\int_0^{+\infty}\frac{t^{z-1}}{e^t}dt

  由「解析延拓原理」可以将这个定义拓展到整个复数域上，非正整数除外。
  在scipy.special中使用scipy.special.gamma()实现\gamma函数的计算，如果对于精度有更高要求，可以使用采用对数坐标的scipy.special.gammaln()函数进行计算。

• 贝塞尔函数
  一般来说，将形如

  \displaystyle x^2\frac{d^2y}{dx^2}+x\frac{dy}{dx}+(x^2-\alpha^2)y=0

  的二阶线性常微分方程称为贝塞尔方程，二贝塞尔方程的标准解函数y(x)就是贝塞尔函数。其中，参数\alpha被称为对应贝塞尔函数的阶。贝塞尔方程是一个二阶线性齐次常微分方程，必然存在两个线性无关的解，然而通常情况下它的解无法用初等函数表示，但是当\alpha=n时，注意到x=0是贝塞尔方程的「正则奇点」，则由常微分方程的广义幂级数解法可以得出贝塞尔方程的两个广义幂级数解：

  \displaystyle y_1=J_n(x)=\sum_{k=0}^\infty\frac{(-1)^k}{\gamma(n+k+1)\gamma(k+1)}\left(\frac{x}{2}\right)^{2k+n}\\{}\\ y_2=J_{-n}(x)=\sum_{k=0}^\infty\frac{(-1)^k}{\gamma(-n+k+1)\gamma(k+1)}\left(\frac{x}{2}\right)^{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 \mathrm{Z}\}

构成f(z)的全部周期。
在scipy.special中使用scipy.special.ellipj()函数计算椭圆函数。

• Erf(高斯曲线的面积)

高斯曲线指的是高斯分布（正态分布）

X\sim N(\mu,\sigma^2)

其概率密度为：

\displaystyle 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+10\sin(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()

本作品采用《CC 协议》，转载必须注明作者和本文链接