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

Scipy是一个用于数学、科学、工程领域的常用软件包，可以处理插值、积分、优化、图像处理、常微分方程数值解的求解、信号处理等问题。它用于有效计算Numpy矩阵，使Numpy和Scipy协同工作，高效解决问题。
Scipy是由针对特定任务的子模块组成：
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$ 函数，也称伽玛函数，是由欧拉积分定义的函数，是阶乘函数在实数域与复数域上的扩展，当 $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()函数进行计算。
贝塞尔函数

  在介绍贝塞尔函数之前，先对贝塞尔方程进行说明，一般来说，我们将形如

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

  高斯曲线是指高斯分布也就是我们常说的正态分布

其概率密度函数为：

其概率密度函数曲线就是高斯曲线也叫钟形曲线，如下：

#图像中的最低点函数值
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()
图形输出如下：
结果显示在经过五次迭代之后找到了一个局部最低点-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()
拟合效果如下：
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()