一  簡單介紹

SciPy是基於NumPy開發的高階模組，它提供了許多數學演算法和函式的實現，用於解決科學計算中的一些標準問題。例如數值積分和微分方程求解，擴充套件的矩陣計算，最優化，概率分佈和統計函式，甚至包括訊號處理等。
作為標準科學計算程式庫，SciPy類似於Matlab的工具箱，它是Python科學計算程式的核心包，它用於有效地計算NumPy矩陣，與NumPy矩陣協同工作。
SciPy庫由一些特定功能的子模組構成，如下表所示：

模組	功能
cluster	向量量化 / K-均值
constants	物理和數學常數
fftpack	傅立葉變換
integrate	積分程式
interpolate	插值
io	資料輸入輸出
linalg	線性代數程式
ndimage	n維影象包
odr	正交距離迴歸
optimize	優化
signal	訊號處理
sparse	稀疏矩陣
spatial	空間資料結構和演算法
special	任何特殊數學函式
stats	統計

以上子模組全依賴於NumPy且相互獨立，匯入NumPy和這些SciPy模組的標準方式如下，示例程式碼：

import numpy as np
from scipy import stats 

以上程式碼表示從SciPy模組中匯入stats子模組，SciPy的其他子模組匯入方式與之相同，限於機器學習研究領域及篇幅限制，本章將重點介紹linalg、optimize、interpolate及stats模組。

二 常用庫的介紹 2.1 線性代數linalg模組 linalg是Linear Algebra的縮寫，NumPy和SciPy都提供了線性代數函式庫linalg，SciPy的線性代數庫比NumPy更加全面。 （1）基本運算 linalg包含了許多方陣（包括矩陣）的基本運算函式，scipy.linalg.det()函式計算方陣的行列式，示例程式碼：

>>> from scipy import linalg
>>> arr = np.array([[1, 2], [3, 4]])
>>> linalg.det(arr)
-2.0
>>> arr = np.array([[3, 2],[6, 4]])
>>> linalg.det(arr) 
0.0
>>> linalg.det(np.ones((3, 4)))        #無論行列式還是逆矩陣只適用於n階矩陣的求解
Traceback (most recent call last):
...
ValueError: expected square matrix

scipy.linalg.inv()函式計算方陣的逆，示例程式碼：

>>> arr = np.array([[1, 2], [3, 4]])
>>> iarr = linalg.inv(arr)
>>> iarr
array([[-2. ,  1. ],
       [ 1.5, -0.5]])
>>>np.allclose(np.dot(arr, iarr), np.eye(2))      #numpy.allclose()函式用於比較兩方陣所有對應元素值，如果完全相同返回真(True)，否則返回假(False)
True

以下計算奇異陣（行列式為0）的逆，其結果將會報錯（LinAlgError），示例程式碼：

>>>arr = np.array([[3, 2], [6, 4]])
>>>linalg.inv(arr)
Traceback (most recent call last):
...
...LinAlgError: singular matrix

scipy.linalg.norm()函式計算方陣的範數，示例程式碼：

>>>A = np.matrix(np.random.random((2, 2)))
>>>A
>>>linalg.norm(A) #預設2範數
>>>linalg.norm(A, 1) #1範數
>>>linalg.norm(A, np.inf) #無窮範數

（2）解線性方程組

scipy.linalg.solve(A,b)和numpy.linalg.solve(A,b)可以用來解線性方程組Ax=b，即計算x=A^-1b。這裡，A是mm的方形矩陣，x和b是長為m的向量。有時候A是固定的，需要對多組b進行求解，因此第二個引數也可以是mn的矩陣B。這樣計算出來的X也是m*n的矩陣，相當於計算A^-1B。
在一些矩陣公式中經常會出現類似於A^-1B的運算，它們都可以用solve(A, B)計算，這要比直接逆矩陣然後做矩陣乘法更快捷一些，下面的程式比較solve()和逆矩陣的運算速度，示例程式碼：

>>> import numpy as np
>>> from scipy import linalg

>>> m, n = 500, 50
>>> A = np.random.rand(m, m)
>>> B = np.random.rand(m, n)
>>> X1 = linalg.solve(A, B)
>>> X2 = np.dot(linalg.inv(A), B)
>>> print(np.allclose(X1, X2))

>>> %timeit linalg.solve(A, B)
13.3 ms ± 834 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

>>> %timeit np.dot(linalg.inv(A), B)
22.4 ms ± 1.48 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

(3) 特徵值和特徵向量

n*n的矩陣A可以看作n維空間中的線性變換。若x為n維空間中的一個向量，那麼A與x的矩陣乘積就是對x進行線性變換之後的向量。如果x是線性變換的特徵向量，那麼經過這個線性變換之後，得到的新向量仍然與原來的x保持在同一方向上，但其長度也許會改變。特徵向量的長度在該線性變換下縮放的比例稱為特徵值。即特徵向量x滿足如下等式，λ的值可以是一個任意複數：Ax=λx。
下面以二維平面上的線性變換矩陣為例，演示特徵值和特徵向量的幾何含義。通過linalg.eig(A)計算矩陣A的兩個特徵值evalues和特徵向量evectors，在evectors中，每一列是一個特徵向量。示例程式碼：

>>> A = np.array([[1, -0.3], [-0.1, 0.9]])
>>> evalues, evectors = linalg.eig(A)

2.2 擬合與求解optimize模組

SciPy的optimize模組提供了許多數值優化的演算法，一些經典的優化演算法包括線性迴歸、函式極值和根的求解以及確定兩函式交點的座標等。下面首先介紹簡單的線性迴歸模型，然後逐漸深入解決非線性資料擬合問題。
（1）擬合 curve_fit()函式 線性迴歸有許多擬合數據的方法，我們將使用curve_fit()函式，它利用的是最小二乘演算法。最小二乘演算法是一種數學優化技術，在機器學習領域最有名和有效的演算法之一。它通過最小化誤差的平方和尋找資料的最佳函式匹配。利用最小二乘法可以簡便地求得未知的資料，並使得這些求得的資料與實際資料之間誤差的平方和為最小。
以下示例中，我們首先從已知函式中生成一些帶有噪聲的資料，然後使用curve_fit()函式擬合這些噪聲資料。示例中的已知函式我們使用一個簡單的線性方程式，即f(x)=ax+b。示例程式碼：

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

#建立函式模型用來生成資料
def func(x, a, b):
       return a*x + b

#生成乾淨資料
x = np.linspace(0, 10, 100)
y = func(x, 1, 2)

#對原始資料新增噪聲
yn = y + 0.9 * np.random.normal(size=len(x))

#使用curve_fit函式擬合噪聲資料
popt, pcov = curve_fit(func, x, yn)

#輸出給定函式模型func的最優引數
print(popt)

結果為：

[ 0.99734363  1.96064258]

如果有一個很好的擬合效果，popt返回的是給定模型的最優引數。我們可以使用pcov的值檢測擬合的質量，其對角線元素值代表著每個引數的方差。

>>>print(pcov)
[[ 0.00105056 -0.00525282]
 [-0.00525282  0.03519569]]

通過以下程式碼繪製出了擬合曲線與實際曲線的差異，示例程式碼：

yfit = func(x,popt[0],popt[1]) 

plt.plot(x, y, color="green",label = "actual data")
plt.plot(x, yn, "o", label = "actual data with noise")
plt.plot(x, yfit,color="yellow", label = "fitting data")
plt.legend(loc = "best")
plt.show()

下面做進一步研究，我們可以通過最小二乘擬合高斯分佈（Gaussian profile），一種非線性函式：α*exp(-(x-μ)²/2σ²)
在這裡，α表示一個標量，μ是期望值，而σ是標準差。示例程式碼：

import numpy as np 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plt
#建立一個函式模型用來生成資料
def func(x, a, b, c):
        return (a*np.exp(-(x-b)**2/2*c**2))

#生成原始資料
x = np.linspace(0, 10, 100)
y = func(x, 1, 5, 2)

#對原始資料增加噪聲
yn = y + 0.2*np.random.normal(size=len(x))

#使用curve_fit函式擬合噪聲資料
popt, pcov = curve_fit(func, x, yn)

#popt返回最擬合給定的函式模型func的引數值，如popt[0]=a,popt[1]=b,popt[2]=3
print(popt)

結果為：

[-0.49627942  2.78765808 28.76127826]

通過以下程式碼繪製出了擬合曲線與實際曲線的差異，示例程式碼：

p0=[1.2,4,3] #初步猜測引數，如果沒有，預設全為1，即[1,1,1]
popt, pcov = curve_fit(func, x, yn,p0=p0)

#popt返回最擬合給定的函式模型func的引數值，如popt[0]=a,popt[1]=b,popt[2]=3
print(popt)

yfit = func(x,popt[0],popt[1],popt[2])

plt.plot(x, y, color="green",label = "actual data")
plt.plot(x, yn, "o", label = "actual data with noise")
plt.plot(x, yfit, color="yellow", label = "fitting data")
plt.legend(loc = "best")
plt.show()

結果如下圖所示：

    通過以上繪圖，我們可以看出對高斯分佈函式擬合的效果是可以接受的。
隨著研究的深入，我們可以擬合一個多重高斯分佈的一維資料集。現在將這個函式擴充套件為包含兩個不同輸入值的高斯分佈函式。這是一個擬合線性光譜的經典例項，示例程式碼如下：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def func(x, a0, b0, c0, a1, b1, c1):
      return a0*np.exp(-(x - b0) ** 2/(2 * c0 ** 2)) + a1 * np.exp(-(x-b1) ** 2/(2 * c1 ** 2))

#生成原始資料
x = np.linspace(0, 20, 200)
y = func(x, 1, 3, 1, -2, 15, 0.5)

#對原始資料增加噪聲
yn = y + 0.9 * np.random.normal(size=len(x))

#如果要擬合一個更加複雜的函式，提供一些估值假設對擬合效果更好
guesses = [1, 3, 1, 1, 15, 1]

#使用curve_fit函式擬合噪聲資料
popt, pcov = curve_fit(func, x, yn, p0=guesses)

#popt返回最擬合給定的函式模型func的引數值，如popt[0]=a,popt[1]=b,popt[2]=3
print(popt)

yfit = func(x,popt[0],popt[1],popt[2],popt[3],popt[4],popt[5])

plt.plot(x, y, color="green",label = "actual data")
plt.plot(x, yn, "o", label = "actual data with noise")
plt.plot(x, yfit, color="yellow", label = "fitting data")
plt.legend(loc = "best")
plt.show()

結果如下圖所示：

（2）最小二乘擬合leastsq()函式 假設有一組實驗資料(x[i], y[i])，我們知道它們之間的函式關係:y = f(x)，通過這些已知資訊，需要確定函式中的一些引數項。例如，如果f是一個線型函式f(x) = k*x+b，那麼引數k和b就是我們需要確定的值。如果將這些引數用 p 表示的話，那麼我們就是要找到一組 p 值使得如下公式中的S函式最小： 這種演算法被稱之為最小二乘擬合(Least-square fitting)。optimize模組提供了實現最小二乘擬合算法的函式leastsq()，leastsq是least square的簡寫，即最小二乘法。 下面是用leastsq()對線性函式進行擬合的程式，示例程式碼：

import matplotlib.pylab as plt
import numpy as np
from scipy import optimize    # 從scipy庫引入optimize模組

X = np.array([ 8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78 ])
Y = np.array([ 7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05 ])

def residuals(p):
        #計算以p為引數的直線和原始資料之間的誤差
        k, b = p
        return Y-(k*X+b)

# leastsq()使得residuals()的輸出陣列的平方和最小，引數的初始值為[1, 0]
r = optimize.leastsq(residuals, [1,0])
k, b = r[0]
print("k=", k, "b=", b)

結果為：

k = 0.613495349193  b = 1.79409254326

可以通過通過繪圖對比真實資料和擬合數據的誤差，示例程式碼；

plt.plot(X, Y, "o", label = "actual data")
plt.plot(X, k*X+b, label = "fitting data")
plt.legend(loc = "best")
plt.show()

結果為：

繪圖中的圓點表示真實資料點，實線表示擬合曲線，由此看出擬合引數得到的函式和真實資料大體一致。 接下來，用leastsq()對正弦波資料進行擬合，示例程式碼：

import numpy as np
from scipy.optimize import leastsq   # 從scipy庫的optimize模組引入leastsq函式
import matplotlib.pyplot as plt    # 引入繪圖模組pylab，並重命名為pl

def func(x, p):
    """
    資料擬合所用的函式: A*sin(2*pi*k*x + theta)
    """
    A, k, theta = p
    return A*np.sin(2*np.pi*k*x+theta)   

def residuals(p, y, x):
    """
    實驗資料x, y和擬合函式之間的差，p為擬合需要找到的係數
    """
    return y - func(x, p) 

x = np.linspace(0, -2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6   # 真實資料的函式引數
y0 = func(x, [A, k, theta])   # 真實資料

y1 = y0 + 2 * np.random.randn(len(x))   # 加入噪聲之後的實驗資料，噪聲是服從標準正態分佈的隨機量    

p0 = [7, 0.2, 0]   # 第一次猜測的函式擬合引數

# 呼叫leastsq進行資料擬合
# residuals為計算誤差的函式
# p0為擬合引數的初始值
# args為需要擬合的實驗資料
plsq = leastsq(residuals, p0, args=(y1, x))

print ("actual parameter:", [A, k, theta]) # 真實引數
print ("fitting parameter", plsq[0]) # 實驗資料擬合後的引數

plt.plot(x, y0, label="actual data") # 繪製真實資料
plt.plot(x, y1, label="experimental data with noise")  # 帶噪聲的實驗資料
plt.plot(x, func(x, plsq[0]), label="fitting data")    # 擬合數據
plt.legend()
plt.show()

這個例子中我們要擬合的函式是一個正弦波函式，它有三個引數 A, k, theta ，分別對應振幅、頻率、相角。假設我們的實驗資料是一組包含噪聲的資料 x, y1，其中y1是在真實資料y0的基礎上加入噪聲的到了。 通過leastsq函式對帶噪聲的實驗資料x, y1進行資料擬合，可以找到x和真實資料y0之間的正弦關係的三個引數： A, k, theta。下面是程式的輸出：

>>>actual parameter: [10, 0.34, 0.5235987755982988]
>>>fitting parameter [ 10.12646889   0.33767587   0.48944317]

我們看到擬合引數雖然和真實引數完全不同，但是由於正弦函式具有周期性，實際上擬合引數得到的函式和真實引數對應的函式是一致的。
（3）標量函式極值求解fmin()函式
首先定義以下函式，然後繪製它，示例程式碼：

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
def f(x):
    return x**2 + 10*np.sin(x)  
x = np.arange(-10, 10, 0.1)
plt.plot(x, f(x)) 
plt.show()

結果如下圖所示：

如圖所示，該函式大約在-1.3有個全域性最小值，在3.8有個區域性最小值。找到這個函式最小值一般而有效的方法是從初始點使用梯度下降法。BFGS演算法是做這個的好方法，BFGS演算法被認為是數值效果最好的擬牛頓法，是由Broyden，Fletcher，Goldfarb，Shanno四個人分別提出的，故稱為BFGS校正。具體演算法思想及解釋請查閱相關資料，這裡直接通過optimize.fmin_bfgs()函式求解最小值，示例程式碼：

>>> optimize.fmin_bfgs(f, 0)
Optimization terminated successfully.
         Current function value: -7.945823
         Iterations: 5
         Function evaluations: 24
         Gradient evaluations: 8
array([-1.30644003])

這個方法一個可能的問題在於，如果函式有區域性最小值，演算法會因初始點不同找到這些區域性最小而不是全域性最小，示例程式碼：

>>> optimize.fmin_bfgs(f, 3, disp=0)#disp是布林型資料，如果為1，列印收斂訊息
array([ 3.83746663])

如果我們不知道全域性最小值的鄰近值來選定初始點，我們需要藉助於耗費資源些的全域性優化。為了找到全域性最小點，最簡單的演算法是蠻力演算法，該演算法求出給定格點的每個函式值。示例程式碼：

>>>grid = (-10, 10, 0.1)
>>>xmin_global = optimize.brute(f, (grid, ))
>>>xmin_global
array([-1.30641113])

對於大點的格點，scipy.optimize.brute()變得非常慢。scipy.optimize.anneal()提供了使用模擬退火的替代函式。對已知的不同類別全域性優化問題存在更有效率的演算法，但這已經超出scipy的範圍。為了找到區域性最小，我們把變數限制在(0,10)之間，使用scipy.optimize.fminbound()，示例程式碼：

>>> xmin_local = optimize.fminbound(f, 0, 10)
>>> xmin_local
3.8374671...

下面的程式通過求解卷積的逆運算演示fmin的功能。對於一個離散線性時不變系統h, 如果輸入是x，那麼其輸出y可以用x和h的卷積表示：

現在的問題是如果已知系統的輸入x和輸出y，如何計算系統的傳遞函式h；或者如果已知系統的傳遞函式h和系統的輸出y，如何計算系統的輸入x。這種運算被稱為反捲積運算，是十分困難的，特別是在實際的運用中，測量系統的輸出總是存在誤差的。 下面用fmin計算反捲積，這種方法只能用在很小規模的數列之上，因此沒有很大的實用價值，不過用來評價fmin函式的效能還是不錯的。示例程式碼：

import scipy.optimize as opt 
import numpy as np 

def test_fmin_convolve(fminfunc, x, h, y, yn, x0): 
    """
    x (*) h = y, (*)表示卷積
    yn為在y的基礎上新增一些干擾噪聲的結果
    x0為求解x的初始值
    """
    def convolve_func(h):
        """
        計算 yn - x (*) h 的power
        fmin將通過計算使得此power最小
        """ 
        return np.sum((yn - np.convolve(x, h))**2) 

    # 呼叫fmin函式，以x0為初始值
    h0 = fminfunc(convolve_func, x0) 

    print fminfunc.__name__ 
    print "---------------------" 
    # 輸出 x (*) h0 和 y 之間的相對誤差
    print "error of y:", np.sum((np.convolve(x, h0)-y)**2)/np.sum(y**2) 
    # 輸出 h0 和 h 之間的相對誤差
    print "error of h:", np.sum((h0-h)**2)/np.sum(h**2) 
    print 

def test_n(m, n, nscale): 
    """
    隨機產生x, h, y, yn, x0等數列，呼叫各種fmin函式求解b
    m為x的長度, n為h的長度, nscale為干擾的強度
    """
    x = np.random.rand(m) 
    h = np.random.rand(n) 
    y = np.convolve(x, h) 
    yn = y + np.random.rand(len(y)) * nscale
    x0 = np.random.rand(n) 

    test_fmin_convolve(opt.fmin, x, h, y, yn, x0) 
    test_fmin_convolve(opt.fmin_powell, x, h, y, yn, x0) 
    test_fmin_convolve(opt.fmin_cg, x, h, y, yn, x0)
    test_fmin_convolve(opt.fmin_bfgs, x, h, y, yn, x0)

test_n(200, 20, 0.1)

程式碼

執行結果為：

fmin
---------------------
error of y: 0.000360456186137
error of h: 0.0122264525455
Optimization terminated successfully.
         Current function value: 0.207509
         Iterations: 96
         Function evaluations: 17400
fmin_powell
---------------------
error of y: 0.000129249083036
error of h: 0.000300953639205
Optimization terminated successfully.
         Current function value: 0.207291
         Iterations: 20
         Function evaluations: 880
         Gradient evaluations: 40
fmin_cg
---------------------
error of y: 0.000129697740414
error of h: 0.000292820536053
Optimization terminated successfully.
         Current function value: 0.207291
         Iterations: 31
         Function evaluations: 946
         Gradient evaluations: 43
fmin_bfgs
---------------------
error of y: 0.000129697643272
error of h: 0.000292817401206

結果

（4）函式求解fsolve()

optimize庫中的fsolve函式可以用來對非線性方程組進行求解，其基本呼叫形式是                                                                               fsolve(func, x0)

• func是用於定義需求解的非線性方程組的函式檔名
• x0為未知數向量的初始值。

首先通過一個簡單的示例，利用fsolve()函式求解當線性函式為0時，x的值，示例程式碼：

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

line = lambda x:x+3

solution = fsolve(line, -2)
print(solution)

結果為：

[-3，]

通過以下繪圖函式可以看出當函式等於0時，x軸的座標值為-3，示例程式碼：

x = np.linspace(-5.0, 0, 100)
plt.plot(x,line(x), color="green",label = "function")
plt.plot(solution,line(solution), "o", label = "root")
plt.legend(loc = "best")
plt.show()

結果為：

下面我們通過一個簡單的示例介紹一下兩個方程交點的求解方法，示例程式碼：

from scipy.optimize import fsolve
import numpy as np
import matplotlib.pyplot as plt
# 定義解函式
def findIntersection(func1, func2, x0):
        return fsolve(lambda x: func1(x)-func2(x),x0)

# 定義兩方程
funky = lambda x : np.cos(x / 5) * np.sin(x / 2)
line = lambda x : 0.01 * x - 0.5

# 定義兩方程交點的取值範圍
x = np.linspace(0, 45, 1000)
result = findIntersection(funky, line, [15, 20, 30, 35, 40, 45])

# 輸出結果
print(result, line(result))


plt.plot(x,funky(x), color="green",label = "funky func")
plt.plot(x,line(x), color="yellow",label = "line func")
plt.plot(result,line(result), "o", label = "intersection")
plt.legend(loc = "best")
plt.show()

結果為：

如果要對如下方程組進行求解的話：

• f1(u1,u2,u3) = 0
• f2(u1,u2,u3) = 0
• f3(u1,u2,u3) = 0

那麼func可以如下定義：

def func(x):
    u1,u2,u3 = x
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

下面是一個實際的例子，求解如下方程組的解：

• 5*x1 + 3 = 0
• 4*x0*x0 - 2*sin(x1*x2) = 0
• x1*x2 - 1.5 = 0

示例程式碼：

from scipy.optimize import fsolve
from math import sin,cos

def f(x):
    x0 = float(x[0])
    x1 = float(x[1])
    x2 = float(x[2])
    return [
        5*x1+3,
        4*x0*x0 - 2*sin(x1*x2),
        x1*x2 - 1.5
    ]
result = fsolve(f, [1,1,1])
print (result)

結果為：

[-0.70622057 -0.6        -2.5       ]

2.3 插值interpolate模組

插值是離散函式逼近的重要方法，利用它可通過函式在有限個點處的取值狀況，估算出函式在其他點處的近似值。與擬合不同的是，要求曲線通過所有的已知資料。SciPy的interpolate模組提供了許多對資料進行插值運算的函式，範圍涵蓋簡單的一維插值到複雜多維插值求解。當樣本資料變化歸因於一個獨立的變數時，就使用一維插值；反之樣本資料歸因於多個獨立變數時，使用多維插值。
計算插值有兩種基本的方法，1、對一個完整的資料集去擬合一個函式；2、對資料集的不同部分擬合出不同的函式，而函式之間的曲線平滑對接。第二種方法又叫做仿樣內插法，當資料擬合函式形式非常複雜時，這是一種非常強大的工具。我們首先介紹怎樣對簡單函式進行一維插值運算，然後進一步深入比較複雜的多維插值運算。 （1）一維插值 一維資料的插值運算可以通過函式interp1d()完成。其呼叫形式如下，它實際上不是函式而是一個類：

interp1d(x, y, kind='linear', ...)

其中，x和y引數是一系列已知的資料點，kind引數是插值型別，可以是字串或整數，它給出插值的B樣條曲線的階數，候選值及作用下表所示：

候選值	作用
‘zero’ 、'nearest'	階梯插值，相當於0階B樣條曲線
‘slinear’ 、'linear'	線性插值，用一條直線連線所有的取樣點，相當於一階B樣條曲線
‘quadratic’ 、'cubic'	二階和三階B樣條曲線，更高階的曲線可以直接使用整數值指定

下面的程式演示了通過不同的 kind引數（linear和quadratic），對一個正弦函式進行插值運算。示例程式碼：

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

#建立待插值的資料
x = np.linspace(0, 10*np.pi, 20)
y = np.cos(x)

# 分別用linear和quadratic插值
fl = interp1d(x, y, kind='linear')
fq = interp1d(x, y, kind='quadratic')

#設定x的最大值和最小值以防止插值資料越界
xint = np.linspace(x.min(), x.max(), 1000)
yintl = fl(xint)
yintq = fq(xint)


plt.plot(xint,fl(xint), color="green", label = "Linear")
plt.plot(xint,fq(xint), color="yellow", label ="Quadratic")
plt.legend(loc = "best")
plt.show()

結果如下圖所示：

（2）噪聲資料插值  我們可以通過interpolate模組中UnivariateSpline()函式對含有噪聲的資料進行插值運算，示例程式碼：

import numpy as np
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt

# 通過人工方式新增噪聲資料
sample = 30
x = np.linspace(1, 10*np.pi, sample)
y = np.cos(x) + np.log10(x) + np.random.randn(sample)/10

# 插值，引數s為smoothing factor 
f = UnivariateSpline(x, y, s=1)
xint = np.linspace(x.min(), x.max(), 1000)
yint = f(xint)

plt.plot(xint,f(xint), color="green", label = "Interpolation")
plt.plot(x, y, color="yellow", label ="Original")
plt.legend(loc = "best")
plt.show()

需要說明的是：在UnivariateSpline()函式中，引數s是平滑向量引數，被用來擬合還有噪聲的資料。如果引數s=0，將忽略噪聲對所有點進行插值運算。結果如下圖所示：

  （3）多維插值  多維插值主要用於重構圖片，interpolate模組中的griddata()函式有很強大的處理多維雜湊取樣點進行插值運算的能力。其呼叫形式如下：

griddata(points, values, xi, method='linear', fill_value=nan)

其中points表示K維空間中的座標，它可以是形狀為(N,k)的陣列，也可以是一個有k個數組的序列，N為資料的點數。values是points中每個點對應的值。xi是需要進行插值運算的座標，其形狀為(M,k)。method引數有三個選項：'nearest'、 ‘linear’、 'cubic'，分別對應0階、1階以及3階插值。以下示例利用1000個隨機雜湊點對1000x1000畫素的圖片進行重構，示例程式碼：

import numpy as np
from scipy.interpolate import griddata#定義一個函式
def ripple(x,y):
    return np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

# 生成grid資料，複數定義了生成grid資料的step，若無該複數則step為5 
grid_x, grid_y = np.mgrid[0:5:1000j, 0:5:1000j] 

# 生成待插值的樣本資料 
points = np.random.rand(1000,2) 

value = ripple(points[:,0]*5,points[:,1]*5) 

# 用nearest方法插值
grid_z0 = griddata(points*5,value, (grid_x,grid_y),method='nearest')

我們還可以使用interpolate模組的SmoothBivariateSpline類進行多元仿樣插值運算，對圖片進行重構。示例程式碼：

import numpy as np
from scipy.interpolate import SmoothBivariateSpline as SBS

#定義一個函式
def ripple(x,y):
    return np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

# 生成grid資料，複數定義了生成grid資料的step，若無該複數則step為5 
grid_x, grid_y = np.mgrid[0:5:1000j, 0:5:1000j] 

# 生成待插值的樣本資料 
points = np.random.rand(1000,2)
 
value = ripple(points[:,0]*5,points[:,1]*5) 

# 用nearest方法插值
fit = SBS(points[:,0]*5, points[:,1]*5, value, s=0.01, kx=4, ky=4)
interp = fit(np.linspace(0, 5, 1000), np.linspace(0, 5, 1000))

我們得到了一個與上個示例同樣的結果。整體上SmoothBivariateSpline函式的表現略好於griddata函式。
通過反覆測試，儘管SmoothBivariateSpline表現略好，但其對給定的樣本資料非常敏感，就可能導致忽略一些顯著特徵。而griddata函式有很強的魯棒性，不管給定的資料樣本，能夠合理的進行插值運算。 2.4 統計stats模組
NumPy庫已經提供了一些基本的統計函式，如求期望、方差、中位數、最大值和最小值等。示例程式碼：

import numpy as np

#構建一個1000個隨機變數的陣列
x = np.random.randn(1000)

#對陣列元素的值進行統計
mean = x.mean()
std = x.std()
var = x.var()

print(mean,std,var)

結果為：

(0.02877273942510088, 0.97623362287515114, 0.95303208643194282)

mean是期望值，std是標準差，var是方差，使用numpy.array物件已有的方法獲得統計指標快速有效，而SciPy庫則提供了更高階的統計工具，它的Stats模組包含了多種概率分佈的隨機變數（隨機變數是指概率論中的概念，不是Python中的變數），其中隨機變數又分為連續和離散兩種。所有的連續隨機變數都是rv_continuous的派生類的物件，而所有的離散隨機變數都是rv_discrete的派生類的物件。 （1）連續概率分佈
SciPy的stats模組提供了大約80種連續隨機變數和10多種離散分佈變數，這些分佈都依賴於numpy.random函式。可以通過如下語句獲得stats模組中所有的連續隨機變數，示例程式碼：

from scipy import stats
[k for k, v in stats.__dict__.items() if isinstance(v, stats.rv_continuous)]

執行結果為：

['ksone', 'kstwobign', 'norm', 'alpha', 'anglit', 'arcsine', 'beta', 'betaprime', 'bradford', 'burr', 'burr12', 'fisk', 'cauchy', 'chi', 'chi2', 'cosine', 'dgamma', 'dweibull', 'expon', 'exponnorm', 'exponweib', 'exponpow', 'fatiguelife', 'foldcauchy', 'f', 'foldnorm', 'frechet_r', 'weibull_min', 'frechet_l', 'weibull_max', 'genlogistic', 'genpareto', 'genexpon', 'genextreme', 'gamma', 'erlang', 'gengamma', 'genhalflogistic', 'gompertz', 'gumbel_r', 'gumbel_l', 'halfcauchy', 'halflogistic', 'halfnorm', 'hypsecant', 'gausshyper', 'invgamma', 'invgauss', 'invweibull', 'johnsonsb', 'johnsonsu', 'laplace', 'levy', 'levy_l', 'levy_stable', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'gilbrat', 'maxwell', 'mielke', 'kappa4', 'kappa3', 'nakagami', 'ncx2', 'ncf', 't', 'nct', 'pareto', 'lomax', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rdist', 'rayleigh', 'reciprocal', 'rice', 'recipinvgauss', 'semicircular', 'skewnorm', 'trapz', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'vonmises_line', 'wald', 'wrapcauchy', 'gennorm', 'halfgennorm']

連續隨機變數物件主要使用如下方法，下表所示：

方法名	全稱	功能
rvs	Random Variates of given type	對隨機變數進行隨機取值，通過size引數指定輸出陣列的大小
pdf	Probability Density Function	隨機變數的概率密度函式
cdf	Cumulative Distribution Function	隨機變數的累積分佈函式，它是概率密度函式的積分
sf	Survival function	隨機變數的生存函式，它的值是1-cdf(t)
ppf	Percent point function	累積分佈函式的反函式
stats	statistics	計算隨機變數的期望值和方差
fit	fit	對一組隨機取樣進行擬合，找出最適合取樣資料的概率密度函式的係數

下面以標準正態分佈（函式表示f(x)=(1/√2π)exp(-x^2/2)）為例，簡單介紹隨機變數的用法。示例程式碼：

from scipy import stats
# 設定正態分佈引數，其中loc是期望值引數，scale是標準差引數
X = stats.norm(loc=1.0, scale=2.0)

# 計算隨機變數的期望值和方差
print(X.stats())

結果為：

(array(1.0), array(4.0))

以上程式碼說明，norm可以像函式一樣呼叫，通過loc和scale引數可以指定隨機變數的偏移和縮放參數。對於正態分佈的隨機變數來說，這兩個引數相當於指定其期望值和標準差，標準差是方差的算術平方根。X的stats()方法，可以計算隨機變數X分佈的特徵值，如期望值和方差。
此外，通過呼叫隨機變數X的rvs()方法，可以得到包含一萬次隨機取樣值的陣列x，然後呼叫NumPy的mean()和var()計算此陣列的均值和方差，其結果符合隨機變數X的特性，示例程式碼：

#對隨機變數取10000個值
x = X.rvs(size=10000)
print(np.mean(x), np.var(x))

結果為：

(1.0287787687588861, 3.9944276709242805)

使用fit()方法對隨機取樣序列x進行擬合，它返回的是與隨機取樣值最吻合的隨機變數引數，示例程式碼：

#輸出隨機序列的期望值和標準差
print(stats.norm.fit(x))

結果為：

(1.0287787687588861, 1.998606432223283)

在下面的例子中，計算取樣值x的直方圖統計以及累計分佈，並與隨機變數的概率密度函式和累積分佈函式進行比較。示例程式碼：

pdf, t = np.histogram(x, bins=100, normed=True)
t = (t[:-1]+t[1:])*0.5
cdf = np.cumsum(pdf) * (t[1] - t[0])
p_error = pdf - X.pdf(t)
c_error = cdf - X.cdf(t)
print("max pdf error: {}, max cdf error: {}".format(np.abs(p_error).max(), np.abs(c_error).max()))

執行結果如下所示：

max pdf error: 0.0208405611169, max cdf error: 0.0126874590568

通過繪圖的方式檢視概率密度函式求得的理論值（theory value）和直方圖統計值（statistic value），可以看出二者是一致的，示例程式碼：

import pylab as pl
pl.plot(t, pdf, color="green", label = "statistic value")
pl.plot(t, X.pdf(t), color="yellow", label ="theory value")
pl.legend(loc = "best")
pl.show()

結果見下圖所示：

也可以用同樣的方式顯示隨機變數X的累積分佈和陣列pdf的累加結果，示例程式碼：

import pylab as pl
pl.plot(t, cdf, color="green", label = "statistic value")
pl.plot(t, X.cdf(t), color="yellow", label ="theory value")
pl.legend(loc = "best")
pl.show()

結果為：

（2）離散概率分佈

# 陣列x儲存骰子的所有可能值，陣列p儲存每個值出現的概率
x = range(1, 7)
p = (0.4, 0.2, 0.1, 0.1, 0.1, 0.1)

# 建立表示這個骰子的隨機變數dice，呼叫其rvs()方法投擲此骰子20次，獲得符合概率p的隨機數
dice = stats.rv_discrete(values=(x, p))
print(dice.rvs(size=20))

執行結果：

array([3, 6, 4, 5, 5, 2, 1, 3, 3, 1, 1, 3, 1, 5, 1, 3, 4, 1, 2, 2])

除了自定義離散概率分佈，我們也可以利用stats模組裡的函式定義各種分佈。下面以生成幾何分佈為例，其函式是geom()，示例程式碼：

import numpy as np
from scipy.stats import geom

# 設定幾何分佈的引數
p = 0.5
dist = geom(p)

# 設定樣本區間  
x = np.linspace(0, 5, 1000)  

# 得到幾何分佈的 PMF 和CDF  
pmf = dist.pmf(x) 
cdf = dist.cdf(x)  

# 生成500個隨機數  
sample = dist.rvs(500)

（3）描述與檢驗函式

 SciPy中有超過60種統計函式。stats模組包括了諸如kstest 和normaltest等樣本測試函式，用來檢測樣本是否服從某種分佈。在使用這些工具前，要對資料有較好的理解，否則可能會誤讀它們的結果。樣本分佈檢驗為例，示例程式碼：

import numpy as np 
from scipy import stats 

# 生成包括100個服從正態分佈的隨機數樣本
sample = np.random.randn(100) 

# 用normaltest檢驗原假設
out = stats.normaltest(sample) 
print('normaltest output') 
print('Z-score = ' + str(out[0])) 
print('P-value = ' + str(out[1])) 

# kstest 是檢驗擬合度的Kolmogorov-Smirnov檢驗，這裡針對正態分佈進行檢驗
# D是KS統計量的值，越接近0越好
out = stats.kstest(sample, 'norm') 
print('\nkstest output for the Normal distribution') 
print('D = ' + str(out[0])) 
print('P-value = ' + st