概述

在统计学中，非线性回归是回归分析的一种形式，其中观测数据由函数建模，该函数是模型参数的非线性组合并且取决于一个或多个独立变量。通过逐次逼近的方法拟合数据。
在非线性回归中，形式的统计模型，
f(x,β)=β1xβ2+x.f(x,\beta) = \frac{\beta_1x}{\beta_2 + x}. f(x,β)=β2+xβ1x.
此函数是非线性的，因为它不能表示为两个 β\betaβ的线性组合。
系统误差可能存在于自变量中，但其处理不在回归分析的范围内。如果自变量不是无差错的，那么这是一个变量误差模型，也在此范围之外。
非线性函数的其他示例包括指数函数，对数函数，三角函数，幂函数，高斯函数和洛伦兹曲线。

回归统计

这个过程的基本假设是模型可以用线性函数近似，即一阶泰勒级数：
f(xi,β)≈f(xi,0)+∑jJijβjf(x_i,\beta) \approx f(x_i,0) + \sum_{j} J_{ij}\beta_j f(xi,β)≈f(xi,0)+j∑Jijβj
其中JIJ=∂f(xi,β)∂βjJ_{IJ} = \frac{\partial f(x_i,\beta)}{\partial \beta_j}JIJ=∂βj∂f(xi,β)，由此得出最小二乘估计量由下式给出。
β^≈(JTJ)−1JTy\hat \beta \approx (J_TJ)^{-1}J^Ty β^≈(JTJ)−1JTy
计算非线性回归统计量并将其用作线性回归统计量，但在公式中使用J代替X. 线性近似将偏差引入统计中。因此，在解释从非线性模型得到的统计数据时，需要比平常更多的谨慎。

普通和加权最小二乘法

最佳拟合曲线通常假定应该看起来平方的总和最小化残差。这是普通的最小二乘（OLS）方法。然而，在因变量不具有恒定方差的情况下，可以最小化加权平方残差的总和;看加权最小二乘法。理想情况下，每个权重应等于观察方差的倒数，但是在迭代加权最小二乘算法中，可以在每次迭代时重新计算权重。

LogisticLogisticLogistic回归

对数几率模型（英语：Logit model，也译作“逻辑模型”、“评定模型”、“分类评定模型”）是离散选择法模型之一。

对数几率分布公式

P(Y=1∣X=x)=ex′β1+ex′βP(Y = 1|X = x) = \frac{e^{x'\beta}}{1 + e^{x'\beta}} P(Y=1∣X=x)=1+ex′βex′β
其中参数β\betaβ常用最大似然估计。
下面介绍与LogisticLogisticLogistic回归相似的一个模型，仅是分布公式上的不同，其算法原理一致。

拟合一个行星轨道

问题一

行星遵从椭圆性轨道，在笛卡尔坐标(x,y)(x,y)(x,y)下可以用下列方程来表示：
b0+b1x+b2y+b3xy+b4y2=x2b_0 + b_1x + b_2y + b_3xy + b_4y^2 = x^2 b0+b1x+b2y+b3xy+b4y2=x2
∙\bullet∙ 根据行星轨道在表中10个位置的观测值，运用最小二乘法来拟合5个参数：b0,b1,b2,b3,b4b_0,b_1,b_2,b_3,b_4b0,b1,b2,b3,b4，构造出拟合曲线。
x1.020.950.870.770.670.560.440.300.160.01y0.390.320.270.220.180.150.130.120.130.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,yx,yx,y的数据，通过这个数据进行拟合，值得注意的是这个拟合不像普通的线性回归是：
f(x)=b0+b1x+b2x2+...+bnxnf(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函数进行计算。

本题的函数是：
b0+b1x+b2y+b3xy+b4y2=x2b_0 + b_1x + b_2y + b_3xy + b_4y^2 = x^2 b0+b1x+b2y+b3xy+b4y2=x2
不妨改写为：
b0+b1x+b2y+b3xy+b4y2−x2=0b_0 + b_1x + b_2y + b_3xy + b_4y^2 - x^2 = 0 b0+b1x+b2y+b3xy+b4y2−x2=0
然后将变量x,yx,yx,y看作是已知量，将b0,b1,b2,b3,b4b_0,b_1,b_2,b_3,b_4b0,b1,b2,b3,b4看作是变量参数，那么可以得到一个距离函数（为最小二乘法作准备）：
ϵ(b0,b1,b2,b3,b4)=∑i=110(b0+xb1+yb2+xyb3+y2b4−x2)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,yx,yx,y的10组数据，现在要根据已知数据去拟合出bbb的值，下面就要利用一些高等数学和线性代数的知识（很基础）。

最小二乘法

要求得 min∣∣ϵ∣∣2min||\epsilon||^2min∣∣ϵ∣∣2，即应满足：
∂ϵb0=∂ϵb1=∂ϵb2=∂ϵb3=∂ϵb4=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
下面以 ∂ϵb4=0\frac{\partial\epsilon}{b_4} = 0b4∂ϵ=0 为例：
0=∂ϵb4=∂(∑i=110(b0+b1x+b2y+b3xy+b4y2−x2)2)b4=2∑i=110(b0+b1x+b2y+b3xy+b4y2−x2)⋅y2=2∑i=110(b0y2+b1xy2+b2y3+b3xy3+b4y4−x2y2)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)
然后将含有b0,b1,b2,b3,b4b_0,b_1,b_2,b_3,b_4b0,b1,b2,b3,b4的项放在左侧，其余放在右侧，得到：
∑b0y2+∑b1xy2+∑b2y3+∑b3xy3+∑b4y4=∑x2y2\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
写成矩阵形式如下：
[∑y2∑xy2∑y3∑xy3∑y4][b0b1b2b3b4]=[∑x2y2]\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]
类似的，得到∂ϵb0\frac{\partial\epsilon}{b_0}b0∂ϵ，∂ϵb1\frac{\partial\epsilon}{b_1}b1∂ϵ，∂ϵb2\frac{\partial\epsilon}{b_2}b2∂ϵ，∂ϵb3\frac{\partial\epsilon}{b_3}b3∂ϵ分别如下：
[∑1∑x∑y∑xy∑y2][b0b1b2b3b4]=[∑x2]\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∑x2∑xy∑x2y∑xy2][b0b1b2b3b4]=[∑x3]\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∑xy∑y2∑xy2∑y3][b0b1b2b3b4]=[∑x2y]\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]
[∑xy∑x2y∑xy2∑x2y2∑xy3][b0b1b2b3b4]=[∑x3y]\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∑xy∑y2∑x∑x2∑xy∑x2y∑xy2∑y∑xy∑y2∑xy2∑y3∑xy∑x2y∑xy2∑x2y2∑xy3][b0b1b2b3b4]=[∑x2∑x3∑x2y∑x3y∑x2y2]\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[b0b1b2b3b4]=NM \begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ b_4 \end{bmatrix}=N M⎣⎢⎢⎢⎢⎡b0b1b2b3b4⎦⎥⎥⎥⎥⎤=N
则参数
[b0b1b2b3b4]=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;

经过编程求解，最后的结果：
[b0b1b2b3b4]=[−0.43290.55143.22290.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⎦⎥⎥⎥⎥⎤
SSE=1.4824×10−4SSE = 1.4824\times 10^{-4} SSE=1.4824×10−4

问题二

∙\bullet∙ 若上述函数表中x,yx,yx,y加上扰动Δx,Δy\Delta x,\Delta yΔx,Δy，对新的x,yx,yx,y重新运用最小二乘法计算b0,b1,b2,b3,b4b_0,b_1,b_2,b_3,b_4b0,b1,b2,b3,b4值，构造出新的拟合曲线。

解答

添加扰动和其数据发生了变化,导致曲线、误差、拟合效果都发生了变化。
其中和方差(The sum of squares dueto error)为
SSE=∑i=1nwi(yi−y^i)2SSE = \sum_{i = 1}^n w_i(y_i - \hat y_i)^2 SSE=i=1∑nwi(yi−y^i)2
SSE越接近于0，说明模型选择和拟合更好，数据预测也越成功。

由于我在画图的过程中，是遍历xxx并根据方程计算yyy，在暴力的求根公式下可能会出现复数根的情况，所以我自己给来xxx一个范围，导致了图像出现了某种莫名的“缺失”。
其中我自己编写了一个测试程序，来粗略的给出xxx的范围。

b = [-0.432894270264437;0.551446963140366;3.22294033810576;0.143646182598893;-2.63562548371186];
i = 1;
for x_ = -0.5:0.01:2A = 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

变量xxx_输出的是求得的根yyy满足不是复数的筛选。

完整代码

由于本题涉及到绘图以及矩阵处理问题，我用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:1A = 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 onplot(x_, y_2(i),'.');i = i + 1;
endfor i = 1: 10plot(x(i),y(i),'*')
endsse = 0;
for i = 1:10sse = 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
欢迎来到我的网站进行评论及研讨
个人邮箱zyairelu@gmail.com

非线性回归模型的原理及评估——解决行星轨道的拟合问题相关推荐

java中数组的内存模型_Java如何在内存有限的情况下实现一个超大数组？jvm性能调优+内存模型+虚拟机原理能解决什么样问题...
在上一篇文章中,小编为您详细介绍了关于<变频器调速问题?三星R458更换CPU>相关知识.本篇中小编将再为您讲解标题Java如何在内存有限的情况下实现一个超大数组?jvm性能调优+内存模型 ...
python 靶心_手把手教你使用Python实战反欺诈模型｜原理+代码
原标题:手把手教你使用Python实战反欺诈模型|原理+代码作者 | 萝卜来源 | 早起Python(ID: zaoqi-python) 本文将基于不平衡数据,使用Python进行反欺诈模型数据 ...
非线性回归模型(part3)--K近邻
学习笔记,仅供参考,有错必纠 PS : 本BLOG采用中英混合模式非线性回归模型 k近邻 The KNN approach simply predicts a new sample using th ...
隐马尔可夫模型 HMM 原理及实现
简介隐马尔可夫模型(Hidden Markov Model,HMM)创立于20世纪70年代.主要用于行为识别,语音识别,文字识别等. 原理简述隐马尔可夫模型由五个部分组成:状态空间S,观测空间O, ...
python非线性回归分析_sklearn实现非线性回归模型
sklearn实现非线性回归模型前言: sklearn实现非线性回归模型的本质是通过线性模型实现非线性模型,如何实现呢?sklearn就是先将非线性模型转换为线性模型,再利用线性模型的算法进行训练模 ...
sklearn模型评选择与评估
sklearn模型评选择与评估 1. 数据集划分 1.1 K折交叉验证 1.1.1 K折交叉验证算法原理 a. 将全部训练及S分成K个不相交的子集,假设S中样本个数为M,那么,每一个子集的训练样本个数 ...
NER依存关系模型：原理，建模及代码实现
NER依存关系模型:原理,建模及代码实现 seq2seq理论基础 NER依存关系建模代码实现命名实体识别(Named Entity Recognization, NER)是AI任务中重要的一类,而 ...
JVM内存模型、原理、垃圾回收、调优
JVM内存模型.原理.垃圾回收.调优,这Java语言的基础,作为Java从业人员是必须要掌握的,另外这也是面试经常会问到的知识. ----------------------------------- ...
面试官：说说Redis之I/O多路复用模型实现原理
你知道的越多,不知道的就越多,业余的像一棵小草! 成功路上并不拥挤,因为坚持的人不多. 编辑:业余草 blog.csdn.net/Seky_fei 推荐:https://www.xttblog.com ...

非线性回归模型的原理及评估——解决行星轨道的拟合问题

非线性回归模型——行星轨道

概述