机器学习-最小二乘拟合
0.前言
解决非线性问题的方法:
最小二乘算法:https://www.mathworks.com/help/optim/ug/least-squares-model-fitting-algorithms.html
这里主要记录Levenberg-Marquardt算法学习笔记!
1.Levenberg-Marquardt算法介绍
基本原理:
计算当C最小时对应的函数参数。
- Levenberg-Marquardt(LM)是一种用于解决非线性最小二乘问题的算法,当通过最小化数据点与函数之间误差的平方和来拟合参数化函数到一组测量数据点时,就会出现最小二乘问题。
- Levenberg-Marquardt曲线拟合方法实际上是两种最小化方法的结合:梯度下降法( gradient
descent)和高斯-牛顿法(Gauss-Newton)。 - 在梯度下降法中,通过在最陡下降方向更新参数来减小误差的平方和。在高斯-牛顿法中,通过假设最小二乘函数是局部二次函数,并求出二次函数的最小值来减小误差的平方和。
- 评估指标——卡方误差chi-squared error criterion:
图片来源:https://msulaiman.org/onewebmedia/LM%20Method%20matlab%20codes%20and%20implementation.pdf
2.Levenberg-Marquardt算法实现-python
最小二乘实现方法:
python实现方法:https://mmas.github.io/least-squares-fitting-numpy-scipy
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)
from scipy.optimize import leastsq
def func(x):return 2*(x-3)**2+1
leastsq(func, 0)
non-linear least squares fit:
def residual(p, x, y):return y - f(x, *p)p0 = [1., 8.]
popt, pcov = optimize.leastsq(residual, p0, args=(x, y))print poptyn = f(xn, *popt)plt.plot(x, y, 'or')
plt.plot(xn, yn)
plt.show()
3.Levenberg-Marquardt算法实现-MATLAB
(1)LM算法原理:https://blog.csdn.net/waitingwinter/article/details/106142276
(2)LM算法实现:Matlab
(3)非线性拟合:MATLAB
rng('default');%将 rand、randi 和 randn 使用的随机数生成器的设置重置为其默认值。保证结果重复。
% 生成数据
len = 1000;
x = linspace(0.1, 300, len)';
y = 300*exp(-x/5) + 10;
% 添加噪声
y_noise = y + 20*randn(len, 1);
% 调用 lsqcurvefit
func = @(p, x) p(1)*exp(-x/p(2)) + p(3);
% 选择合适的初始值
p0 = [100, 2 1];
p = lsqcurvefit(func, p0, x, y_noise)
linespace函数:np.linspace(a,b,c)用于创建一个等差序列的向量,向量值是[a,b]之间均匀分布的c个实数。
拟合结果:
% 拟合结果
y_fit_lsq = func(p, x);
% 调用 nlinfit
pp = nlinfit(x, y_noise, func, p0)
nlinfit函数:beta = nlinfit(X,Y,modelfun,beta0)
使用 modelfun 指定的模型,返回一个向量,其中包含 Y 中的响应对 X 中的预测变量的非线性回归的估计系数。它使用迭代最小二乘估计来估计系数,初始值由 beta0 指定。
绘图:
figure;
plot(x, y_noise, '.', 'LineWidth', 2);
hold on
plot(x, y_fit_lsq, 'LineWidth', 2);
plot(x, y_fit_nlin, 'LineWidth', 2);
hold off
legend('原始数据', '拟合数据-lsqcurvefit', '拟合数据-nlinfit');
结果展示:
LM算法实现:
p1 = LSFittingExpFree(x, y_noise)
fit1 = func(p1, x);
plot(x, fit1, 'LineWidth', 2);%%Levenberg_Marquardt LS Fitting
function result = LSFittingExpFree(x, y)
% Levenberg_Marquardt LS Fitting
% 按列排
x = x(:);
y = y(:);
M = length(x);
% Least Squares Minimization
rng('default');
% 初始值
a = sqrt(rand);
b = sqrt(rand);
c = sqrt(rand);
% 变量个数
nParam = 3;
% 残差
r = y - (a^2 * exp(-x / b^2) + c^2);
% 目标函数
f = r'*r;
% 正则化因子初值
lambda = 1;
% 迭代次数
it = 0;
% 更新标记
updateFlag = true; % 最大迭代次数
maxIter = 100;
% 迭代计算
while it < maxIterit = it + 1;if updateFlagJa = 2 * a * exp(-x / b^2);Jb = 2 * a^2 * x .* exp(-x / b^2) / b^3;Jc = 2 * c * ones(M, 1);J = [Ja, Jb, Jc];g = -2 * J' * r;H = 2 * (J'*J);endHess = H + lambda * eye(nParam);s = -Hess \ g;a1 = a + s(1);b1 = b + s(2);c1 = c + s(3);r1 = y - (a1^2 * exp(-x / b1^2) + c1^2);f1 = r1'*r1;fdr = (f - f1) / f;if fdr > 0a = a1;b = b1;c = c1;f = f1;r = r1;lambda = 0.1 * lambda;updateFlag = true;elselambda = 10 * lambda;updateFlag = false;endif max(abs(s)) < 1e-6disp(['达到收敛条件,迭代次数:' num2str(it)]);break;elseif it == maxIterdisp('达到最高迭代次数,可能不收敛;');end
end
result = [a^2, b^2, c^2];
end
代码来源:https://zhuanlan.zhihu.com/p/388676905
完整代码:lsqcurvefit、LSFittingExpFree、nlinfit、lsqnonlin
% 保证结果可重复
rng('default');
% 生成数据
len = 1000;
x = linspace(0.1, 300, len)';
y = 300*exp(-x/5) + 10;
% 添加噪声
y_noise = y + 20*randn(len, 1);
% 调用 lsqcurvefit
func = @(p, x) p(1)*exp(-x/p(2)) + p(3);
%lsqnonlin
funl = @(p,x,y_noise) p(1)*exp(-x/p(2)) + p(3)- y_noise;
% 选择合适的初始值
p0 = [100, 2 1];
p = lsqcurvefit(func, p0, x, y_noise)
options.Algorithm = 'levenberg-marquardt';
% options = optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective');
[p2,resnorm,residual] = lsqnonlin(funl,p0,[],[],options,x,y_noise)
% 拟合结果
y_fit_lsq = func(p, x);
% 调用 nlinfit
pp = nlinfit(x, y_noise, func, p0)y_fit_nlin = func(pp, x);
p1 = LSFittingExpFree(x,y_noise)
fit1 = func(p1, x);
fit2 = funl(p2,x,y_noise);figure;
plot(x, y_noise, '.', 'LineWidth', 2);
hold on
plot(x, y_fit_lsq, 'LineWidth', 2);
plot(x, fit1, 'LineWidth', 2);
plot(x, y_fit_nlin, 'LineWidth', 2);
plot(x, fit2+y_noise, 'LineWidth', 2);
hold off
legend('原始数据', '拟合数据-lsqcurvefit','拟合数据-LSFittingExpFree','拟合数据-nlinfit','拟合数据-lsqnonlin');%%Levenberg_Marquardt LS Fitting
function result = LSFittingExpFree(x, y)
% Levenberg_Marquardt LS Fitting
% 按列排
x = x(:);
y = y(:);
M = length(x);
% Least Squares Minimization
rng('default');
% 初始值
a = sqrt(rand);
b = sqrt(rand);
c = sqrt(rand);
% 变量个数
nParam = 3;
% 残差
r = y - (a^2 * exp(-x / b^2) + c^2);
% 目标函数
f = r'*r;
% 正则化因子初值
lambda = 1;
% 迭代次数
it = 0;
% 更新标记
updateFlag = true; % 最大迭代次数
maxIter = 100;
% 迭代计算
while it < maxIterit = it + 1;if updateFlagJa = 2 * a * exp(-x / b^2);Jb = 2 * a^2 * x .* exp(-x / b^2) / b^3;Jc = 2 * c * ones(M, 1);J = [Ja, Jb, Jc];g = -2 * J' * r;H = 2 * (J'*J);endHess = H + lambda * eye(nParam);s = -Hess \ g;a1 = a + s(1);b1 = b + s(2);c1 = c + s(3);r1 = y - (a1^2 * exp(-x / b1^2) + c1^2);f1 = r1'*r1;fdr = (f - f1) / f;if fdr > 0a = a1;b = b1;c = c1;f = f1;r = r1;lambda = 0.1 * lambda;updateFlag = true;elselambda = 10 * lambda;updateFlag = false;endif max(abs(s)) < 1e-6disp(['达到收敛条件,迭代次数:' num2str(it)]);break;elseif it == maxIterdisp('达到最高迭代次数,可能不收敛;');end
end
result = [a^2, b^2, c^2];
end
结果对比:
4.Levenberg-Marquardt模型评估
(1)模型评估方法:如lsqnonlin函数中,可以输出resnorm — Squared norm of the residual
或residual — Value of objective function at solution
。
(2)利用卡方检验Chi-Square Test方法来评估模型:
- 卡方检验,也称χ 2检验,是检验两个变量之间有没有关系。
- 卡方检验基本思想:以卡方分布为基础,计算观察值和期望值之间的偏离程度。
卡方检验、t检验和方差分析的区别:https://blog.csdn.net/symoriaty/article/details/100178446
(3)卡方检验(x2检验chi-square test)计算方法:
方法一:四格表资料的x2检验
检验的基本公式为:
A为实际数;T为理论数;
通过查x2值表求P值,判断两者相关性。
方法二:近似法
计算公式:
(4)误差分析(来自论文:The Levenberg-Marquardt method for nonlinear least squares
curve-fitting problems)
参考资料:
论文:The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems
https://msulaiman.org/onewebmedia/LM%20Method%20matlab%20codes%20and%20implementation.pdf
PPT:Data Modeling and Least Squares Fitting 2
https://www.cs.princeton.edu/courses/archive/fall11/cos323/notes/cos323_f11_lecture08_lsq2
PPT:Numerical Optimization using the Levenberg-Marquardt Algorithm
https://mads.lanl.gov/presentations/Leif_LM_presentation_m.pdf
python实现方法-Least squares fitting with Numpy and Scipy:https://mmas.github.io/least-squares-fitting-numpy-scipy
非线性回归中的Levenberg-Marquardt算法理论和代码实现:https://zhuanlan.zhihu.com/p/349953779;代码:https://github.com/manfrezord/MediumArticles/blob/main/LM-Algorithm.ipynb
机器学习-最小二乘拟合相关推荐
- python非线性最小二乘拟合_非线性函数的最小二乘拟合——兼论Jupyter notebook中使用公式 [原创]...
突然有个想法,利用机器学习的基本方法--线性回归方法,来学习一阶rc电路的阶跃响应,从而得到rc电路的结构特征--时间常数τ(即r*c).回答无疑是肯定的,但问题是怎样通过最小二乘法.正规方程,以更多 ...
- 最小二乘拟合,L1、L2正则化约束--转
原文地址:http://blog.csdn.net/u013164528/article/details/45042895 最小二乘法(又称最小平方法)是一种数学优化技术.它通过最小化误差的平方和寻找 ...
- matlab 中最小二乘拟合,matlab 最小二乘拟合
matlab 最小二乘拟合 2018-1-25来自ip:12.12.148.103的网友咨询 浏览量:157 问题补充: matlab 最小二乘拟合 这是m文件 function y=nihe4(p, ...
- 迭代最小二乘拟合椭圆
迭代最小二乘拟合椭圆 //二维点坐标 struct P2d {double x = 0.0;//横坐标double y = 0.0;//纵坐标double angle = 0.0;//角度int nu ...
- 数值分析:数据的最小二乘拟合
1 实验目的 在已知某天在不同时间的前温度高低,借用最小二乘法确定这一天的气温变化规律.通过MATLAB编程,选取适当函数进行求解绘图. 2 实验内容 假定某天的气温变化记录如下表所示: 时间(t) ...
- 数值计算——最小二乘拟合二元一次多项式
数值计算--最小二乘拟合二元一次多项式 最小二乘拟合: 就是根据一系列给定的数据点,求一条曲线使得数据点到曲线的某些(水平.竖直.垂直)距离最短. 推导过程: 1. 设拟合多项式为: 2. ...
- Python之建模数值逼近篇–最小二乘拟合
Python之建模数值逼近篇–最小二乘拟合 介绍 系数ak的确定 函数rk(x)r_k(x)rk(x)的选取 理解和区别 样例 介绍 曲线拟合问题的提法是,已知一组(二维)数据,即平面上的n个点(x ...
- python画抛物线_在python中利用最小二乘拟合二次抛物线函数的方法
1.最小二乘也可以拟合二次函数 我们都知道用最小二乘拟合线性函数没有问题,那么能不能拟合二次函数甚至更高次的函数呢?答案当然是可以的.下面我们就来试试用最小二乘来拟合抛物线形状的的图像. 对于二次函数 ...
- 【多项式最小二乘拟合实验】
一.实验目的 编制以函数{xk}k=0,-,n:为基的多项式最小二乘拟合程序.了解在matlab环境下曲线拟合问题的思想,掌握最小二乘拟合多项式拟合法和曲线拟合法,区分插值于拟合的不同之处. 二.实验 ...
- 最小二乘拟合n阶多项式【Matlab】
最小二乘拟合问题的求解 附录4 最小二乘拟合n阶多项式程序Ⅰ 附录5 最小二乘拟合n阶多项式程序Ⅱ 附录5所示程序运行得到的拟合曲线等结果,与运行附录4所示程序结果相同,也相互验证了附录4.5所示 ...
最新文章
- 给演讲增色的10种简单方法
- Maven问题总结 - 2
- uCOS-II任务建立示例
- 银行系普惠和小贷系普惠,哪个贷款更靠谱?
- pov-inc_yourself劳自己-懒惰的设计师的POV和一些Figma
- Oracle入门(十四.6)之使用标量数据类型
- 输入url到页面返回的过程
- 在Oracle数据库启动时提示没有权限 ora-01031:insufficient privileges
- 【语音判别】基于matlab双门限法判别语音信号【含Matlab源码 1720期】
- 格式工厂 wav 比特率_Mac音乐格式转换工具
- 上海黑马python培训
- python面板数据模型操作步骤_面板数据模型估计一般要做哪些步骤?
- 关于服务器托管,你了解多少?
- 华硕ac68u最佳设置_华硕AC68U路由器APP远程控制设置教程
- css锚点定位不准确问题
- hwclock -w的应用
- 成人世界的规则,越早了解,越早受益
- 如何把立创EDA上导出的原理图和封装导入AD的元件库
- Linux哪个压缩命令可以在window上解压的
- 免格式化转换U盘格式
热门文章
- windows 安装metis_Win10 VS2013 suitesparse-metis-for-windows 1.3.1
- webpack入门+路由配置
- 单片机歌曲代码大全_对于 51 单片机的四大误区!
- verilog设计一个补码加减法运算器_一文搞懂:计算机中为什么用补码来存储数据?...
- 机械自动化算不算计算机相关专业,机械设计制造及其自动化专业属于什么类别...
- CSS:布局——左右两个DIV,左侧宽度固定,右侧占满剩余部分
- CSS:实现跳动小球蒙版效果
- JavaScript:对象转换为字符串、字符串转换为对象
- linux mysql tomcat_Linux下安装Tomcat,Linux下安装Mysql
- 研究自动驾驶技术的算法需要哪些知识?