非線性回歸模型的原理及評估——解決行星軌道的拟合問題

非線性回歸模型——行星軌道

• 概述
• 回歸統計
• 普通和權重最小二乘法
• L o g i s t i c Logistic Logistic回歸
• • 對數幾率分布公式
• 拟合一個行星軌道
• • 問題一
  • 思路及解答
  • • 最小二乘法
  • 問題二
  • • 解答
  • 完整代碼

概述

  在統計學中， 非線性回歸是回歸分析的一種形式，其中觀測資料由函數模組化，該函數是模型參數的非線性組合并且取決于一個或多個獨立變量。 通過逐次逼近的方法拟合資料。

  在非線性回歸中，形式的統計模型 ，

f ( x , β ) = β 1 x β 2 + x . f(x,\beta) = \frac{\beta_1x}{\beta_2 + x}. f(x,β)=β2​+xβ1​x​.

  此函數是非線性的，因為它不能表示為兩個 β \beta β的線性組合。

  系統誤差可能存在于自變量中，但其處理不在回歸分析的範圍内。 如果自變量不是無差錯的，那麼這是一個變量誤差模型 ，也在此範圍之外。

  非線性函數的其他示例包括指數函數 ， 對數函數 ， 三角函數 ， 幂函數 ， 高斯函數和洛倫茲曲線 。

回歸統計

  這個過程的基本假設是模型可以用線性函數近似，即一階泰勒級數 ：

f ( x i , β ) ≈ f ( x i , 0 ) + ∑ j J i j β j f(x_i,\beta) \approx f(x_i,0) + \sum_{j} J_{ij}\beta_j f(xi​,β)≈f(xi​,0)+j∑​Jij​βj​

  其中 J I J = ∂ f ( x i , β ) ∂ β j J_{IJ} = \frac{\partial f(x_i,\beta)}{\partial \beta_j} JIJ​=∂βj​∂f(xi​,β)​，由此得出最小二乘估計量由下式給出。

β ^ ≈ ( J T J ) − 1 J T y \hat \beta \approx (J_TJ)^{-1}J^Ty β^​≈(JT​J)−1JTy

  計算非線性回歸統計量并将其用作線性回歸統計量，但在公式中使用J代替X. 線性近似将偏差引入統計中。 是以，在解釋從非線性模型得到的統計資料時，需要比平常更多的謹慎。

普通和權重最小二乘法

  最佳拟合曲線通常假定應該看起來平方的總和最小化殘差 。 這是普通的最小二乘 （OLS）方法。 然而，在因變量不具有恒定方差的情況下，可以最小化權重平方殘差的總和;看權重最小二乘法 。 理想情況下，每個權重應等于觀察方差的倒數，但是在疊代權重最小二乘算法中，可以在每次疊代時重新計算權重。

L o g i s t i c Logistic Logistic回歸

  對數幾率模型（英語：Logit model，也譯作“邏輯模型”、“評定模型”、“分類評定模型”）是離散選擇法模型之一。

對數幾率分布公式

P ( Y = 1 ∣ X = x ) = e x ′ β 1 + e x ′ β P(Y = 1|X = x) = \frac{e^{x'\beta}}{1 + e^{x'\beta}} P(Y=1∣X=x)=1+ex′βex′β​

  其中參數 β \beta β常用最大似然估計。

  下面介紹與 L o g i s t i c Logistic Logistic回歸相似的一個模型，僅是分布公式上的不同，其算法原理一緻。

拟合一個行星軌道

問題一

  行星遵從橢圓性軌道，在笛卡爾坐标 ( x , y ) (x,y) (x,y)下可以用下列方程來表示：

b 0 + b 1 x + b 2 y + b 3 x y + b 4 y 2 = x 2 b_0 + b_1x + b_2y + b_3xy + b_4y^2 = x^2 b0​+b1​x+b2​y+b3​xy+b4​y2=x2

∙ \bullet ∙ 根據行星軌道在表中10個位置的觀測值，運用最小二乘法來拟合5個參數： b 0 , b 1 , b 2 , b 3 , b 4 b_0,b_1,b_2,b_3,b_4 b0​,b1​,b2​,b3​,b4​，構造出拟合曲線。

x 1.02 0.95 0.87 0.77 0.67 0.56 0.44 0.30 0.16 0.01 y 0.39 0.32 0.27 0.22 0.18 0.15 0.13 0.12 0.13 0.15 \begin{array}{c|lcr} x & 1.02 & 0.95 & 0.87 & 0.77 & 0.67 & 0.56 & 0.44 & 0.30 & 0.16 & 0.01\\ \hline y & 0.39 & 0.32 & 0.27 & 0.22 & 0.18 & 0.15 & 0.13 & 0.12 & 0.13 & 0.15 \\ \end{array} xy​1.020.39​0.950.32​0.870.27​0.770.22​0.670.18​0.560.15​0.440.13​0.300.12​0.160.13​0.010.15​​

在同一幅圖上畫出拟合的橢圓性軌道（連續的實線）及觀測結果，并計算誤差的平方和，評估拟合效果。

思路及解答

  本題是給出了10對 x , y x,y x,y的資料，通過這個資料進行拟合，值得注意的是這個拟合不像普通的線性回歸是：

f ( x ) = b 0 + b 1 x + b 2 x 2 + . . . + b n x n f(x) = b_0 + b_1x + b_2x^2 + ... + b_nx^n f(x)=b0​+b1​x+b2​x2+...+bn​xn

其中 n = 1 , 2 , 3 , . . . n = 1,2,3,... n=1,2,3,...

  如是這種形式，可以直接利用MATLAB中的regress函數進行計算。

  本題的函數是：

b 0 + b 1 x + b 2 y + b 3 x y + b 4 y 2 = x 2 b_0 + b_1x + b_2y + b_3xy + b_4y^2 = x^2 b0​+b1​x+b2​y+b3​xy+b4​y2=x2

  不妨改寫為：

b 0 + b 1 x + b 2 y + b 3 x y + b 4 y 2 − x 2 = 0 b_0 + b_1x + b_2y + b_3xy + b_4y^2 - x^2 = 0 b0​+b1​x+b2​y+b3​xy+b4​y2−x2=0

  然後将變量 x , y x,y x,y看作是已知量，将 b 0 , b 1 , b 2 , b 3 , b 4 b_0,b_1,b_2,b_3,b_4 b0​,b1​,b2​,b3​,b4​看作是變量參數，那麼可以得到一個距離函數（為最小二乘法作準備）：

ϵ ( b 0 , b 1 , b 2 , b 3 , b 4 ) = ∑ i = 1 10 ( b 0 + x b 1 + y b 2 + x y b 3 + y 2 b 4 − x 2 ) 2 \epsilon(b_0,b_1,b_2,b_3,b_4) = \sum_{i=1}^{10} (b_0 + xb_1 + yb_2 + xyb_3 + y^2b_4 - x^2)^2 ϵ(b0​,b1​,b2​,b3​,b4​)=i=1∑10​(b0​+xb1​+yb2​+xyb3​+y2b4​−x2)2

  已知的是 x , y x,y x,y的10組資料，現在要根據已知資料去拟合出 b b b的值，下面就要利用一些高等數學和線性代數的知識（很基礎）。

最小二乘法

  要求得 m i n ∣ ∣ ϵ ∣ ∣ 2 min||\epsilon||^2 min∣∣ϵ∣∣2，即應滿足：

∂ ϵ b 0 = ∂ ϵ b 1 = ∂ ϵ b 2 = ∂ ϵ b 3 = ∂ ϵ b 4 = 0 \frac{\partial\epsilon}{b_0} = \frac{\partial\epsilon}{b_1} = \frac{\partial\epsilon}{b_2} = \frac{\partial\epsilon}{b_3} = \frac{\partial\epsilon}{b_4} = 0 b0​∂ϵ​=b1​∂ϵ​=b2​∂ϵ​=b3​∂ϵ​=b4​∂ϵ​=0

下面以 ∂ ϵ b 4 = 0 \frac{\partial\epsilon}{b_4} = 0 b4​∂ϵ​=0 為例：

0 = ∂ ϵ b 4 = ∂ ( ∑ i = 1 10 ( b 0 + b 1 x + b 2 y + b 3 x y + b 4 y 2 − x 2 ) 2 ) b 4 = 2 ∑ i = 1 10 ( b 0 + b 1 x + b 2 y + b 3 x y + b 4 y 2 − x 2 ) ⋅ y 2 = 2 ∑ i = 1 10 ( b 0 y 2 + b 1 x y 2 + b 2 y 3 + b 3 x y 3 + b 4 y 4 − x 2 y 2 ) 0 = \frac{\partial\epsilon}{b_4} = \frac{\partial(\sum_{i=1}^{10} (b_0 + b_1x + b_2y + b_3xy + b_4y^2 - x^2)^2)}{b_4} \\ =2\sum_{i=1}^{10}(b_0 + b_1x + b_2y + b_3xy + b_4y^2 - x^2)\cdot y^2\\ =2\sum_{i=1}^{10}(b_0y^2 + b_1xy^2 + b_2y^3 + b_3xy^3 + b_4y^4 - x^2y^2) 0=b4​∂ϵ​=b4​∂(∑i=110​(b0​+b1​x+b2​y+b3​xy+b4​y2−x2)2)​=2i=1∑10​(b0​+b1​x+b2​y+b3​xy+b4​y2−x2)⋅y2=2i=1∑10​(b0​y2+b1​xy2+b2​y3+b3​xy3+b4​y4−x2y2)

然後将含有 b 0 , b 1 , b 2 , b 3 , b 4 b_0,b_1,b_2,b_3,b_4 b0​,b1​,b2​,b3​,b4​的項放在左側，其餘放在右側，得到：

∑ b 0 y 2 + ∑ b 1 x y 2 + ∑ b 2 y 3 + ∑ b 3 x y 3 + ∑ b 4 y 4 = ∑ x 2 y 2 \sum b_0y^2 + \sum b_1xy^2 + \sum b_2y^3 + \sum b_3 xy^3 + \sum b_4y^4 = \sum x^2y^2 ∑b0​y2+∑b1​xy2+∑b2​y3+∑b3​xy3+∑b4​y4=∑x2y2

寫成矩陣形式如下：

[ ∑ y 2 ∑ x y 2 ∑ y 3 ∑ x y 3 ∑ y 4 ] [ b 0 b 1 b 2 b 3 b 4 ] = [ ∑ x 2 y 2 ] \begin{bmatrix} \sum y^2 & \sum xy^2 & \sum y^3 & \sum xy^3 & \sum y^4 \\ \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}= \begin{bmatrix} \sum x^2y^2 \end{bmatrix} [∑y2​∑xy2​∑y3​∑xy3​∑y4​]⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=[∑x2y2​]

類似的，得到 ∂ ϵ b 0 \frac{\partial\epsilon}{b_0} b0​∂ϵ​， ∂ ϵ b 1 \frac{\partial\epsilon}{b_1} b1​∂ϵ​， ∂ ϵ b 2 \frac{\partial\epsilon}{b_2} b2​∂ϵ​， ∂ ϵ b 3 \frac{\partial\epsilon}{b_3} b3​∂ϵ​分别如下：

[ ∑ 1 ∑ x ∑ y ∑ x y ∑ y 2 ] [ b 0 b 1 b 2 b 3 b 4 ] = [ ∑ x 2 ] \begin{bmatrix} \sum 1 & \sum x & \sum y & \sum xy & \sum y^2 \\ \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}= \begin{bmatrix} \sum x^2 \end{bmatrix} [∑1​∑x​∑y​∑xy​∑y2​]⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=[∑x2​]

[ ∑ x ∑ x 2 ∑ x y ∑ x 2 y ∑ x y 2 ] [ b 0 b 1 b 2 b 3 b 4 ] = [ ∑ x 3 ] \begin{bmatrix} \sum x & \sum x^2 & \sum xy & \sum x^2y & \sum xy^2 \\ \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}= \begin{bmatrix} \sum x^3 \end{bmatrix} [∑x​∑x2​∑xy​∑x2y​∑xy2​]⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=[∑x3​]

[ ∑ y ∑ x y ∑ y 2 ∑ x y 2 ∑ y 3 ] [ b 0 b 1 b 2 b 3 b 4 ] = [ ∑ x 2 y ] \begin{bmatrix} \sum y & \sum xy & \sum y^2 & \sum xy^2 & \sum y^3 \\ \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}= \begin{bmatrix} \sum x^2y \end{bmatrix} [∑y​∑xy​∑y2​∑xy2​∑y3​]⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=[∑x2y​]

[ ∑ x y ∑ x 2 y ∑ x y 2 ∑ x 2 y 2 ∑ x y 3 ] [ b 0 b 1 b 2 b 3 b 4 ] = [ ∑ x 3 y ] \begin{bmatrix} \sum xy & \sum x^2y & \sum xy^2 & \sum x^2y^2 & \sum xy^3 \\ \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}= \begin{bmatrix} \sum x^3y \end{bmatrix} [∑xy​∑x2y​∑xy2​∑x2y2​∑xy3​]⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=[∑x3y​]

其中， ∑ 1 = 10 \sum1 = 10 ∑1=10.

  将五個矩陣方程合并成一個，用來表示總的方程組，如下：

[ ∑ 1 ∑ x ∑ y ∑ x y ∑ y 2 ∑ x ∑ x 2 ∑ x y ∑ x 2 y ∑ x y 2 ∑ y ∑ x y ∑ y 2 ∑ x y 2 ∑ y 3 ∑ x y ∑ x 2 y ∑ x y 2 ∑ x 2 y 2 ∑ x y 3 ] [ b 0 b 1 b 2 b 3 b 4 ] = [ ∑ x 2 ∑ x 3 ∑ x 2 y ∑ x 3 y ∑ x 2 y 2 ] \begin{bmatrix} \sum 1 & \sum x & \sum y & \sum xy & \sum y^2 \\ \sum x & \sum x^2 & \sum xy & \sum x^2y & \sum xy^2 \\ \sum y & \sum xy & \sum y^2 & \sum xy^2 & \sum y^3 \\ \sum xy & \sum x^2y & \sum xy^2 & \sum x^2y^2 & \sum xy^3 \\ \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}= \begin{bmatrix} \sum x^2 \\ \sum x^3 \\ \sum x^2y \\ \sum x^3y \\ \sum x^2y^2 \end{bmatrix} ⎣⎢⎢⎡​∑1∑x∑y∑xy​∑x∑x2∑xy∑x2y​∑y∑xy∑y2∑xy2​∑xy∑x2y∑xy2∑x2y2​∑y2∑xy2∑y3∑xy3​⎦⎥⎥⎤​⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎡​∑x2∑x3∑x2y∑x3y∑x2y2​⎦⎥⎥⎥⎥⎤​

在程式中是：

M = [10, sum(x), sum(y), sum(x.*y), sum(y.^2);
     sum(x), sum(x.^2), sum(x.*y), sum(x.*x.*y), sum(x.*y.*y);
     sum(y), sum(x.*y), sum(y.^2), sum(x.*y.*y),sum(y.^3);
     sum(x.*y), sum(x.*x.*y), sum(x.*y.*y),sum(x.*x.*y.*y),sum(x.*y.*y.*y);
     sum(y.^2), sum(x.*y.*y), sum(y.^3), sum(x.*y.*y.*y), sum(y.^4)];
N = [sum(x.^2);
     sum(x.^3);
     sum(x.*x.*y);
     sum(x.*x.*x.*y);
     sum(x.*x.*y.*y)];
           

可以寫成：

M [ b 0 b 1 b 2 b 3 b 4 ] = N M \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}=N M⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=N

則參數

[ b 0 b 1 b 2 b 3 b 4 ] = M − 1 ⋅ N \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}=M^{-1}\cdot N ⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=M−1⋅N

即：

b = M \ N;
           

  經過程式設計求解，最後的結果：

[ b 0 b 1 b 2 b 3 b 4 ] = [ − 0.4329 0.5514 3.2229 0.1436 − 2.6356 ] \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}= \begin{bmatrix} -0.4329 \\ 0.5514 \\ 3.2229 \\ 0.1436 \\ -2.6356 \end{bmatrix} ⎣⎢⎢⎢⎢⎡​b0​b1​b2​b3​b4​​⎦⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎡​−0.43290.55143.22290.1436−2.6356​⎦⎥⎥⎥⎥⎤​

S S E = 1.4824 × 1 0 − 4 SSE = 1.4824\times 10^{-4} SSE=1.4824×10−4

問題二

∙ \bullet ∙ 若上述函數表中 x , y x,y x,y加上擾動 Δ x , Δ y \Delta x,\Delta y Δx,Δy，對新的 x , y x,y x,y重新運用最小二乘法計算 b 0 , b 1 , b 2 , b 3 , b 4 b_0,b_1,b_2,b_3,b_4 b0​,b1​,b2​,b3​,b4​值，構造出新的拟合曲線。

解答

  添加擾動和其資料發生了變化,導緻曲線、誤差、拟合效果都發生了變化。

  其中和方差(The sum of squares dueto error)為

S S E = ∑ i = 1 n w i ( y i − y ^ i ) 2 SSE = \sum_{i = 1}^n w_i(y_i - \hat y_i)^2 SSE=i=1∑n​wi​(yi​−y^​i​)2

  SSE越接近于0，說明模型選擇和拟合更好，資料預測也越成功。

  由于我在畫圖的過程中，是周遊 x x x并根據方程計算 y y y，在暴力的求根公式下可能會出現複數根的情況，是以我自己給來 x x x一個範圍，導緻了圖像出現了某種莫名的“缺失”。

  其中我自己編寫了一個測試程式，來粗略的給出 x x x的範圍。

b = [-0.432894270264437;0.551446963140366;3.22294033810576;0.143646182598893;-2.63562548371186];
i = 1;
for x_ = -0.5:0.01:2
    A = b(5);
    B = b(4) * x_ + b(3);
    C = -x_^2 + b(2) * x_ + b(1);
    delta(i)= (B^2 - 4 * A * C)^0.5;
    if isreal(delta(i))
        x_
    end
end
           

變量 x x x_輸出的是求得的根 y y y滿足不是複數的篩選。

完整代碼

  由于本題涉及到繪圖以及矩陣處理問題，我用MATLAB比較順手，故先給出MATLAB的代碼部分，python的代碼以後會貼出。

% main.m
clear
clc
% polyfit
x = [1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.3, 0.16, 0.01];
y = [0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.13, 0.15];

delta_x = [-0.029, 0.0007, -0.0082, -0.0038, -0.0041,...
    0.0026, -0.0001, -0.0058, -0.0005, -0.0034];
delta_y = [-0.0033, 0.0043, 0.0006, 0.0020, 0.0044,...
    0.0009, 0.0028, 0.0034, 0.0059, 0.0024];

[b1, sse1] = my_polyfit(x,y);
[b2, sse2] = my_polyfit(x+delta_x, y+delta_y);

% my_polyfit.m
function [b,sse] = my_polyfit(x,y)
M = [10, sum(x), sum(y), sum(x.*y), sum(y.^2);
     sum(x), sum(x.^2), sum(x.*y), sum(x.*x.*y), sum(x.*y.*y);
     sum(y), sum(x.*y), sum(y.^2), sum(x.*y.*y),sum(y.^3);
     sum(x.*y), sum(x.*x.*y), sum(x.*y.*y),sum(x.*x.*y.*y),sum(x.*y.*y.*y);
     sum(y.^2), sum(x.*y.*y), sum(y.^3), sum(x.*y.*y.*y), sum(y.^4)];

N = [sum(x.^2);
     sum(x.^3);
     sum(x.*x.*y);
     sum(x.*x.*x.*y);
     sum(x.*x.*y.*y)];
 
b = M \ N;

% drow and compare
i = 1;
for x_ = -0.48:0.001:1
    A = b(5);
    B = b(4) * x_ + b(3);
    C = -x_^2 + b(2) * x_ + b(1);
    y_1(i) = (-B + (B^2 - 4 * A * C)^0.5) / (2 * A);
    y_2(i) = (-B - (B^2 - 4 * A * C)^0.5) / (2 * A);
    plot(x_, y_1(i),'.');
    hold on
    plot(x_, y_2(i),'.');
    i = i + 1;
end



for i = 1: 10
    plot(x(i),y(i),'*')
end

sse = 0;
for i = 1:10
    sse = sse + (b(1)+b(2)*x(i)+b(3)*y(i)+b(4)*x(i)*y(i)+b(5)*y(i)*y(i)-x(i)^2)^2;
end
           

該文章首發于 zyairelu.cn

歡迎來到我的網站進行評論及研讨

個人郵箱[email protected]