首先，是通常的導入工作：
In [1]: import numpy as np
from pylab import plt, mpl
In [2]: plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
%matplotlib inline 本節使用的主函數示例如下，由一個三角函數項和一個線性項組成：
In [3]: def f(x):
return np.sin(x) + 0.5 * x 重點是在給定區間内通過回歸和插值求取該函數的近似值。首先，生成該函數的圖形，以便更好地觀察逼近法的效果。我們感興趣的區間是[−2π，2π]。圖11-1顯示了該函數在np.linspace()函數定義的固定區間上的圖像。create_plot()是一個助手函數，可以建立本章多次要使用的同類圖表：
圖11-1 示例函數圖表
In [4]: def create_plot(x, y, styles, labels, axlabels):
plt.figure(figsize=(10, 6))
for i in range(len(x)):
plt.plot(x[i], y[i], styles[i], label=labels[i])
plt.xlabel(axlabels[0])
plt.ylabel(axlabels[1])
plt.legend(loc=0)
In [5]: x = np.linspace(-2 * np.pi, 2 * np.pi, 50) ❶
In [6]: create_plot([x], [f(x)], ['b'], ['f(x)'], ['x', 'f(x)']) ❶ 用于繪圖和計算的x值。
11.1.1 回歸
回歸是相當高效的函數近似值計算工具。它不僅适用于求取一維函數的近似值，在更高次元上也很有效。得出回歸結果所需要的數值化方法很容易實作，執行也很快速。本質上，回歸的任務是在給定一組所謂“基函數”bd，d∈{1，…，D}的情況下，根據公式11-1找出最優參數
，…
，其中對于i∈{1,…I}觀察點，yi ≡ f(xi)。xi可以視為自變量觀測值，yi可視為因變量觀測值（從函數或者統計的意義上說）。
公式11-1. 最小化回歸問題
1．作為基函數的單項式
最簡單的情況是以單項式作為基函數——也就是說，b1=1，b2=x，b3=x2，b4=x3，…在這種情況下，NumPy有可以确定最優參數（np.polyfit()）和通過一組輸入值求取近似值（np.polyval()）的内建函數。
表11-1列出了np.polyfit()函數的參數。在np.polyfit()傳回的最優回歸相關系數ρ基礎上，np.polyval(ρ,x)傳回x坐标的回歸值。
典型向量化風格的np.polyfit()和np.polyval()線性回歸（deg=1）的應用方式如下。由于回歸估算值儲存在ry數組中，是以我們可以如圖11-2那樣比較回歸結果和原始函數。當然，線性回歸無法處理示例函數的sin部分：
In [7]: res = np.polyfit(x, f(x), deg=1, full=True) ❶
In [8]: res ❷
Out[8]: (array([ 4.28841952e-01, -1.31499950e-16]),
array([21.03238686]),
2,
array([1., 1.]),
1.1102230246251565e-14)
In [9]: ry = np.polyval(res[0], x) ❸
In [10]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
['f(x)', 'regression'], ['x', 'f(x)']) ❶ 線性回歸步驟。
❷ 完整的結果：回歸參數、殘差、有效秩、奇異值和相對條件數。
❸ 使用回歸參數求值。
圖11-2 線性回歸
為了處理示例函數的sin部分，必須使用更高次的單項式。下一個回歸試圖使用5次單項式作為基函數。果不其然，回歸結果（如圖11-3所示）看上去更接近原始函數。但是，它還遠稱不上完美：
In [11]: reg = np.polyfit(x, f(x), deg=5)
ry = np.polyval(reg, x)
In [12]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
['f(x)', 'regression'], ['x', 'f(x)']) 圖11-3 使用最高5次的單項式進行回歸
最後一次嘗試使用7次的單項式作為基函數來計算示例函數的近似值。這次的結果如圖11-4所示，相當有說服力：
圖11-4 7次單項式回歸
In [13]: reg = np.polyfit(x, f(x), 7)
ry = np.polyval(reg, x)
In [14]: np.allclose(f(x), ry)❶
Out[14]: False
In [15]: np.mean((f(x) - ry) ** 2) ❷
Out[15]: 0.0017769134759517689
In [16]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
['f(x)', 'regression'], ['x', 'f(x)']) ❶ 檢查函數和回歸值是否相同（至少接近）。
❷ 根據函數值計算回歸值均方差（MSE）。
2．單獨的基函數
一般來說，當您選擇更好的基函數組時，可以得到更好的回歸結果，例如利用對函數的認識進行近似值計算。在這種情況下，單獨的基函數必須通過一個矩陣方法定義（也就是使用NumPy的ndarray對象）。首先，例子中的多項式最高為3次（圖11-5）。本例的核心函數是np.linalg.lstsq()：
In [17]: matrix = np.zeros((3 + 1, len(x))) ❶
matrix[3, :] = x ** 3 ❷
matrix[2, :] = x ** 2 ❷
matrix[1, :] = x ❷
matrix[0, :] = 1 ❷
In [18]: reg = np.linalg.lstsq(matrix.T, f(x), rcond=None)[0] ❸
In [19]: reg.round(4) ❹
Out[19]: array([ 0. , 0.5628, -0. , -0.0054])
In [20]: ry = np.dot(reg, matrix) ❺
In [21]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
['f(x)', 'regression'], ['x', 'f(x)']) ❶ 基函數值（矩陣）所用的ndarray對象。
❷ 從常數到三次基函數值。
❸ 回歸步驟。
❹ 最優回歸參數。
❺ 函數值的回歸估算。
圖11-5 有單獨基函數的回歸
根據前面單項式的經驗，圖11-5中的結果并不真的如預期那麼好。使用更通用的方法可以讓我們利用對示例函數的認識。我們知道函數中有一個sin部分。是以，在基函數中包含一個正弦函數是有意義的。簡單起見，我們替換最高次的單項式。現在的拟合很完美，如圖11-6所示：
In [22]: matrix[3, :] = np.sin(x) ❶
In [23]: reg = np.linalg.lstsq(matrix.T, f(x), rcond=None)[0]
In [24]: reg.round(4) ❷
Out[24]: array([0. , 0.5, 0. , 1. ])
In [25]: ry = np.dot(reg, matrix)
In [26]: np.allclose(f(x), ry) ❸
Out[26]: True
In [27]: np.mean((f(x) - ry) ** 2) ❸
Out[27]: 3.404735992885531e-31
In [28]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
['f(x)', 'regression'], ['x', 'f(x)']) ❶ 新的基函數利用關于示例函數的知識。
❷ 最優回歸參數恢複原始參數。
❸ 現在，回歸産生了完美的拟合。
圖11-6 使用正弦基函數的回歸
3．有噪聲的資料
回歸對于有噪聲的資料同樣能夠很好的處理，這種資料來自于模拟或者（不完善的）測量。為了闡述這個要點，我們生成同樣具有噪聲的自變量觀測值和因變量觀測值。圖11-7表明，回歸結果比有噪聲的資料點更接近原始函數。在某種意義上，回歸在一定程度上平均了噪聲：
圖11-7 使用有噪聲資料的回歸
In [29]: xn = np.linspace(-2 * np.pi, 2 * np.pi, 50) ❶
xn = xn + 0.15 * np.random.standard_normal(len(xn)) ❷
yn = f(xn) + 0.25 * np.random.standard_normal(len(xn)) ❸
In [30]: reg = np.polyfit(xn, yn, 7)
ry = np.polyval(reg, xn)
In [31]: create_plot([x, x], [f(x), ry], ['b', 'r.'],
['f(x)', 'regression'], ['x', 'f(x)']) ❶ 新的x确定值。
❷ 在x值中引入噪聲。
❸ 在y值中引入噪聲。
4．未排序資料
回歸的另一個重要特點是，它可以無縫地處理未排序資料。前面的例子都依賴于經過排序的x資料，情況并不總是這樣的。為了說明這一點，我們随機生成自變量資料點。在這種情況下，僅靠從視覺上檢查原始資料很難識别出任何結構：
In [32]: xu = np.random.rand(50) * 4 * np.pi - 2 * np.pi ❶
yu = f(xu)
In [33]: print(xu[:10].round(2)) ❶
print(yu[:10].round(2)) ❶
[-4.17 -0.11 -1.91 2.33 3.34 -0.96 5.81 4.92 -4.56 -5.42]
[-1.23 -0.17 -1.9 1.89 1.47 -1.29 2.45 1.48 -1.29 -1.95]
In [34]: reg = np.polyfit(xu, yu, 5)
ry = np.polyval(reg, xu)
In [35]: create_plot([xu, xu], [yu, ry], ['b.', 'ro'],
['f(x)', 'regression'], ['x', 'f(x)']) ❶ 随機化x值。
和有噪聲資料一樣，回歸方法不關心觀測點的順序。這在研究公式11-1所示的最小化問題的結構時很明顯。從圖11-8中顯示的結果來看也很明顯。
圖11-8 使用未排序資料的回歸
5．多元
最小二乘回歸方法的另一個優點是，不需要太多的修改就可以用于多元的情況。接下來以fm()函數為例進行講解：
In [36]: def fm(p):
x, y = p
return np.sin(x) + 0.25 * x + np.sqrt(y) + 0.05 * y ** 2 為了正确地可視化這個函數，我們需要自變量資料點的網格（在兩個次元上）。圖11-9根據以x、y、z表示的自變量和因變量二維資料點網格，顯示了fm()函數的形狀：
In [37]: x = np.linspace(0, 10, 20)
y = np.linspace(0, 10, 20)
X, Y = np.meshgrid(x, y) ❶
In [38]: Z = fm((X, Y))
x = X.flatten() ❷
y = Y.flatten() ❷
In [39]: from mpl_toolkits.mplot3d import Axes3D ❸
In [40]: fig = plt.figure(figsize=(10, 6))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
cmap='coolwarm', linewidth=0.5,
antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
fig.colorbar(surf, shrink=0.5, aspect=5) ❶ 從一維ndarray對象生成二維ndarray對象（網格）。
❷ 從二維ndarray對象得到一維ndarray對象。
❸ 必要時從matplotlib導入3D繪圖功能。
圖11-9 使用兩個參數的函數
為了獲得好的回歸結果，我們将應用基本函數集，包括fm()函數、np.sin()和np.sqrt()函數。圖11-10直覺展示了完美的回歸結果：
In [41]: matrix = np.zeros((len(x), 6 + 1))
matrix[:, 6] = np.sqrt(y) ❶
matrix[:, 5] = np.sin(x) ❷
matrix[:, 4] = y ** 2
matrix[:, 3] = x ** 2
matrix[:, 2] = y
matrix[:, 1] = x
matrix[:, 0] = 1
In [42]: reg = np.linalg.lstsq(matrix, fm((x, y)), rcond=None)[0]
In [43]: RZ = np.dot(matrix, reg).reshape((20, 20)) ❸
In [44]: fig = plt.figure(figsize=(10, 6))
ax = fig.gca(projection='3d')
surf1 = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
cmap=mpl.cm.coolwarm, linewidth=0.5,
antialiased=True) ❹
surf2 = ax.plot_wireframe(X, Y, RZ, rstride=2, cstride=2,
label='regression') ❺
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
ax.legend()
fig.colorbar(surf, shrink=0.5, aspect=5) ❶ 用于y參數的np.sqrt()函數。
❷ 用于x參數的np.sin()函數。
❸ 将回歸結果轉化為網格結構。
❹ 繪制原始函數曲面。
❺ 繪制回歸曲面。
圖11-10 雙參數函數的回歸曲面
回歸
最小二乘回歸方法有多種應用領域，包括簡單的函數逼近和基于有噪聲或者未排序資料的函數逼近。這些方法可以應用于一維問題，也可以應用于多元問題。由于這種方法的基礎數學理論，是以它在一維問題和多元問題上的應用總是“幾乎相同”。
11.1.2 插值
與回歸相比，插值（例如，3次樣條插值）在數學上更為複雜。它還被限制在低次元問題上。給定一組有序的觀測點（按照x維排序），基本的思路是在兩個相鄰資料點之間進行回歸，這樣做不僅産生的分段插值函數完全比對資料點，而且函數在資料點上連續可微分。連續可微分性需要至少3階插值——也就是3次樣條插值。然而，這種方法一般也适用于4次或者線性樣條插值。
下面的代碼可實作線性樣條插值，結果如圖11-11所示：
In [45]: import scipy.interpolate as spi ❶
In [46]: x = np.linspace(-2 * np.pi, 2 * np.pi, 25)
In [47]: def f(x):
return np.sin(x) + 0.5 * x
In [48]: ipo = spi.splrep(x, f(x), k=1) ❷
In [49]: iy = spi.splev(x, ipo) ❸
In [50]: np.allclose(f(x), iy) ❹
Out[50]: True
In [51]: create_plot([x, x], [f(x), iy], ['b', 'ro'],
['f(x)', 'interpolation'], ['x', 'f(x)']) ❶ 從SciPy導入必要的子庫。
❷ 實作線性樣條插值。
❸ 得出内插值。
❹ 檢查内插值是否（足夠）接近函數值。
圖11-11 線性樣條插值（完整資料集）
如果有按照x值排序的一組資料點，那麼應用本身也和使用np.polyfit()和np.polyval()函數一樣簡單。在本例中，它們對應的函數是sci.splrep()和sci.splev()。表11-2列出了sci.splrep()函數的主要參數。
表11-3列出sci.splev()函數的參數。
樣條插值在金融學中往往用于估算未包含在原始觀測點中的自變量資料點的因變量值。為此，在下個例子中選擇一個更小的區間，仔細觀察一次樣條插入的值。圖11-12說明，插值函數确實可以線性地在兩個觀測點之間插值。對于某些應用，這可能不夠精确。此外，很明顯函數在原始資料點上不是連續可微分的——這是另一個不足：
In [52]: xd = np.linspace(1.0, 3.0, 50) ❶
iyd = spi.splev(xd, ipo)
In [53]: create_plot([xd, xd], [f(xd), iyd], ['b', 'ro'],
['f(x)', 'interpolation'], ['x', 'f(x)']) ❶ 具有更多資料點的較小區間。
圖11-12 線性樣條插值（資料子集）
重複整個練習，這次使用3次樣條，結果明顯改善（見圖11-13）：
In [54]: ipo = spi.splrep(x, f(x), k=3) ❶
iyd = spi.splev(xd, ipo) ❷
In [55]: np.allclose(f(xd), iyd) ❸
Out[55]: False
In [56]: np.mean((f(xd) - iyd) ** 2) ❹
Out[56]: 1.1349319851436892e-08
In [57]: create_plot([xd, xd], [f(xd), iyd], ['b', 'ro'],
['f(x)', 'interpolation'], ['x', 'f(x)']) ❶ 完整資料集上的3次樣條插值。
❷ 結果應用到更小的時間間隔。
❸ 插值仍然不完美。
❹ 但好于從前。
插值
在可以應用樣條插值的情況下，可以期望得到比最小二乘回歸方法更好的近似結果。但是要記住，必須有排序（且“無噪聲”）的資料，該方法僅限于低次元問題。樣條插值的計算要求也更高，在某些用例中可能導緻花費的時間比回歸方法的長得多。
圖11-13 三次樣條插值（資料子集）
本文摘自《Python金融大資料分析》（第2版）