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 協議》，轉載必須註明作者和本文連結