- 擬合是用一個連續函數（曲線）靠近給定的離散數據，使其與給定的數據相吻合。
- 數據擬合的演算法相對比較簡單，但調用不同工具和方法時的函數定義和參數設置有所差異，往往使小白感到困惑。
- 本文基於 Scipy 工具包，對單變數、多變數線性最小二乘擬合，指數函數、多項式函數、樣條函數的非線性擬合，單變數、多變數的自定義函數擬合問題進行分析、給出完整常式和結果，數據擬合從此無憂。
1. 數據擬合
在科學研究和工程應用中經常通過測量、採樣、實驗等方法獲得各種數據。對一組已知數據點集，通過調整擬合函數（曲線）的參數，使該函數與已知數據點集相吻合，這個過程稱為數據擬合，又稱曲線擬合。
插值和擬合都是根據一組已知數據點，求變化規律和特徵相似的近似曲線的過程。但是插值要求近似曲線完全經過所有的給定數據點，而擬合只要求近似曲線在整體上儘可能接近數據點，並反映數據的變化規律和發展趨勢。因此插值可以看作是一種特殊的擬合，是要求誤差函數為 0 的擬合。
1.1 數據擬合問題的分類
數據擬合問題，可以從不同角度進行分類：
- 按照擬合函數分類，分為線性函數和非線性函數。非線性函數用於數據擬合，常用的有多項式函數、樣條函數、指數函數和冪函數，針對具體問題還有自定義的特殊函數顯示。
- 按照變數個數分類，分為單變數函數和多變數函數。
- 按照擬合模型分類，分為基於模型的數據擬合和無模型的函數擬合。基於模型的數據擬合，是通過建立數學模型描述輸入輸出變數之間的關係，擬合曲線不僅能擬合觀測數據，擬合模型的參數通常具有明確的物理意義。而無模型的函數擬合，是指難以建立描述變數關係的數學模型，只能採用通用的函數和曲線擬合觀測數據，例如多項式函數擬合、樣條函數擬合，也包括機器學習和神經網路模型，這種模型的參數通常沒有明確的意義。
1.2 數據擬合的原理和方法
數據擬合通過調整擬合函數中的待定參數，從整體上接近已知的數據點集。
這是一個優化問題，決策變數是擬合函數的待定參數，優化目標是觀測數據與擬合函數的函數值之間的某種誤差指標。典型的優化目標是擬合函數值與觀測值的誤差平方和；當觀測數據的重要性不同或分布不均勻時，也可以使用加權誤差平方和作為優化目標。
數據擬合的基本方法是最小二乘法。對於觀測數據 ( x i , y i ) , i = 1 , . . n，將觀測值 y i與擬合函數 y = f ( x , p ) 的計算值 f ( x i ) 的誤差平方和最小作為優化問題的目標函數：
( p 1 , ⋯ p m ) 是擬合函數中的待定參數。
對於線性擬合問題，設擬合函數為直線 f ( x ) = p 0 + p 1 ∗ x ， 由極值的必要條件 $ partial J/partial p_j = 0,; (j=0,1)$ 可以解出係數 p 0 , p 1 ：
對於多變數線性最小二乘問題，設擬合函數為直線 f ( x ) = p 0 + p 1 ∗ x 1 + ⋯ + p m ∗ x m ， 類似地，可以解出係數 p 0 , p 1 , ⋯ p m 。
對於非線性函數的擬合問題，通常也是按照最小二乘法的思路，求解上述誤差平方和最小化這個非線性優化問題，常用的具體演算法有搜索演算法和迭代演算法兩類。
1.3 Python 數據擬合方法
數據擬合是常用演算法，Python 語言的很多工具包都提供了數據擬合方法，常用的如 Scipy、Numpy、Statsmodel、Scikit-learn 工具包都帶有數據擬合的函數與應用。
Scipy 是最常用的 Python 工具包，本系列中非線性規劃、插值方法也都是使用 Scipy 工具包實現，因此仍以 Scipy 工具包講解數據擬合問題。
Scipy 工具包對於不同類型的數據擬合問題，提供了不同的函數或類。由於 Scipy 工具包是多個團隊合作完成，而且經過了不斷更新，因此調用不同函數和方法時的函數定義和參數設置有所差異，往往使小白感到困惑。
本文對單變數、多變數線性最小二乘擬合，指數函數、多項式函數、樣條函數的非線性擬合，單變數、多變數的自定義函數擬合問題進行分析、給出完整常式和結果，數據擬合從此無憂。
2. 線性最小二乘擬合
2.1 線性最小二乘擬合函數說明
線性最小二乘擬合是最簡單和最常用的擬合方法。scipy.optimize 工具箱中的 leastsq()、lsq_linear()，scipy.stats 工具箱中的 linregress()，都可以實現線性最小二乘擬合。
2.1.1 scipy.optimize.leastsq 函數說明
leastsq() 根據觀測數據進行最小二乘擬合計算，只需要觀測值與擬合函數值的誤差函數和待定參數 的初值，返回擬合函數中的待定參數 ( p 1 , ⋯ p m ) ，但不能提供參數估計的統計信息。leastsq() 可以進行單變數或多變數線性最小二乘擬合，對變數進行預處理後也可以進行多項式函數擬合。
scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)
主要參數：
- func：可調用的函數，描述擬合函數的函數值與觀測值的誤差，形式為 error(p,x,y)，具有一個或多個待定參數 p。誤差函數的參數必須按照 (p,x,y) 的順序排列，不能改變。
- x0：一維數組，待定參數 ( p 1 , ⋯ p m ) 的初值。
- args：元組，線性擬合時提供觀測數據值 (xdata, ydata)，觀測數據 xdata 可以是一維數組（單變數問題），也可以是多維數組（多變數問題）。
返回值：
- x：一維數組，待定參數 ( p 1 , ⋯ p m )的最小二乘估計值。
2.1.2 scipy.stats.linregress 函數說明
linregress() 根據兩組觀測數據 (x,y) 進行線性最小二乘回歸，不僅返回擬合函數中的待定參數 ( p 1 , p 1 )，而且可以提供參數估計的各種統計信息，但只能進行單變數線性擬合。
scipy.stats.linregress(x, y=None, alternative=『two-sided』)
主要參數：
- x, y：x, y 是長度相同的一維數組。或者 x 是二維數組，且 y=none，則二維數組 x 相當於 長度相同的一維數組 x, y。
返回值：
- slope：斜率，直線 f ( x ) = p 0 + p 1 ∗ x中的 p 1。
- intercept：截距，直線 f ( x ) = p 0 + p 1 ∗ x中的 p 0 。
- rvalue：r^2 值，統計量。
- pvalue：p 值，P檢驗的統計量。
- stderr：標準差，統計量。
2.2 Python 常式：單變數線性擬合
程序說明：
- scipy.optimize.leastsq() 與 scipy.stats.linregress() 都可以進行單變數線性擬合。leastsq() 既可以用於單變數也可以用於多變數問題；linregress() 只能用於單變數問題，但可以給出很多參數估計的統計結果。
- leastsq() 要以子函數來定義觀測值與擬合函數值的誤差函數，常式中分別定義了擬合函數 fitfunc1(p, x) 與誤差函數error1(p, x, y) ，是為了方便調用擬合函數計算擬合曲線在數據點的函數值。注意 p 為數組 。
- leastsq() 中誤差函數的函數名可以任意定義，但誤差函數的參數必須按照 (p,x,y) 的順序排列，不能改變次序。
- leastsq() 中觀測數據 (x, yObs) 是以動態參數 args 的方式進行傳遞的。這種處理方式非常獨特，沒有為什麼， leastsq() 就是這樣定義的。
- linregress() 只要將觀測數據 (x,yObs) 作為參數，默認單變數線性擬合，不需要定義子函數。
- leastsq() 與 linregress() 進行線性擬合，得到的參數估計結果是相同的。
Python 常式：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 1. 單變數線性擬合：最小二乘法 scipy.optimize.leastsq
import numpy as np
import matplotlib.pyplot as plt # 導入 Matplotlib 工具包
from scipy.optimize import leastsq # 導入 scipy 中的最小二乘法擬合工具
from scipy.stats import linregress # 導入 scipy 中的線性回歸工具
def fitfunc1(p, x): # 定義擬合函數為直線
p0, p1 = p # 擬合函數的參數
y = p0 + p1*x # 擬合函數的表達式
return y
def error1(p, x, y): # 定義觀測值與擬合函數值的誤差函數
err = fitfunc1(p,x) - y # 誤差
return err
# 創建給定數據點集 (x,yObs)
p = [2.5, 1.5] # y = p[0] + p[1] * x
x = np.array([0., 0.5, 1.5, 2.5, 4.5, 5.5, 7.5, 8.0, 8.5, 9.0, 10.0])
y = p[0] + p[1] * x # 理論值 y
np.random.seed(1)
yObs = y + np.random.randn(x.shape[-1]) # 生成帶有雜訊的觀測數據
# print(x.shape, y.shape, yObs.shape)
# 由給定數據點集 (x,y) 求擬合函數的參數 pFit
p0 = [1, 1] # 設置擬合函數的參數初值
pFit, info = leastsq(error1, p0, args=(x,yObs)) # 最小二乘法求擬合參數
print("Data fitting with Scipy.optimize.leastsq")
print("y = p[0] + p[1] * x")
print("p[0] = {:.4f}np[1] = {:.4f}".format(pFit[0], pFit[1]))
# 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
yFit = fitfunc1(pFit,x)
# 比較：線性回歸，可以返回斜率，截距，r 值，p 值，標準誤差
slope, intercept, r_value, p_value, std = linregress(x, yObs)
print("nLinear regress with Scipy.stats.linregress")
print("y = p[0] + p[1] * x")
print("p[0] = {:.4f}".format(intercept)) # 輸出截距 intercept
print("p[1] = {:.4f}".format(slope)) # 輸出斜率 slope
print("r^2_value: {:.4f}".format(r_value**2)) # 輸出 r^2 值
print("p_value: {:.4f}".format(p_value)) # 輸出 p 值
print("std: {:.4f}".format(std)) # 輸出標準差 std
# 繪圖
fig, ax = plt.subplots(figsize=(8,6))
ax.text(8,3,"youcans-xupt",color='gainsboro')
ax.set_title("Data fitting with linear least squares")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'r--', label="theoretical curve")
plt.plot(x, yFit, 'b-', label="fitting curve")
plt.legend(loc="best")
plt.show()程序運行結果：
Data fitting with Scipy.optimize.leastsq
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528
Linear regress with Scipy.stats.linregress
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528
r^2_value: 0.9521
p_value: 0.0000
std: 0.1161
2.3 Python 常式：多變數線性擬合
程序說明：
- scipy.optimize.leastsq() 既可以用於單變數也可以用於多變數問題，本常式求解一個二元線性擬合問題：y = p[0] + p[1] * x1 + p[2] * x2。
- leastsq() 求解多變數問題的方法與單變數問題類似，以子函數 error2(p, x1, x2, y) 來定義觀測值與擬合函數值的誤差函數，以動態參數 args 的方式傳遞觀測數據 (x, yObs) 。
- leastsq() 中誤差函數的函數名可以任意定義，但誤差函數的參數必須按照 (p,x1,x2,y) 的順序排列，不能改變次序。
- scipy 只能做一元線性回歸，常式中通過調用 statsmodels.api 進行多元線性回歸，可以得到各種統計參數，供讀者參考。
Python 常式：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 2. 多變數線性擬合：最小二乘法 scipy.optimize.leastsq
import numpy as np
import matplotlib.pyplot as plt # 導入 Matplotlib 工具包
from scipy.optimize import leastsq # 導入 scipy 中的最小二乘法工具
def fitfunc2(p, x1, x2): # 定義擬合函數為直線
p0, p1, p2 = p # 擬合函數的參數
y = p0 + p1*x1 + p2*x2 # 擬合函數的表達式
return y
def error2(p, x1, x2, y): # 定義觀測值與擬合函數值的誤差函數
err = fitfunc2(p, x1, x2) - y # 計算殘差
return err
# 創建給定數據點集 (x,yObs)
np.random.seed(1)
p = [2.5, 1.5, -0.5] # y = p[0] + p[1] * x1 + p[2] * x2
x1 = np.array([0., 0.5, 1.5, 2.5, 4.5, 5.5, 7.5, 8.0, 8.5, 9.0, 10.0])
x2 = np.array([0., 1.0, 1.5, 2.2, 4.8, 5.0, 6.3, 6.8, 7.1, 7.5, 8.0])
z = p[0] + p[1]*x1 + p[2]*x2 # 理論值 z
zObs = z + np.random.randn(x1.shape[-1]) # 生成帶有雜訊的觀測數據
print(x1.shape, z.shape, zObs.shape)
# 由給定數據點集 (x,z) 求擬合函數的參數 pFit
p0 = [1, 1, 1] # 設置擬合函數的參數初值
pFit, info = leastsq(error2, p0, args=(x1,x2,zObs)) # 最小二乘法求擬合參數
print("Data fitting with Scipy.optimize.leastsq:")
print("z = p[0] + p[1]*x1 + p[1]*x2")
print("p[0]={:.4f}np[1]={:.4f}np[2]={:.4f}".format(pFit[0], pFit[1], pFit[2]))
# 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
zFit = fitfunc2(pFit, x1, x2)
# 多元線性回歸：最小二乘法(OLS)
import statsmodels.api as sm
x0 = np.ones(x1.shape[-1]) # 截距列 x0=[1,...1]
X = np.column_stack((x0, x1, x2)) # (nSample,3): [x0, x1, x2]
model = sm.OLS(zObs, X) # 建立 OLS 模型: y = b0*x0 + b1*x1 + b2*x2 + e
results = model.fit() # 返回模型擬合結果
zFit = results.fittedvalues # 模型擬合的 y值
print(results.summary()) # 輸出回歸分析的摘要
print("nOLS model: y = b0*x0 + b1*x1 + b2*x2")
print('Parameters: ', results.params) # 輸出：擬合模型的係數程序運行結果：
Data fitting with Scipy.optimize.leastsq
z = p[0] + p[1]*x1 + p[1]*x2
p[0]=2.6463
p[1]=2.2410
p[2]=-1.3710
OLS Regression Results
OLS model: y = b0*x0 + b1*x1 + b2*x2
Parameters: [ 2.64628055 2.24100973 -1.37104475]
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 2.6463 0.942 2.808 0.023 0.473 4.820
x1 2.2410 1.043 2.148 0.064 -0.165 4.647
x2 -1.3710 1.312 -1.045 0.326 -4.396 1.654
==============================================================================3. 非線性函數數據擬合
3.1 非線性擬合函數說明
非線性函數是非常廣泛的概念。本節討論指數函數、多項式函數和樣條函數三種常用的通用形式的非線性函數擬合問題，分別使用了 Scipy 工具包中的 scipy.optimize.leastsq()、scipy.linalg.lstsq() 和
scipy.interpolate.UnivariateSpline() 函數。
scipy.optimize.leastsq() 的使用方法已在本文 2.1 中進行了介紹，
scipy.interpolate.UnivariateSpline() 的使用方法在《22. 插值方法》文中進行了介紹，以下介紹 scipy.linalg.lstsq() 函數。
lstsq() 函數只要傳入觀測數據 (x,yObs)，並將 x 按多項式階數轉換為 X，即可求出多項式函數的係數，不需要定義擬合函數或誤差函數，非常適合比較不同階數的多項式函數擬合的效果。
scipy.linalg.lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None)
主要參數：
- a：(m,n) 數組，表示方程 Ax=b 的左側。
- b：(m,) 數組，表示方程 Ax=b 的右側。
返回值：
- x：最小二乘的解，指多項式函數的係數。
注意：lstsq() 函數中求解方程 Ax=b，A 是指由觀測數據 x i x_ix 按多項式階數轉換為矩陣 ( x i 0 , x i 1 , . . . x i m ) , i = 1 , ,n，b 是指 y i , i = 1 ,,n，而 x 是指多項式函數的係數，詳見常式。
3.2 Python 常式：指數函數擬合
程序說明：
- scipy.optimize.leastsq() 本質上是求解帶有待定參數的誤差函數最小化問題，因此可以用於指數函數的最小二乘擬合。類似地，原理上 leastsq() 也可以用於其它形式非線性函數的擬合問題，但對於某些函數擬合誤差可能會比較大。
- leastsq() 以子函數 error3(p, x, y) 來定義觀測值與擬合函數值的誤差函數，以動態參數 args 的方式傳遞觀測數據 (x, yObs) 。
- leastsq() 中誤差函數的函數名可以任意定義，但誤差函數的參數必須按照 (p,x1,x2,y) 的順序排列，不能改變次序。
Python 常式：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 3. 非線性函數擬合：指數函數擬合(exponential function)
import numpy as np
import matplotlib.pyplot as plt # 導入 Matplotlib 工具包
from scipy.optimize import leastsq # 導入 scipy 中的最小二乘法工具
def fitfunc3(p, x): # 定義擬合函數
p0, p1, p2 = p # 擬合函數的參數
y = p0 + p1 * np.exp(-p2*x) # 擬合函數的表達式
return y
def error3(p, x, y): # 定義殘差函數
err = fitfunc3(p,x) - y # 殘差
return err
# 創建給定數據點集 (x,yObs)
p = [0.5, 2.5, 1.5] # y = p0 + p1 * np.exp(-p2*x)
x = np.linspace(0, 5, 50)
y = fitfunc3(p, x)
np.random.seed(1)
yObs = y + 0.2*np.random.randn(x.shape[-1]) # 生成帶有雜訊的觀測數據
# print(x.shape, y.shape, yObs.shape)
# 由給定數據點集 (x,y) 求擬合函數的參數 pFit
p0 = [1, 1, 1] # 設置擬合函數的參數初值
pFit, info = leastsq(error3, p0, args=(x,yObs)) # 最小二乘法求擬合參數
print("Data fitting of exponential function")
print("y = p0 + p1 * np.exp(-p2*x)")
print("p[0] = {:.4f}np[1] = {:.4f}np[2] = {:.4f}".format(pFit[0], pFit[1], pFit[2]))
# 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
yFit = fitfunc3(pFit, x)
# 繪圖
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Data fitting of exponential function")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'r--', label="theoretical curve")
plt.plot(x, yFit, 'b-', label="fitting curve")
plt.legend(loc="best")
plt.show()程序運行結果：
Data fitting of exponential function
y = p0 + p1 * np.exp(-p2*x)
p[0] = 0.5216
p[1] = 2.5742
p[2] = 1.6875
3.3 Python 常式：多項式函數擬合
程序說明：
- scipy.optimize.leastsq() 本質上是求解帶有待定參數的誤差函數最小化問題，因此可以用於多項式函數的最小二乘擬合，使用方法與線性擬合、指數擬合類似。
- 由於 leastsq() 要以子函數 error(p, x, y) 來定義觀測值與擬合函數值的誤差函數，在比較不同階數的多項式函數擬合時需要定義多個對應的誤差函數，比較繁瑣。
- scipy.linalg.lstsq() 只要傳入觀測數據 (x,yObs)，並將 x 按多項式階數轉換為 X，即可求出多項式函數的係數，不需要定義擬合函數或誤差函數，非常適合比較不同階數的多項式函數擬合的效果。
- 對於相同階數的多項式函數，leastsq() 與 lstsq() 的參數估計和擬合結果是相同的。
- 增大多項式的階數，可以減小擬合曲線與觀測數據的誤差平方和，但也更容易導致過擬合，雖然能更好地擬合訓練數據，但並不能真實反映數據的總體規律，因而對於訓練數據以外的測試數據的擬合效果反而降低了。
Python 常式：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 4. 非線性函數擬合：多項式函數擬合(Polynomial function)
import numpy as np
import matplotlib.pyplot as plt # 導入 Matplotlib 工具包
from scipy.linalg import lstsq # 導入 scipy 中的 linalg, stats 函數庫
from scipy.optimize import leastsq # 導入 scipy 中的最小二乘法工具
def fitfunc4(p, x): # 定義擬合函數
p0, p1, p2, p3 = p # 擬合函數的參數
y = p0 + p1*x + p2*x*x + p3*x*x*x # 擬合函數的表達式
return y
def error4(p, x, y): # 定義觀測值與擬合函數值的誤差函數
err = fitfunc4(p,x) - y # 殘差
return err
# 創建給定數據點集 (x,yObs)
p = [1.0, 1.2, 0.5, 0.8] # y = p0 + ((x*x-p1)**2+p2) * np.sin(x*p3)
func = lambda x: p[0]+((x*x-p[1])**2+p[2])*np.sin(x*p[3]) # 定義 y=f(x)
x = np.linspace(-1, 2, 30)
y = func(x) # 計算已知數據點的理論值 y =f(x)
np.random.seed(1)
yObs = y + 0.5*np.random.randn(x.shape[-1]) # 生成帶有雜訊的觀測數據 yObs
# 繪圖
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Polynomial fitting with least squares")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'c--', label="theoretical curve")
# 用 scipy.optimize.leastsq() 進行多項式函數擬合
p0 = [1, 1, 1, 1] # 設置擬合函數的參數初值
pFit, info = leastsq(error4, p0, args=(x,yObs)) # 最小二乘法求擬合參數
yFit = fitfunc4(pFit, x) # 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
ax.plot(x, yFit, '-', label='leastsq')
print("Polynomial fitting by scipy.optimize.leastsq")
print("y = p[0] + p[1]*x + p[2]*x^2 +p[3]*x^3") # 擬合函數的表達式
print("p[0] = {:.4f}np[1] = {:.4f}np[2] = {:.4f}np[3] = {:.4f}".format(pFit[0], pFit[1], pFit[2], pFit[3]))
# 用 scipy.linalg.lstsq() 進行多項式函數擬合
print("nPolynomial fitting by scipy.linalg.lstsq")
print("y = w[0] + w[1]*x + w[2]*x^2 + ... + w[m]*x^m") # 擬合函數的表達式
# 最小二乘法多項式數據擬合，求解多項式函數的係數 W=[w[0],...w[m]]
for order in range(1,5):
# order = 3 # 多項式階數, m
X = np.array([[(xi ** i) for i in range(order + 1)] for xi in x])
Y = np.array(yObs).reshape((-1, 1))
W, res, rnk, s = lstsq(X, Y) # 最小二乘法求解 A*W = b
print("order={:d}".format(order))
for i in range(order+1): # W = [w[0],...w[order]]
print("tw[{:d}] = {:.4f}".format(i, W[i,0]))
# 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
yFit = X.dot(W) # y = w[0] + w[1]*x + w[2]*x^2 + ... + w[m]*x^m
# 繪圖：n 次多項式函數擬合曲線
ax.plot(x, yFit, '-', label='order={}'.format(order))
plt.legend(loc="best")
plt.show()程序運行結果：
Polynomial fitting by scipy.optimize.leastsq
y = p[0] + p[1]*x + p[2]*x^2 +p[3]*x^3
p[0] = 0.9586
p[1] = -0.4745
p[2] = -0.3581
p[3] = 1.1721
Polynomial fitting by scipy.linalg.lstsq
y = w[0] + w[1]*x + w[2]*x^2 + ... + w[m]*x^m
order=1
w[0] = 1.0331
w[1] = 1.7354
order=2
w[0] = 0.2607
w[1] = 0.3354
w[2] = 1.4000
order=3
w[0] = 0.9586
w[1] = -0.4745
w[2] = -0.3581
w[3] = 1.1721
order=4
w[0] = 1.0131
w[1] = 1.6178
w[2] = -1.1039
w[3] = -1.5209
w[4] = 1.3465
3.4 Python 常式：樣條曲線擬合
程序說明：
scipy.interpolate.UnivariateSpline() 類是一種基於固定數據點創建函數的方法，使用樣條曲線擬合到給定的數據點集。
UnivariateSpline 類由已知數據點集生成樣條插值函數 y=spl(x)，通過調用樣條插值函數可以計算指定 x 的函數值 f(x)。
UnivariateSpline 類既可以進行數據插值，也可以進行擬合。參數 s=0 表示數據插值，樣條曲線必須通過所有數據點；s>0 表示數據擬合，默認 s= len(w)。
通過 set_smoothing_factor(sf) 設置光滑因子，可以對樣條擬合函數進行調節，使擬合曲線更好地反映觀測數據特徵，避免過擬合。
Python 常式：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 5 非線性函數擬合：樣條函數擬合(Spline function)
import numpy as np
import matplotlib.pyplot as plt # 導入 Matplotlib 工具包
from scipy.interpolate import UnivariateSpline # 導入 scipy 中的樣條插值工具
# 創建給定數據點集 (x,yObs)
# y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
np.random.seed(1)
p0, p1, p2, p3 = [1.0, 1.2, 0.5, 0.8] # 擬合函數的參數
x = np.linspace(-1, 2, 30) # 生成已知數據點集的 x
y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3) # 生成理論值 y
yObs = y + 0.5*np.random.randn(x.shape[-1]) # 生成帶有雜訊的觀測數據
# 由給定數據點集 (x,y) 求擬合函數的參數 fSpl
fSpl = UnivariateSpline(x, yObs) # 三次樣條插值，默認 s= len(w)
coeffs = fSpl.get_coeffs() # Return spline coefficients
print("Data fitting with spline function")
print("coeffs of 3rd spline function:n ", coeffs)
# 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
yFit = fSpl(x) # 由插值函數 fSpl1 計算插值點的函數值 yFit
# 對擬合函數 fitfunc 進行平滑處理
fSpl.set_smoothing_factor(10) # 設置光滑因子 sf
yS = fSpl(x) # 由插值函數 fSpl(sf=10) 計算插值點的函數值 ys
coeffs = fSpl.get_coeffs() # 平滑處理後的參數
print("coeffs of 3rd spline function (sf=10):n ", coeffs)
# 繪圖
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Data fitting with spline function")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'r--', label="theoretical curve")
plt.plot(x, yFit, 'b-', label="3rd spline fitting")
plt.plot(x, yS, 'm-', label="smoothing spline")
plt.legend(loc="best")
plt.show()程序運行結果：
Data fitting with spline function
coeffs of 3rd spline function:
[-0.09707885 3.66083026 -4.20416235 7.95385344]
coeffs of 3rd spline function (sf=10):
[0.41218039 0.52795588 1.6248287 0.76540737 8.49462738]
4. 自定義函數曲線擬合
4.1 scipy.optimize.curve_fit() 函數說明
curve_fit() 使用非線性最小二乘法將自定義的擬合函數擬合到觀測數據，不僅可以用於直線、二次曲線、三次曲線的擬合，而且可以適用於任意形式的自定義函數的擬合，使用非常方便。curve_fit() 允許進行單變數或多變數的自定義函數擬合。
**scipy.optimize.curve_fit(f,xdata,ydata,p0=None,sigma=None,absolute_sigma=False,check_finite=True,bounds=(-inf,inf),method=None,jac=None,kwargs) **
主要參數：
- f：可調用的函數，自定義的擬合函數，具有一個或多個待定參數。擬合函數的形式為 func(x,p1,p2,…)，其中參數必須按照 (x,p1,p2,…) 的順序排列，p1, p2,… 是標量不能表達為數組。
- xdata：n*m數組，n 為觀測數據長度，m為變數個數。觀測數據 xdata 可以是一維數組（單變數問題），也可以是多維數組（多變數問題）。
- ydata：數組，長度為觀測數據長度 n。
- p0：可選項，待定參數 [p1,p2,…] 的初值，默認值無。
返回值：
- popt：待定參數 ( p 1 , ⋯ p m ) 的最小二乘估計值。
- pcov：參數 ( p 1 , ⋯ p m )的估計值 popt 的協方差，其對角線是各參數的方差。
4.2 Python 常式：單變數自定義函數曲線擬合
程序說明：
- 不同於 leastsq() 定義觀測值與擬合函數值的誤差函數，scipy.optimize.curve_fit() 直接定義一個自定義的擬合函數，更為直觀和便於理解。
- curve_fit() 定義一個擬合函數，函數名可以任意定義，但擬合函數的參數必須按照 (x,p1,p2,…) 的順序排列，不能改變次序。p1, p2,… 是標量，不能寫成數組。注意 leastsq() 中誤差函數的參數必須按照 (p,x,y) 的順序排列，與 curve_fit() 不同。
- leastsq() 也可以對自定義的擬合函數進行最小二乘擬合。
- 由於本常式中自定義擬合函數使用了觀測數據的實際模型，而不是通用的多項式函數或樣條函數，因此擬合結果不僅能很好的擬合觀測數據，而且能更準確地反映實際模型的趨勢。
Python 常式：
# 6. 自定義函數曲線擬合：單變數
import numpy as np
import matplotlib.pyplot as plt # 導入 Matplotlib 工具包
from scipy.optimize import leastsq, curve_fit # 導入 scipy 中的曲線擬合工具
def fitfunc6(x, p0, p1, p2, p3): # 定義擬合函數為自定義函數
# p0, p1, p2, p3 = p # 擬合函數的參數
y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
return y
# def error6(p, x, y): # 定義觀測值與擬合函數值的誤差函數
# p0, p1, p2, p3 = p
# err = fitfunc6(x, p0, p1, p2, p3) - y # 計算殘差
# return err
# 創建給定數據點集 (x, yObs)
p0, p1, p2, p3 = [1.0, 1.2, 0.5, 0.8] # y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
x = np.linspace(-1, 2, 30)
y = fitfunc6(x, p0, p1, p2, p3)
np.random.seed(1)
yObs = y + 0.5*np.random.randn(x.shape[-1]) # 生成帶有雜訊的觀測數據
# # 用 scipy.optimize.leastsq() 進行函數擬合
# pIni = [1, 1, 1, 1] # 設置擬合函數的參數初值
# pFit, info = leastsq(error6, pIni, args=(x, yObs)) # 最小二乘法求擬合參數
# print("Data fitting of custom function by leastsq")
# print("y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)")
# print("p[0] = {:.4f}np[1] = {:.4f}np[2] = {:.4f}np[3] = {:.4f}"
# .format(pFit[0], pFit[1], pFit[2], pFit[3]))
# 用 scipy.optimize.curve_fit() 進行自定義函數擬合(單變數)
pFit, pcov = curve_fit(fitfunc6, x, yObs)
print("Data fitting of custom function by curve_fit:")
print("y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)")
print("p[0] = {:.4f}np[1] = {:.4f}np[2] = {:.4f}np[3] = {:.4f}"
.format(pFit[0], pFit[1], pFit[2], pFit[3]))
print("estimated covariancepcov:n",pcov)
# 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
yFit = fitfunc6(x, pFit[0], pFit[1], pFit[2], pFit[3])
# 繪圖
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Data fitting of custom function")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'r--', label="theoretical curve")
plt.plot(x, yFit, 'b-', label="fitting curve")
plt.legend(loc="best")
plt.show()程序運行結果：
Data fitting of custom function by curve_fit：
y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
p[0] = 0.9460
p[1] = 1.1465
p[2] = 0.8291
p[3] = 0.6008
estimated covariancepcov:
[[ 0.01341654 0.00523061 -0.01645431 0.00455901]
[ 0.00523061 0.02648836 -0.04442234 0.02821206]
[-0.01645431 -0.04442234 0.20326672 -0.07482843]
[ 0.00455901 0.02821206 -0.07482843 0.0388316 ]]結果分析：
4.3 Python 常式：多變數自定義函數曲線擬合
程序說明：
scipy.optimize.curve_fit() 既可以用於單變數也可以用於多變數問題，本常式求解一個二元非線性擬合問題。
curve_fit() 定義一個擬合函數 fitfunc7(X, p0, p1, p2, p3)，函數名可以任意定義，但擬合函數的參數必須按照 (x,p1,p2,…) 的順序排列，不能改變次序。p1, p2,… 是標量，不能寫成數組。
curve_fit(fitfunc7, X, yObs) 中的 X 是 (n,m) 數組，n 是觀測數據點集的長度，m 是變數個數。
Python 常式：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 7. 自定義函數曲線擬合：多變數
import numpy as np
import matplotlib.pyplot as plt # 導入 Matplotlib 工具包
from scipy.optimize import curve_fit # 導入 scipy 中的曲線擬合工具
def fitfunc7(X, p0, p1, p2, p3): # 定義多變數擬合函數, X 是向量
# p0, p1, p2, p3 = p # 擬合函數的參數
y = p0 + p1*X[0,:] + p2*X[1,:] + p3*np.sin(X[0,:]+X[1,:]+X[0,:]**2+X[1,:]**2)
return y
# 創建給定數據點集 (x,yObs)
p = [1.0, 0.5, -0.5, 5.0] # 自定義函數的參數
p0, p1, p2, p3 = p # y = p0 + p1*x1 + p2*x2 + p3*np.sin(x1+x2+x1^2+x2^2)
np.random.seed(1)
x1 = 2.0 * np.random.rand(8) # 生成隨機數組，長度為 8
x2 = 3.0 * np.random.rand(5) # 生成隨機數組，取值範圍 (0,3.0)
xmesh1, xmesh2 = np.meshgrid(x1, x2) # 生成網格點的坐標 xx,yy (二維數組)
xx1= xmesh1.reshape(xmesh1.shape[0]*xmesh1.shape[1], ) # 將網格點展平為一維數組
xx2= xmesh2.reshape(xmesh2.shape[0]*xmesh2.shape[1], ) # 將網格點展平為一維數組
X = np.vstack((xx1,xx2)) # 生成多變數數組，行數為變數個數
y = fitfunc7(X, p0, p1, p2, p3) # 理論計算值 y=f(X,p)
yObs = y + 0.2*np.random.randn(y.shape[-1]) # 生成帶有雜訊的觀測數據
print(x1.shape,x2.shape,xmesh1.shape,xx1.shape,X.shape)
# 用 scipy.optimize.curve_fit() 進行自定義函數擬合(多變數)
pFit, pcov = curve_fit(fitfunc7, X, yObs) # 非線性最小二乘法曲線擬合
print("Data fitting of multivariable custom function")
print("y = p0 + p1*x1 + p2*x2 + p3*np.sin(x1+x2+x1^2+x2^2)")
for i in range(4):
print("p[{:d}] = {:.4f}tp[{:d}]_fit = {:.4f}".format(i, p[i], i, pFit[i]))
# 由擬合函數 fitfunc 計算擬合曲線在數據點的函數值
yFit = fitfunc7(X, pFit[0], pFit[1], pFit[2], pFit[3])程序運行結果：
Data fitting of multivariable custom function:
y = p0 + p1*x1 + p2*x2 + p3*np.sin(x1+x2+x1^2+x2^2)
p[0] = 1.0000 p[0]_fit = 1.1316
p[1] = 0.5000 p[1]_fit = 0.5020
p[2] = -0.5000 p[2]_fit = -0.5906
p[3] = 5.0000 p[3]_fit = 5.0061
estimated covariancepcov:
[[ 9.51937904e-03 -2.82863223e-03 -5.26393413e-03 -8.51457970e-04]
[-2.82863223e-03 4.88275894e-03 9.39281331e-05 3.73832161e-04]
[-5.26393413e-03 9.39281331e-05 3.86701646e-03 4.65766686e-04]
[-8.51457970e-04 3.73832161e-04 4.65766686e-04 1.85374067e-03]]
原創文章，作者：投稿專員，如若轉載，請註明出處：https://www.506064.com/zh-tw/n/233091.html
微信掃一掃 
支付寶掃一掃 