非线性回归模型——行星轨道
- 概述
- 回归统计
- 普通和加权最小二乘法
- 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β1x.
此函数是非线性的,因为它不能表示为两个 β \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 β^≈(JTJ)−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+b1x+b2y+b3xy+b4y2=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} xy1.020.390.950.320.870.270.770.220.670.180.560.150.440.130.300.120.160.130.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+b1x+b2x2+...+bnxn
其中 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+b1x+b2y+b3xy+b4y2=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+b1x+b2y+b3xy+b4y2−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+b1x+b2y+b3xy+b4y2−x2)2)=2i=1∑10(b0+b1x+b2y+b3xy+b4y2−x2)⋅y2=2i=1∑10(b0y2+b1xy2+b2y3+b3xy3+b4y4−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 ∑b0y2+∑b1xy2+∑b2y3+∑b3xy3+∑b4y4=∑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]⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=[∑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]⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=[∑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]⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=[∑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]⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=[∑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]⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=[∑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⎦⎥⎥⎤⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡∑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⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=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 ⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=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} ⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡−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∑nwi(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]