• 本文是对 Matlab 中 polyfit 函数进行原理解析,并没介绍该如何具体使用 polyfit 函数,具体方法请自行查阅官方文档。
  • 主要涉及数值分析的相关内容。
  • 简单介绍了数据标准化(Z-score 标准化)、QR 分解、Matlab 中的求逆运算符。
  • 相对详细介绍了线性方程组求解的稳定性问题,并引出条件数的定义。
  • 最后根据 polyfit 的源码对它进行计算流程解析,并分析相关的警告该如何处理。

预备知识

标准化

对向量

进行标准化(Standardization,Z-score normalization)定义如下:

其中

表示
的均值,即
表示标准差,即
  • 标准分数(standard score)也叫 Z分数(Z-score);

  • Z分数表示原始数据与平均值的差,单位为标准差,也就时给定值距离均值多少个标准差。

QR 分解定理

阶实非奇异矩阵,则存在正交矩阵
和上三角矩阵
使得
,当
的对角元均为正时,分解唯一。

注:正交矩阵


求逆运算

Y = inv(X) 计算方阵

的逆矩阵。
  • X^(-1) 等效于 inv(X)
  • x = Ab 的计算方式与 x = inv(A)*b 不同,建议用于求解线性方程组。

方程组求解的稳定性

设线性方程组

,其中
有小扰动
,得到解
,那么:

,可得

考虑相对误差:

条件数

  • 的条件数;
  • 表示
    变化时解的相对误差灵敏度的度量;
  • 较小时,解对
    的扰动不敏感。

注:

表示矩阵范数,满足正定性、齐次性以及三角不等式。常见的范数定义有:

polyfit 函数原理

polyfit 使用一维向量

构造具有
列和

m = length(x) 行的范德蒙(Vandermonde)矩阵

(见源码55-59行)并生成线性方程组:

拟合多项式

转换为求解线性方程组(2)(非上述方程组),具体地,

polyfit 先对

进行

QR分解,然后再使用反斜杠运算进行求解(见源码61-69行)。

具体过程如下:

转化为 Matlab 运算即为

p = R(Q'*y)

异常处理逻辑

共以下4种情况

  • MATLAB:polyfit:XYSizeMismatch
  • MATLAB:polyfit:PolyNotUnique
  • MATLAB:polyfit:RepeatedPoints
  • MATLAB:polyfit:RepeatedPointsOrRescale

输入判断

  • polyfit函数只用于拟合(x,y)型二维数据,因此输入数据应是两两对应的,即向量

    长度相等(见源码43-45行)。
  • 若提示MATLAB:polyfit:XYSizeMismatch,请检查上述问题。

欠定方程

  • 若提示MATLAB:polyfit:PolyNotUnique 表示线性方程组

    欠定(矩阵
    的列数大于行数)

考虑

  1. 添加更多的不同的拟合点;
  2. 减少多项式的次数.

条件数过大

由于 Vandermonde 矩阵中的列是向量

的幂,因此

条件数对于高阶拟合来说通常较大,生成一个奇异系数矩阵。在这些情况下,中心化和缩放可改善系统的数值属性以产生更可靠的拟合。而数据的标准化由函数调用时输出情况决定,当且仅当调用[P,S,MU] = POLYFIT(X,Y,N) 才执行数据标准化(见源码50-53行)。

过大时才会提示以下其中一个警告:
  • 若提示 MATLAB:polyfit:RepeatedPoints,考虑添加更多的不同的拟合点
  • 若提示 MATLAB:polyfit:RepeatedPointsOrRescale,考虑1、添加更多的不同的拟合点;2、减少多项式的次数

源码阅读

function [p,S,mu] = polyfit(x,y,n)
%POLYFIT Fit polynomial to data.
%   P = POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
%   degree N that fits the data Y best in a least-squares sense. P is a
%   row vector of length N+1 containing the polynomial coefficients in
%   descending powers, P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).
%
%   [P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and a
%   structure S for use with POLYVAL to obtain error estimates for
%   predictions.  S contains fields for the triangular factor (R) from a QR
%   decomposition of the Vandermonde matrix of X, the degrees of freedom
%   (df), and the norm of the residuals (normr).  If the data Y are random,
%   an estimate of the covariance matrix of P is (Rinv*Rinv')*normr^2/df,
%   where Rinv is the inverse of R.
%
%   [P,S,MU] = POLYFIT(X,Y,N) finds the coefficients of a polynomial in
%   XHAT = (X-MU(1))/MU(2) where MU(1) = MEAN(X) and MU(2) = STD(X). This
%   centering and scaling transformation improves the numerical properties
%   of both the polynomial and the fitting algorithm.
%
%   Warning messages result if N is >= length(X), if X has repeated, or
%   nearly repeated, points, or if X might need centering and scaling.
%
%   Example: simple linear regression with polyfit
%
%     % Fit a polynomial p of degree 1 to the (x,y) data:
%       x = 1:50;
%       y = -0.3*x + 2*randn(1,50);
%       p = polyfit(x,y,1);
%
%     % Evaluate the fitted polynomial p and plot:
%       f = polyval(p,x);
%       plot(x,y,'o',x,f,'-')
%       legend('data','linear fit')
%
%   Class support for inputs X,Y:
%      float: double, single
%
%   See also POLY, POLYVAL, ROOTS, LSCOV.%   Copyright 1984-2018 The MathWorks, Inc.if ~isequal(size(x),size(y))error(message('MATLAB:polyfit:XYSizeMismatch'))
endx = x(:);
y = y(:);if nargout > 2mu = [mean(x); std(x)];x = (x - mu(1))/mu(2);
end% Construct the Vandermonde matrix V = [x.^n ... x.^2 x ones(size(x))]
V(:,n+1) = ones(length(x),1,class(x));
for j = n:-1:1V(:,j) = x.*V(:,j+1);
end% Solve least squares problem p = Vy to get polynomial coefficients p.
[Q,R] = qr(V,0);
oldws = warning('off','all');   % Turn all warnings off before solving
tryp = R(Q'*y);               % Same as p = Vy
catch MEwarning(oldws);             % Restore initial warning statethrow(ME);
end
warning(oldws);                 % Restore initial warning state% Issue warnings.
if size(R,2) > size(R,1)warning(message('MATLAB:polyfit:PolyNotUnique'))
elseif warnIfLargeConditionNumber(R)if nargout > 2warning(message('MATLAB:polyfit:RepeatedPoints'));elsewarning(message('MATLAB:polyfit:RepeatedPointsOrRescale'));end
endif nargout > 1r = y - V*p;% S is a structure containing three elements: the triangular factor% from a QR decomposition of the Vandermonde matrix, the degrees of% freedom and the norm of the residuals.S.R = R;S.df = max(0,length(y) - (n+1));S.normr = norm(r);
endp = p.'; % Polynomial coefficients are row vectors by convention.function flag = warnIfLargeConditionNumber(R)
if isa(R, 'single')flag = (condest(R) > 1e+05);
elseflag = (condest(R) > 1e+10);
end

参考资料

  1. Matlab inv
  2. Matlab @polyfit
  3. Stoer J, Bulirsch R. Introduction to Numerical Analysis[M]. Third Edition. New York, NY: Springer, 2002.
  4. 我的原始语雀文档

java 多项式拟合最多的项数_Matlab polyfit 详解 | 方程组求解的稳定性 | 条件数相关推荐

  1. java 多项式拟合最多的项数_Matlab概率统计与曲线拟合

    一.二项分布 二项分布来源于伯努利试验 (事件发生概率 ) : 含义为独立重复N次试验后, 事件总共发生k次的概率 分布函数 二项分布记为 binopdf 获得事件共发生次的概率 binocdf 为事 ...

  2. java 多项式拟合最多的项数_matlab 多项式拟合EXCEL中复杂数据

    有如下原始数据x,y,它对应的曲线图形为: -9552 -2036.81 -9328 -2025.62 -9168 -2014.43 -9024 -2003.25 -8928 -1992.06 -88 ...

  3. java 多项式拟合最多的项数_MATLAB绘制带置信区间的拟合曲线

    曲线拟合是已知离散点上的数据集,构造一个解析函数(其图形为一曲线),使在原离散点上尽可能接近给定的值.MATLAB中与曲线拟合有关的函数主要有polyfit.polyval和polyconf. 01p ...

  4. java 多项式拟合最多的项数_机器学习(1)--线性回归和多项式拟合

    机器学习(1)--线性回归和多项式拟合 机器学习(2)逻辑回归 (数学推导及代码实现) 机器学习(3)softmax实现Fashion-MNIST分类 一 线性回归 线性回归,顾名思义是利用线性模型对 ...

  5. java 多项式拟合最多的项数_python实现2019nCoV确诊数据拟合与预测

    在获取2019-nCoV疫情实时追踪数据后,接下来就要着手探索这些数据的规律和内在联系了,本期主要介绍如何利用logistic模型来拟合和预测确诊人数.主要分为下面三个步骤来进行第一步,画出现有数据的 ...

  6. java 多项式拟合最多的项数_牛顿插值法、曲线拟合、多项式拟合

    2020年10月4日研究了一下牛顿插值法,其用途是使用x.y两组数值,根据新的x值返回对应的y值,与TREND.FORECAST函数不同,这种方法可应对非线性数据.其作用类似于图表中的曲线拟合或LIN ...

  7. 如何确定matlab多项式拟合的阶数,基于多项式拟合函数趋势项与阶数估计加速度、速度、位移的方法与流程...

    本发明属于信号处理领域,尤其涉及一种基于多项式拟合函数趋势项与阶数估计加速度.速度.位移的方法. 背景技术: 目前信号处理领域常用的加速度积分方法主要有时域积分和频域积分两种.时域积分常数项经积分会产 ...

  8. 不愧是京东大牛!用Java实现黄金分割数的示例详解(附代码)

    这篇文章主要介绍了java 实现黄金分割数的示例详解,具有很好的参考价值,希望对大家有所帮助.一起跟随小编过来看看吧. 黄金分割数 0.618 与美学有重要的关系.舞台上报幕员所站的位置大约就是舞台宽 ...

  9. 农夫过河算法java,Java农夫过河问题的继承与多态实现详解

    Java农夫过河问题的继承与多态实现详解 发布时间:2020-08-22 06:04:29 来源:脚本之家 阅读:61 作者:小任性嘛 题目描述: 一个农夫带着一匹狼.一只羊.一颗白菜要过河,只有一条 ...

  10. Java中print、printf、println的区别 详解

    Java中print.printf.println的区别详解 printf主要是继承了C语言的printf的一些特性,可以进行格式化输出 print就是一般的标准输出,但是不换行 println和pr ...

最新文章

  1. iptables配置-Linux系统安全防火墙
  2. 【selenium 3】 Mac 下测试环境搭建 Firefox 47+ gecko driver Mac
  3. 互联网或将进入泡沫2.0时代
  4. [Swift]LeetCode522. 最长特殊序列 II | Longest Uncommon Subsequence II
  5. 数学建模——TOPSIS综合评价模型Python代码
  6. java byte 判断相等_转发收藏 | 史上最全Java面试题+面试网站推荐!(含答案)
  7. MySQL中变量的定义和变量的赋值使用(转)
  8. 使用HTML5、CSS3和jQuery增强网站用户体验
  9. Hadoop的学习路线图
  10. SpringBoot笔记整理(二)
  11. ubuntu java 1.6 安装,ubuntu 中安装java jdk 1.6
  12. php面试php数组变ahp,关于PHP字符串的一道面试题
  13. Servlet实现图片读取显示
  14. tp5ajax即点即改,TP5中即点即改,json分页,单删
  15. java音乐网站源码_Vue + SpringBoot + MyBatis 音乐网站
  16. 2017.10.17笔记
  17. html5 天地图,天地图API
  18. 先面对现实,再寻找理想
  19. echart 地图 某个地区_一个让echarts中国地图包含省市轮廓的技巧
  20. 机器人学习--移动机器人定位导航性能评估规范

热门文章

  1. [Bochs]Bochs调试技术
  2. Tomcat优化之配置线程池高并发连接
  3. 39. Element compareDocumentPosition() 方法
  4. css3中的过度transition与动画animation以及字体@font-face
  5. 如何遍历JTree的每一个节点
  6. [AtCoder Beginner Contest 133]F - Colorful Tree
  7. 工具使用-curl/wget
  8. layui select下拉框选项不显示
  9. 00048_this关键字
  10. Linux环境下配置JDK,java环境