图片缩略图

polyfitn

版本1.3 (56.7 KB) 约翰D 'Errico
1维或n维多项式建模
4.9
78年评级

235下载

更新2016年4月27日

查看版本历史

查看许可协议

Polyfitn是polyfit的扩展,允许用户创建具有多个自变量的模型。它还允许用户指定一个通用模型,例如,一个二次模型,有常数项和二次项,但没有线性项。
例如,为了将多项式模型拟合到从余弦曲线中选择的点,我们只需要偶数有序项。
X = -2:.1:2;
Y = cos(x);
P = polyfitn(x,y,'constant x^2 x^4 x^6');
p.Coefficients
ans =
[0.99996 -0.49968 0.041242 -0.0012079]
当然,系数不会是精确的,因为我只使用了有限数量的项,这本质上是一个截断的泰勒级数,而且我只有有限数量的点来构建模型。余弦级数的前4个系数是:
>> [1 -1/2 1/24 -1/720]
ans =
1 -0.5 0.041667 -0.0013889

我们得到了预期的结果。

当然,polyfitn适用于高维空间,因为它的设计初衷就是要解决这个问题。

X =兰特(100,1);
Y = rand(100,1);
Z = exp(x+y) + randn(100,1)/100;
P = polyfitn([x,y],z,3);

结果可以转换为符号形式,以便更简单地查看模型。这里我将使用我的sympoly工具箱,但也提供了一个polyn2sym函数。

polyn2sympoly (p)
ans =
0.43896* x1 ^3 + 1.4919* x1 ^2* x2 + 0.041084* x1 ^2 + 1.4615* x1 * x2 ^2 - 0.095977* x1 * x2 + 1.2799* x1 + 0.56912* x2 ^3 - 0.15306* x2 ^2 + 1.361* x2 + 0.94819

当然,参数误差估计是为那些希望确定所生成的项的重要性的人生成的。

我还提供了用于评估这些模型和区分模型的工具。

警告-小心使用高阶多项式来拟合你的数据。仅仅因为低阶模型有效,高阶模型就不一定更好。高阶多项式在数据点之间经常遭受严重的振铃。总是绘制你的数据。考虑一下您将要构建的模型。然后绘制得到的模型。用你的眼睛来验证结果,而不仅仅是看一个r平方系数(尽管我也会返回这个参数)。

如果你确实发现高阶多项式模式是必要的,因为你的曲线太复杂了,考虑使用回归或平滑样条模型代替。

引用作为

John D'Errico(2021)。polyfitn(//www.tatmou.com/matlabcentral/fileexchange/34765-polyfitn), MATLAB中央文件交换。检索

意见及评分(135

埃Parodi

你好,功能很好,谢谢!
不过我有个小问题…当你使用polyn2sympoly(p)提取的系数时,这些系数有几个小数,使用它们的拟合不会很好。我怎么解决这个问题呢?请

Jihwan嘘

杰克Souweha

Daniel H

当两个自变量的大小不同时,是否有可能使它工作?我有一个多维数组,我想适合,但出于各种目的,我的自变量向量有不同的长度!

托马斯gerrit

托马斯gerrit

非常好的工具箱!我有一个问题;有没有办法为每个indepvar的ModelTerms指定一个特定的顺序?假设我将函数调用为:
p = polyfitn([X(:),Y(:),Q(:)], Z(:), 2);
p.ModelTerms变成:
2 0 0
11 10 0
10 0 1
1 0 0
0 2 0
0 1 1
0 1 0
0 0 2
0 0 1
0 0 0

但如果我想让Y和Q的顺序都是1而不是2,我将不得不手动构建modelTerms矩阵等于:
2 0 0
11 10 0
10 0 1
1 0 0
0 1 1
0 1 0
0 0 1
0 0 0

有没有办法让函数自动解析这个?

林王

当我尝试用三个输入拟合多项式模型时,我得到了这个错误“错误使用betainc
Z必须是实数
非负。

在polyfitn (line
266)
polymodel。p =
betainc (polymodel.DoF. / (t ^ 2
+
polymodel.DoF)、polymodel.DoF / 2 1/2);”。但实际上,所有的输入都是非负的实数。有人知道哪里出了问题,怎么解决吗?

Nikolaos Efkarpidis

蒂姆

托马斯。

Sushant Powar

John D'Errico的精彩分享。但有一个小问题,大多数时候当我们使用高阶多项式来建模大型数据集时,自变量矩阵的列变得线性相关,矩阵变得奇异。
i.在polyfitn的:232-241行中,我们可以用伪逆代替法线吗?
2但是,这会导致一些系数的方差变为零,t统计量变为NaN,如何克服这个问题?

Yihan张

看起来很棒,似乎非常准确(与我手动计算的几个相比)。但这只是很好的包装,你得到所有你需要在一个点击!伟大的工作!

Heehang金

托拜厄斯Kistler

阿里sameer

亲爱的,约翰D 'Errico
感谢您分享您的有益代码
你能解释一下你用它来求项数的方程吗?实际上,我根据下面的公式计算了二次模型的项(系数)的数量

No = (n+2)(n+1)/ 2% n为变量数
目前,如果我使用3度或以上,我在计算系数数量时面临一个问题。你能帮帮我吗?

问候

莱奥波尔多巴西

尤金Pashinov

有人知道如何在gpu或显著加速计算上运行这个函数吗?

西蒙•施密德

谢谢你的工作。你知道如何使用多项式作为被接受的优化函数(例如fminimax)吗?多变量的输入必须是向量化的。

Suhail Najm Abdullah

它只是工作得很好。拇指。

约翰D 'Errico

术语[11 11 10 0]是一个三阶项,至少在polyfitn看来是这样的。这是一个三方交互的术语。如果你想要在你的模型中使用这个项,你总是可以指定这个模型。

Andreas Erkinger

阿图罗冈萨雷斯

伟大的工作!我只有一个问题。当我看modelterms是如何生成的,让我们说:
Order = 2;
P = 4;
我们得到:
modelterms =
2 0 0 0
1 1 0 0
1 0 1 0
1 0 0 1
1 0 0 0
0 2 0 0
0 1 1 0
0 1 1
0 1 0 0
0 0 2 0
0 0 1 1
0 0 1 0
0 0 0 2
0 0 0 1
0 0 0 0
varlist = {'X1'} {'X2'} {'X3'} {'X4'}
我的问题是为什么有些组合缺失了,比如:
1 1 1 0
1 0 1 1

阿图罗冈萨雷斯

丽萃Hemsley

卢卡斯Theisen

我仍然是这个套餐的忠实粉丝!

康巴罗斯

嗨,约翰!了不起的工作!我是MatLab的初学者,我试着用你们的工具来得到曲线的方程。我有两组点,x和y,其中命令图(x,y)产生所需的曲线(https://imgur.com/a/Nq7pQXA)。我想要的是这条曲线的方程。可以用你们的工具来做吗,还是我全错了?非常感谢:)

塞德里克Devivier

请在下面的地址找到一个函数的修正版本,该函数采用Massimo Ciacci解决方案,并在系统确定时修复betainc错误(Ndata == Nterms)。
http://ecco-qua.modexpeng.com/?fi=IVAp1nE6iR&fn=polyfitn.m

希望能帮助到别人

马西莫Ciacci

在第419行发现了一个简单的错误:
------------
如果没有这个修复程序,它会触发一个解析错误,其中包含“Y^2*X^2”这样的术语。
K > > pTot = polyfitn ([YY (:), XX (:)), Zfunc (:), Y Y ^ 2 Y * X ^ 2 Y ^ 2 * X ^ 2”)
使用polyfitn>解析模型时出现错误(第379行)
变量不是一个有效的名称:'2'
------------
解决办法:
Term = Term (1+k1:end);修复解析错误
%term(k1(1)) = ' ';

------------

通过这个修正,指数部分也从术语字符串中消耗,而不会留下类似
“X ^ 2 * X ^ 2”
“^ 2 * X ^ 2”
' 2x ^2' <——错
'X^2' <——对

马西莫Ciacci

这真是太棒了!这招很管用。对于那些有一维向量x y和矩阵Z的人,这样做可能会有帮助
(X, Y) = ndgrid (X, Y);
p = polyfitn ([X (:), Y(:))、Z(:),“恒X1 X1 ^ 2 X2 ^ X1 * X2”)%
感谢这个伟大的工具!

桑托什

ZuzkaT

请问如何解决“betainc”错误?

杰夫·拉皮埃尔

发现这个工具非常有用!
幂律也适用吗?如果是,它的语法是什么?提前谢谢!

Catsui

陈显乐

劳伦Bocognano

不好意思,刚发完消息我终于找到了“下载”…——“无论如何,还是非常感谢你

亚历山大Jelzow

在这里,在文件交换中,我找到了“permn”工具。
permn (0:2, 2)
给出了一个完整二次模型的全部9项。

亚历山大Jelzow

伟大的工具!
然而,我注意到如果我使用polyfitn([x1 x2], y, 2),模型项不包含混合项[2 1]和[1 2]。只使用混合术语[11 1]。这是有意的吗?还是一个bug?

何塞•巴尔加斯

本杰明。博得纳

@Jacob如果你的数据被组织在一个10x10的网格中,你的x, y和z向量每个都有100个值。X和y将有每个不同的网格值重复10次。这与使用scatter3绘图函数的格式相同。这样做更有用,因为它不需要将x和y值组织成网格。

雅各Mannhardt

@Chad:
谢谢乍得!是的,我最终弄明白了,但我仍然不明白,当z的大小不是n_x x n_y时,它怎么可能是x - y平面上的数据分布?

乍得Stucki

@Jacob:
polyfitn的语法要求在polyfitn(x,y,z…)是列向量或行向量,而不是网格。

雅各Mannhardt

嘿,约翰,
谢谢这个功能!我有一个问题:由于我在Matlab中绝对缺乏知识,我无法理解以下内容:我有一个(x,y)平面/网格-让我们说10x10 -和z=f(x,y)的测量数据点。z是一个有100个元素的10x10矩阵。每当我尝试polyfitn([x,y], z, 3) Matlab告诉我z必须是一个向量,而不是一个矩阵。你的“手册”说z的大小必须是nx1,其中n是x或y的大小。但我不明白为什么z只能有10个条目,而不是10^2。如果你能帮我就太好了。雅各干杯

Umberto Maria Ciucani

非常感谢!

甘伊萨克

秦许

托拜厄斯Eberhorn

嗨,约翰,
谢谢你的工具箱,这是一个很大的帮助!

我正在处理一个有4-5个输入参数(从excel文件导入,一个列向量)和其他参数依赖的一个参数的问题。
我用带有两个参数的CF工具箱验证了我的函数。
当我尝试使用第三个输入参数时,我得到了一个警告:

警告:矩阵接近奇异或缩放不良。结果可能不准确。RCOND = 1.794064e-17。

此外,我检查了R2(=0.9495),结果在我看来调整得太好了。

你能想出一个可能的方法来解决这个问题吗还是我遗漏了什么?
我以为工具箱是专门为处理多维问题而设计的。

如果你能想到任何可以帮助我的东西,我将非常感激!

非常感谢,
欢呼声托拜厄斯

这是我的代码。

x = Durchmesser3;
y = Plattendicke3;
w = EModulNiet3;
z = kraft t3;

P = polyfit([x,y,w],z,2);

如果存在('sympoly') == 2
polyn2sympoly (p)
结束
如果存在('sym') == 2
polyn2sym (p)
结束

约翰D 'Errico

@Ramon:

A + b*√(x):
polyfitn (x, y) (0; 1/2))

A + b/x:
polyfitn (x, y, [0, 1])

Jaime

约翰! !谢谢。不用再说了,我要给它编程。谢谢!;)

雷蒙Dalmau

它是多拟合的,能够拟合这样的函数:

A + b*√(x)或者A + b/x

?

莫里茨踢

非常感谢这个工具箱,非常有用!

Andrejs Fedjajevs

它是否能够替代统计工具箱中的逐步eglm功能?

塞巴斯蒂安·巴哈蒙德

松鸦

是否有具体的会议论文解释了polyfitn的整个过程?
这将有助于深入研究算法,而不是试图理解编码。
谢谢

Per-Johan

苏西阿曼

有人有在3d中如何使用这个的例子吗?

谢谢!

直到Huelnhagen

我正在寻找一种方法来拟合三维领域使用多项式。这个函数完美而快速地解决了任务。

谢谢分享!

丹尼尔Cipriano

嗨,约翰!我是Matlab的初学者。我需要从一个有x和y两个向量的曲面找到一个多项式拟合。

X是1x191 y是1x51。

然后是平面z,它的尺寸是51x191。

我可以用:
表面(x, y, z)

但是我不明白我怎样才能用你的文件来找到我的方程。

谢谢你的回答!

AndBu

罗伯特·威利

我已经克服了这个缺点,成功地使用了它。完整的魔法!谢谢你!

罗伯特·威利

嗨,约翰,
我可以开始说,我是一个严重的matlab新手。我有三组记录的数据,我用scatteredInterpolant和meshgrid创建一个曲面图。现在我想求出这个曲面的公式。但是当输入三个数据集时,我得到“错误使用betainc
Z必须是实数
非负。

在polyfitn (line
266)
polymodel。p =
betainc (polymodel.DoF. / (t ^ 2
+
polymodel.DoF)、polymodel.DoF / 2 1/2);“

Z变量都是25到26的正数。我不确定这是否是使用三个特定数据集而不是从函数中获取数据的问题。

多谢。

塞Scanziani

约翰D 'Errico

Mimo -我不知道你在问什么。

你是在问是否只有一个变量x,但有多个y向量?Polyfitn不是为了解决这种形式的问题而设置的,但它只是Polyfitn内部的一个循环。我想我可以写polyfitn来遍历y的列。这是我下次发布版本时的一个想法。:)

问题是你是否使用polyfitn来简单地计算回归系数,或者你是否想要整个9码,即标准误差,R^2等。

例如,如果你的问题像你描述的那样简单……

假设x是数据的一个nx1列向量,自变量向量,对所有问题都是固定的。Y是一个NXP数组,所以有p个独立的问题需要解决。那么回归系数很容易得到:

Ab = [ones(n,1),x(:)]\y;

数组ab将是一个2xp数组,其中ab的每一列都是一组回归系数。所以ab的第一行是每个问题的常数项。

A = ab(1,:);
B = ab(2,:);

同样,这没有提供任何统计数据,R^2等。但这是一个简单的解决方案,我认为你提出的问题。

米姆

嗨,约翰,
这个函数很神奇。我刚刚遇到了一个问题,是否也有可能找到一个多元线性回归,形式为“y(I)=a(I)+b(I)”。X”,只有一个变量“X”,通过找到向量a(i)和b(i)?

约翰D 'Errico

假设和是变量,那么就这样做:

MDL = polyfitn([phi(:),theta(:)],z(:),2);

您可能需要使用meshgrid或ndgrid(视情况而定)来生成phi和theta。

嗨,伙计们!我一直在三维图中工作,(从- /3到/3的点)和相同,每个点都有一个值(z),所以最后我得到了和的每个可能组合的矩阵。我想用和的函数多项式地近似矩阵的点。我已经看了多fitn但你必须放的点(z)不可能是一个矩阵。Polyfit也不行因为它只适用于2D (,z)或(teta,z)

乔治Varverakis

不错的工作!

约翰D 'Errico

这里没有超出基本的东西。(我承认,基本的东西是因人而异的。)就是经典的多元线性回归。可以参考Draper和Smith关于回归分析的文章。所以stats工具箱中的任何工具都可以做同样的事情,只要你为你的模型正确地准备了设计矩阵。

亲爱的约翰

谢谢你的好工具。我正在积极地使用它,它工作得非常好。我想学习你使用的方法。我想知道你是否知道我在哪里可以找到相关的理论。

欢呼,

塞巴斯蒂安。

netorha

你好约翰,
非常感谢你的精彩剧本。它工作得非常非常好!

但是,我有一个关于一个有点难的问题的问题…实际上,我通过使用matlab fit()与以下代码执行数据拟合(其中arr_node_coords_y, _x和arr_temperaturefield是33x5-Matrices):

[xData, yData, zData] = prepareSurfaceData(arr_node_coords_y, arr_node_coords_x, arr_temperaturefield);

设置fit类型和选项。
Ft = fittype('poly24');
Opts = fitoptions(ft);

拟合模型到数据。
[fitresult, gof] = fit([xData, yData], zData, ft, opts);

%的阴谋
Subplot (2,2,2);
h = plot(fitresult, [xData, yData], zData);
轴([0 2.5 0 8 20 60]);
xlabel('Width [m]');
ylabel('高度[m]');
zlabel('温度[°C]');
标题(“合身的温度”);
网格;
Daspect ([8 8 60]);
colorbar;

...怎么可能用polyfitn()执行同样的事情?...-我尝试了几种方法,但不幸的是,我没有找到一个有效的解决方案。

也许有人能帮我…?

rihab

Balasundaram汉

你好,约翰
自变量是100x1100x1因变量是100x100也是复数。你能帮我做一下多项式拟合吗?
而我一直试图做多边形对角线元素的因变量,它显示'错误使用betainc
输入必须是实数、满值和双精度或
单”。
这意味着什么?我该如何克服这个问题?

约翰D 'Errico

@Sagar -但这不是一个多项式模型。所以你不能使用多项式模型的工具来拟合任何一般的非线性模型。使用非线性回归工具。曲线拟合工具箱是最简单的选择,但对于初学者来说,在统计或优化工具箱中还有许多其他选择。或者使用像fminsearch或fminbod这样的优化器,但你需要最小化残差的平方和。

Sagar

嗨,功能很好。我知道我们可以这样做:poly_fit = polyfitn(x, y, 3);用于测试输入变量的所有组合(例如x1和x2)。但是如果我想构造自定义的输入指数,我该怎么做呢?我想拟合这样的东西:
Poly_fit = 2.3x1^3 + 1.5e^-x2

我怎么做呢?

阿卡什

约翰D 'Errico

修正了NaN问题。很抱歉。

乍得格林

John D'Errico的众多优秀作品之一。

如果因变量包含任何nan,则会产生令人困惑的betainc错误。我建议在代码早期添加这个错误检查:

assert(sum(isnan(depvar))==0,'因变量不能包含任何nan .')

卢卡斯Theisen

JinSun

难以置信的谢谢你的贡献!

阿卡什

工作很好,应该包括在标准包装

Christoph

嗨,约翰,

我使用“polyfitn”来寻找基于三个输入变量的非线性系统的变化参数的多项式。插入这些多项式在三个方程,然后我尝试“解”他们,但我不能收敛到一个解决方案插入不同阶的多项式。

我一直在查看你的文件交换列表,寻找其他相关代码。我想知道我是否应该尝试使用一些样条方法而不是“polyfitn”,但我不确定我如何能将它插入方程并解决它们。同时求解三个非线性方程的其它方法也是如此。

谢谢

约翰

约翰D 'Errico

抱歉,但我只能假设Gautam和Mitch一样,使用的是一个相当老的MATLAB版本。

如果你真的想解决这个问题,你应该提供足够的信息让我知道是不是这样。由于您所描述的测试用例工作得非常完美,正如我在这里展示的那样,我的建议是使用一个合理的最近的版本,因为Mitch说他使用的是一个7年前的版本,并且得到了同样的错误,所以我假设情况就是这样。

如果您确实在使用最新的版本,并且仍然存在问题,那么我只能假设您使用不当,或者您安装了其他导致问题的东西。例如,当我执行这个测试时,它运行得很好,但我没有7年前的MATLAB版本,我甚至不能运行那么旧的东西。

Indvar = rand(12,2);
Depvar = rand(12,1);
M = [0 0;1 0;0 1;2 0;1 1;0 2;3 0;2 1;1 - 2;0 3];
P = polyfitn(indvar,depvar,m)
p =

ModelTerms: [10x2 double]
系数:[-0.0610 -4.6966 5.5355 27.4400 -24.9658 -0.4228 -32.4825 42.3187 -18.4250 5.9917]
参数var: [12.5177 249.5518 129.5089 850.2142 511.0641 282.0906 504.6530 953.0042 477.2848 209.7347]
ParameterStd: [3.5380 15.7972 11.3802 29.1584 22.6067 16.7956 22.4645 30.8708 21.8468 14.4822]
景深:2
P: [0.9878 0.7943 0.6748 0.4460 0.3845 0.9822 0.2851 0.3040 0.4878 0.7192]
R2: 0.6497
AdjustedR2: -0.9265
RMSE: 0.1327
VarNames: {'' ''}

没有你提供的有用信息,我猜不出是什么问题。我在这里提到的信息包括告诉我你有什么MATLAB版本。这将包括告诉我完整的错误消息,而不仅仅是错误消息的片段。

顺便说一下,试图将一个只有12个数据点的数据集与一个有10个项的模型拟合,将会产生一个几乎毫无意义的结果。我认为这是数据过度拟合的一个例子。

Gautam

D'Errico先生,非常感谢您的建议。它看起来很有帮助,除了当我给出一个12 × 2的矩阵作为自变量矩阵和一个12 × 1的列矩阵作为因变量时,它显示了以下结果:
未定义函数或变量“varlist”。

polyfitn错误(第233行)
polymodel。VarNames = varlist;

我使用的模型是这样的:
M = [0 0;1 0;0 1;2 0;1 1;0 2;3 0;2 1;1 - 2;0 3];

任何关于如何让你的功能为我工作的建议将是非常有用的,非常感激。非常感谢。
Gautam。

约翰D 'Errico

米奇,

您认为它在您的测试中被正确定义是不相关的。我只是不知道它的定义是什么,因为你给了我所有的东西,除了这个变量。所以除非你能给我一个失败的例子,否则我不知道你在做什么。我使用不同的变量名是完全无关的。我可以很容易地将传入的变量命名为Fred和Barney,它就可以运行了。

我也不知道旧的MATLAB版本是否有问题。该网站告诉我,我使用R2007b编写了这个工具,所以它应该适合你,但也可以想象,多年来对它进行的某个mod现在使它在旧版本中失败了。这是我无法测试的东西,因为我不能使用旧的MATLAB版本。

所以我的第一个建议是运行我给出的测试。验证它是否正常运行。如果它没有运行,那么您需要看看您的测试中的输入参数有什么不同。或者你可以把你的资料发给我,如果你不想在网上公布的话。

如果问题是你的旧版本确实是问题所在,那么你将需要学习使用反斜杠(不是那么困难的任务)来解决问题,或者可能获得一个stats工具箱的副本。

嗨,约翰,
谢谢你的快速回放。
我很抱歉在我的问题中缺少数据,但我的变量err_func是你的outpar,它是定义的。
我检查了所有变量的尺寸,他们是正确的。
然后我尝试复制和过去你的例子,出现同样的错误:

???未定义函数或变量“varlist”。

在==> polyfitn的232错误
polymodel。VarNames = varlist;

可能是我的Matlab版本(R2007b)的问题?

提前谢谢你
米奇

约翰D 'Errico

米奇,

我不知道你在做什么。我最初的猜测是您试图从编辑器或类似的地方“运行”该函数。您不会运行这样的函数,尽管许多人尝试这样做。

我确实看到你在你的例子中没有定义变量err_func,所以我甚至不能真正猜出你的错误是什么。

Inpar = rand(20,2);
超额=兰特(20,1);
Model_exp = [1 0;2 0;0 1;1 1;0 2];
P = polyfitn(inpar,outpar,model_exp);

这个测试用例运行得很好。

嗨,约翰,
谢谢你的投稿。
我试着用它来处理我自己的数据。但是我不能正确设置modelterms参数:

即。

in_par =(ζθ);
Model_exp = [1 0;2 0;0 1;1 1;0 2];
P = polyfitn(in_par,err_func,model_exp)

运行后出现以下错误:

???未定义函数或变量“varlist”。

在==> polyfitn的232错误
polymodel。VarNames = varlist;

==> launch_regr在14时出现错误
P = polyfitn(in_par,err_func,model_exp)

你能帮帮我吗?

弗朗西斯科

约翰D 'Errico

我想这仍然回避了一个问题,polyfitn是否需要有能力拟合多个右手边,即同一个模型,但有多个因变量?

就我个人而言,我发现这种功能很少需要,而且当我们在我以前的工作中使用这种功能支持一大组颜色科学家时,它实际上只对构建监测模型有用。金宝app从那时起,监控器变得更加复杂,因此不那么线性。我猜想,使用更复杂的模型可以更好地处理较新的显示器。对于打印机等来说,这当然是真的,最好使用3d查找表来处理,就像我自己的gridfit工具一样。

但事实是,我经常惊讶地发现其他人以我意想不到的方式使用这些工具。

因此,如果我看到至少有一些人要求这种功能,我会考虑添加它。在此之前,我的偏好是在这方面代码的简单性。(嘿,我应该退休了!)

马克Mikofski

我完全同意约翰的观点;我没有仔细考虑,就轻率地提出了那个建议。如果我把你引入歧途,我很抱歉。感谢约翰的建议和高超的数学工具——它们是无价的!

约翰D 'Errico

我很抱歉,但是我完全不同意马克给安德森的建议。

具有多个输出变量的系统只需要一个简单的polyfitn循环。该循环的成本根本不足以增加代码的复杂性和接口的复杂性。

建议人们使用诸如Gauss-Newton、lsqnonlin或lsqcurvefit等迭代方案,这些都是荒谬的方法。你仍然需要循环,或者你需要生成一个大问题本质上是一个块对角雅可比矩阵。所有这些都是为了避免一个简单的循环?或者你会写一个循环,只是为了避免循环?注意,在Mark的建议中,您将需要编写一些代码来评估模型。如果你使用lsqnonlin或lsqcurvefit这样的工具,那么你要么需要提供雅可比矩阵,要么你需要让这个工具多次调用你的目标函数来计算有限差分来生成雅可比矩阵。

不要为了避免简单的循环而使用迭代工具来解决真正可以用基本线性代数解决的问题。请注意,即使迭代工具将快速收敛,这些工具实际上将尝试进行多余的迭代,直到它有效地意识到它确实收敛了。这需要一个备用的雅可比矩阵求值。

马克Mikofski

@Anderson,你可以用newtonraphson (//www.tatmou.com/matlabcentral/fileexchange/43097-newton-raphson-solver)或者如果你有优化工具箱lsq曲率或lsqnonlin (//www.tatmou.com/help/optim/nonlinear-least-squares-curve-fitting.html

安德森

这个工具看起来很有用。但是我是否正确地理解了它不适合有多个输出变量的系统?如果是这样,有人知道有一个工具可以对有多个输出变量的系统进行多项式拟合吗?

约翰D 'Errico

Eran -这对一些人来说可能很有趣,但我认为这会鼓励人们去适应那些要求太高的模型,这是我经常被迫提醒人们不要做的事情。

因为您的想法需要在代码中添加更多的选项,这将使代码更加复杂。同样,使用该工具也不那么简单,因为如果生成并返回许多模型,那么像polyvaln和polyn2sym这样的工具就不能直接工作,或者也需要修改。

就效率的显著提高而言,请记住,当代码中有更多选项时,对于基本用户来说效率就会降低。复杂的代码也更有可能出现意想不到的错误。

所以我自己不会做这种改进。你自己能做到吗?是的,理论上可以。为了提高效率,需要围绕建模过程进行循环,可能使用qrupdate函数从模型中添加或删除术语。与在整个代码中重复调用polyfitn相比,我不确定这是否值得付出努力,除非您广泛使用该工具。

伊兰

你好再次,

非常感谢你上次的回答!
当然,这是一个愚蠢的问题:)

我还有一个问题,可能有点棘手
假设我有n个点的集合。
我想检查对于每一个k当然我可以做一个for循环得到我的答案…但这样效率很低。

有什么方法可以稍微修改一下您的代码以使其正常工作吗?

我认为主要问题在于你定义“设计矩阵”的部分。

再次感谢你的快速帮助!

约翰D 'Errico

作为延续,让我这样解释它,因为这是回归中常见的问题。

当您创建一个y = x/2的数据集时,您的数据位于一个不填充(x,y)平面区域的区域中。相反,它完全在直线y = x/2上。实际上,x和y完全依赖于彼此,给定x,我们知道y的值。

问题是,我们不知道当你的点脱离这条线时z会发生什么,所以没有一个完整模型的系数估计是令人满意的。

为了让你相信发生了什么,让我们从你创建的原始函数z开始。

z = 4 + x + x, y ^ 2 + x。*;

将已知的y与z的关系代入z,得到它只是关于x的函数。然后我们会得到:

Z = 4 + x + 1.5*x^2

我们可以把z建模为x的函数。

= polyfitn (x, z, 2)
了=
ModelTerms: [3x1 double]
系数:[1.5 1 4]
ParameterVar: [1.5362e-32 1.6407e-28 7.6808e-26]
ParameterStd: [1.2394 -16 1.2809e-14 2.7714e-13]
R2: 1
AdjustedR2: 1
RMSE: 9.327 e-13
VarNames: {X1的}

这里没有任何警告,而且拟合基本上是完美的,产生了我们所期望的系数值。我可以很容易地把z的拟合作为y的函数。

约翰D 'Errico

嗨,伊兰,

但是想想你在做什么。

乍一看,我迅速地估计x是中等大小的,因此这只是一个缩放问题。然后我看到了你的问题。就是这条线出了问题!

Y = x/2

你的模型是x和y的二次方程,也就是说它有x ^2 x *y和y ^2。当然,你希望polyfitn会告诉你x ^2的系数是0,就像你创建数据时隐式的那样。

但是你创建的y是x的简单倍数,你应该看到x *y = 2*y。^2 == x ^2/2。

世界上没有线性代数方案可以智能地恢复这些项的系数,并知道你从哪组系数开始。解决方案不是唯一的。这就是为什么你会得到关于奇点的警告。你会得到一个答案,但不是唯一的答案。

这里的故事就是垃圾进,垃圾出。您的数据不足以支持您希望拟合的模型的估计。金宝app

伊兰

嗨,乔恩,

非常感谢你做了这个非常需要的功能。可惜matlab不支持这个…金宝app

我尝试使用这个函数,但遇到了一个小问题。
我使用:
x = 0:100;
y = x / 2;
z = 4 + x + x ^ 2 + x.y;
在= (x, y);

= polyfitn (z, 2);

我得到的是一个关于由函数qr产生的矩阵R的奇异性的警告。当然,当我使用“\”时,实际上是在使用“R”的倒数。
在使用函数“qr”之后,这发生了两次。

你对如何解决这个问题有什么建议吗?

我问这个问题,当然是因为我知道“polyfitn”产生的系数与你所期望的非常不同……

非常感谢!

约翰D 'Errico

嗨,亚历克斯,

我觉得你想用ezsurf,而不是ezplot。例如,这对我来说很管用:

Xy = rand(100,2);
z = 3 + 2 * (xy, 2)总和+总和(xy。^ 2,2)+ 4 *刺激(xy, 2) + randn (100 1) / 10;
模型= polyfitn(xy,z,2);
Fun = polyn2sym(模型);
ezsurf(有趣)

ECOSat团队

嗨,约翰,

这些功能太棒了!

我用它来模拟两个垂直锥的交点,这为我提供了一个三维空间中的抛物线。我想知道,有没有一种内在的方法来绘制合适尺寸的曲线?我尝试使用你的polyn2sym工具,然后使用ezplot,但它不能正确地显示3D曲线。

谢谢,
亚历克斯
Ecosat团队

约翰D 'Errico

你好雷蒙德,

导出了这些参数,并将其应用于有噪声的模型。例如,假设我在一个点上对一个给定过程进行了10次测量。常数模型的最小二乘最佳估计,即y=c,由样本的平均值给出。同时,我们可以计算不确定性的估计值,用它来计算参数周围的近似置信区间。我们可以用更复杂的回归模型做类似的事情,根据数据中感知到的噪声,对模型中的参数和这些参数周围的不确定性进行估计。噪声被认为基本上是模型估计后剩下的东西。这样我们就能估计出参数本身的标准差。(将提供的标准偏差乘以学生的t统计量,以获得所需置信度水平的近似置信区间。使用大约+/- 2或3西格玛是这里常见的事情,以获得置信区间。)

事实上,大多数时候,当你有比参数更多的数据点来拟合一个模型时,我们也可以在这些参数中形成对不确定性的估计。当然,你拥有的数据越多越好,因为它可以让我们更好地估计参数。一件好事是我们甚至不需要纯粹的复制来形成参数中不确定性的估计。不确定性是基于模型无法在数据中预测的内容。

传统上,对于多项式模型,这些标准偏差估计是决定在模型中删除项的一种方法。如果置信区间包括零作为一种可能性,那么我们可能得出结论,我们无法将给定参数的估计值与零区分开来。

至少在传统模型中是这样的,我们试图将正确的模型适合我们的数据,而所有阻碍这一过程的是随机(高斯)噪声的存在。上述近似也假定参数是独立的,因此根据数据的不同,x和x^2通常不是正交变量。正交性的缺乏导致回归估计中的相关性,不幸的是,在用于报告的参数标准差的近似值中忽略了这一点。

当模型缺乏拟合时,这种方法还存在其他问题。在这里,我们试图让一个模型适合数据,但实际上,这个过程并不是真正基于这个模型的。我们选择的模型可能只是接近现实,也可能完全是错误的。假设你有一个曲面,你试图估计一个线性模型。在这种情况下,会缺乏拟合,因为模型没有很好地表示数据。

当然,线性代数不理解缺乏拟合不只是纯粹的噪声。数学在这方面有点笨。它不能像我们一样看到残差的模式。它简单地认为残差是纯噪声,并使用标准统计技术使用该噪声来估计参数中的不确定性。

最后,我想说的是,这些标准偏差粗略地衡量了软件认为参数估计的好坏,基于你提供的数据,基于几个近似值,基于你的模型对这些数据的有效性。

雷蒙德曾

嗨,约翰,

这是一个很棒的工具!我用它来拟合三维晶格上磁场的二次函数。

我感兴趣的是拟合参数的不确定性。将“ParameterStd”解释为这种不确定性正确吗?(如果不是,那是什么?)鉴于我没有向代码提供测量不确定度,我对如何确定“ParameterStd”感到有点困惑。你能解释一下这个“ParameterStd”是如何计算的吗?

谢谢你!
雷蒙德

约翰D 'Errico

对系数的约束不是选项之一。相反,您可以使用lsqnonneg之类的工具来对系数进行基本的非负约束。例如,如果x和y是列向量,那么为了得到a和b作为向量的非负系数,你会这样做:

Ab = lsqnonneg([x, x.^2],y);

对于系数上更复杂的约束,通常需要优化工具箱中的lsqlin之类的工具。例如,如果您只想约束其中一个系数,或者希望约束系数的总和,那么lsqlin将很有用。如果没有这个工具箱,还有其他工具,比如我自己的来自文件交换的lee,可以按照这些方式为您做一些事情,并且通过一些努力,您总是可以将许多一般问题转换为lsqnonneg可以处理的问题。

yonatan

太棒了!
有没有办法限制函数的系数为正?
例如,f(x) = a*x + b*x^2对于正a和b?
提前表示感谢。

安德鲁

这太棒了!

Syeda

陈云锅

约翰D 'Errico

Martin -你需要列出模型中的术语。例如,要指定一个总阶数不高于3,但y中没有三次项的模型,您可以将modelterms参数指定为

Modelterms = [3 0 .
2 1;
2 0;
1 - 2;
1 1;
1 0;
0 2;
0 1;
0 0]

或者,你也可以为相同的模型提供一个字符串的术语列表:

modelterms = ' x ^ 3, x ^ 2 * y, x ^ 2, x * y ^ 2, x * y, x, y ^ 2, y,不变的

马丁

嗨,约翰,首先,干得好!
我的问题是,多项式有可能有两种阶数吗?因为我设法用阶6 (n=6)的多项式来拟合我的函数(3D)。但如果用一个多项式来拟合我的函数就更好了,比如在一个方向上是5,在另一个方向上是6 (m和n)。
z (x, y) = (n) * x ^ n + b (m) * y ^ m + (n - 1) * x ^ (n - 1) + b (m - 1) * y ^ (m - 1)…

约翰D 'Errico

Jeff - polyfitn返回结构体的原因是因为模型可能有些复杂,用户在几个变量中指定了通用术语。1-d多项式更容易描述为系数列表,因为1-d多拟合总是假设一个非常特定的模型形式。

然而,从返回的结构中提取系数完全不是问题。例如,考虑这个简单的模型z(x,y)。

一些简单的数据…
X = rand(10,1);
Y = rand(10,1);
Z = x + 2*y + 3 + randn(size(x)))/10;

生成一个模型
MDL = polyfitn([x,y],z,{'x','y','constant'})
mdl =
ModelTerms: [3x2 double]
系数:[0.9706 2.0162 3.0302]
ParameterVar: [0.0214 0.0266 0.0093]
ParameterStd: [0.1462 0.1631 0.0966]
R2: 0.9787
AdjustedR2: 0.9726
RMSE: 0.1180
VarNames: {'x' 'y'}

当然,系数只是结构域中的数字。所以我们可以把系数提取成一个向量:

格式长
C = mdl。系数
C =
0.970599545975740 2.016179608143575 3.030225230230983

它们并不完全是我输入的,但数据上的随机噪声有点损坏了东西。

或者我们甚至可以把结构体转换成一个符号多项式,或者是一个符号,或者是我自己的符号形式。

polyn2sympoly (mdl)
ans =
0.97059954597574*x + 2.01617960814358*y + 3.03022523023098

我们可以自己求出这个多项式而不需要多项式,只要你知道模型是什么。(因为每个项的指数都存储在结构体中,这也很容易。)

C(1) *。5 + c(2)*。25 + c (3)
ans =
4.019569905254746

比较:
polyvaln (mdl [5, .25])
ans =
4.019569905254746

杰夫

约翰,干得好。问题:我想我读到过用polyfitn解出的系数只能用polyvaln解释。对吗?如果我想使用polyfitn来寻找数据拟合练习的系数,然后将它们应用到matlab之外的未来数据集(因此没有polyvaln),这是不可能的吗?还是我误解了什么?

埃里克·H。

丹尼斯·伍特斯
谢谢你提供这个解决方案。它的工作原理!

丹尼斯·伍特斯

约翰,干得好。
我遇到了一个小故障。
当将modelterms指定为指数数组时,变量'varlist'没有赋值。
我通过在第160行插入以下代码来修复这个问题:
Else varlist = repmat({"},1,p);

约翰D 'Errico

尽管我写过polyfitn,但高阶多项式(尤其是多维的)很少能很好地替代具有适当行为的非线性模型。也就是说,如果你有这样一个模型。这些多项式往往会对你产生不好的影响,尤其是在推断中。除非您有足够的数据,否则即使在您的数据中,模型也可能表现得很奇怪。我称之为插值问题。

基本上,你需要非常小心高阶多项式拟合。

顺便说一下,nlinfit确实可以处理多个自变量。你说没有,但这是没有限制的。问题是在多个维度上选择一个智能模型。一般来说,这是非常困难的。一个很好的原因是,像w = f(x,y,z)这样的函数形式很难形象化。

就我个人而言,我总是倾向于用样条模型来解决这些问题。我经常无法凭空得出一个模型,我更倾向于避免高阶多项式。我自己的网格拟合就是一个简单的例子,有两个自变量。

卡兰

你好约翰,

谢谢你的意见。非常感谢你的帮助。从我的理解是,数据是不充分的,不足以使用这样一个高阶多项式拟合。实际上,我正在考虑将数据拟合到我自己的自定义方程中,但还没有找到它的模型-类似于nlin拟合,但有多个自变量。如果你有任何建议,请告诉我。

谢谢你!

卡兰

约翰D 'Errico

嗨,卡兰,

在回归建模中,您陷入了一个常见的谬误。它是这样的:“如果n是好的,那么n+1更好,n+2一定非常棒。”

这种谬论很容易产生。我们在数据上尝试线性模型。它确实有效,但不是很好。毕竟,must问题只是局部线性的。所以我们尝试二次模型。它更合身。减少了不匹配。但这个模型仍然存在问题。它不足以符合我们数据的形状。

既然我们用二次函数做得很好,试试三阶模型。是的,它做得更好。哇。继续推进,直到我们得到一个六阶模型,或者十阶模型。多项式模型就像泰勒级数近似一样,所以项越多,结果就越好。但最终,这一切都变成了我们眼前的数字垃圾。到底发生了什么?

这里有3个变量。一阶模型有一个常数项,以及x, y和z项。这个模型中有4项。你有74个数据点,估计相应的系数肯定没有问题。

二次模型增加了6项,x^2, y^2, z^2,还有叉乘项x*y, x*z和y*z,所以总共有10个系数需要估计。

谬误的诅咒在于,当我们按顺序往上移动时,项的数量增长得越来越快。在3阶时,我们有20项要估计,在4阶时有35项,在5阶时有56项,在6阶时有84项。十阶模型有286项。

你的数据呢?你只有74个数据点。(在我早期做统计咨询的时候,我的客户总是会问他们需要多少数据。我的回答几乎总是“比你拥有的更多。”)

每个不同的数据点最多是一个额外的信息。(复制点在这里不计算。它们所做的就是减少你估算的不确定性。因此,复制点很好,但在某些方面不如独特点有价值。)

更糟糕的是,新的数据点并不总是提供新的信息。例如,两点决定一条直线。该线上的第三个点对该线的斜率或截距没有增加新的信息。就像复制一样,它所做的就是减少你估计系数时的不确定性。

那么,当你试图估计一个只有74个数据点的六阶模型时会发生什么呢?您的数据中信息太少,无法支持如此大的模型。金宝app返回的消息来自MATLAB,告诉你你的矩阵是数值奇异的(秩亏的)。

说了这么多,情况可能会变得更糟,使问题更加困难。有些数据比其他数据更糟糕。例如,如果你把1e6加到你所有的自变量上,一个给定顺序的多项式在数学上仍然是可估计的(理论上),然而它可能会得到一个数值上奇异的矩阵。这是因为线性代数是使用双精度浮点算法编写的。这些数字的动态范围可能太大,以至于代数运算无法成功。(当你把1e6数量级的数依次提高到更高次幂时会发生什么?)

所以当你有更多的系数估计比数据点,你将总是有一个秩亏矩阵,但取决于数据本身,你可能仍然有秩亏问题。

如果您将变量居中并缩放,一些多项式模型将变得更易于估计。因此,使每个变量的均值为零,并且大致位于区间[-1,1]。请注意,将正/-1数量级的数字提高到更高的幂并不会产生动态范围问题。这将迫使您转换您的问题,但它通常允许您获得对其他棘手问题的可行估计。(使用正交多项式是另一种经常有用的方法。)

我希望这能有所帮助。说真的,这是我看到的一个学期的课程,所以我只接触到表面。

约翰

卡兰

非常感谢。我得到这个警告-
警告:Rank不足,Rank = 64, tol = 5.8216e-010。
我使用了3个自变量,74个数据点,使用度为6。你能告诉我这是什么意思吗?谢谢你!

约翰·马奥尼

约翰,

和往常一样,代码很棒!

当只给出一个数据点时,它似乎失败了。

N = 1;
X = rand(n,2);
Y = exp(sum(x,2)) + randn(n,1)/100;
P = polyfitn(x,y,3);

虽然在某些情况下失败可能是可取的,但它实际上应该抛出一个错误。您的代码确实为n-1个输入提供了n个参数的输出,因此处理这种情况不会不一致。

我认为它在第120行之后的逻辑中-检查输入的大小并相应地转置。

谢谢!

斯科特DeWolf

约翰D 'Errico

这里有两个可能的问题。首先,如果你对所谓的多重右手边有问题,那么是的,它是非常可向量化的。我的意思是,问题的自变量集相同,但因变量集不同。然后反斜杠是一种方法,允许你分解方程一次,节省大量的时间。当然,你需要建立设计矩阵,但这并不难做。

然而,如果你有多个问题,每个问题都有不同的自变量集,那么向量化将要求你构建一个稀疏块对角设计矩阵。因此,您可以在这种努力中获得一些好处,但没有第一种情况的好处那么大。您需要将设计矩阵构建为稀疏矩阵。这里的技巧是将块作为单独的稀疏矩阵提供给blkdiag。因此,在单元格数组中创建单个(稀疏)设计矩阵,然后将其作为逗号分隔的列表传递给blkdiag。你能有效地做到这一点,使反斜杠求解比循环调用更有效吗?是的,但就像我说的,这需要一些努力。

您可能会从我的回答中了解到,对这段代码的结果进行向量化不是一件简单的事情,但是如果投入一些时间,您可以更有效地完成工作。

约翰

很棒的函数集!

你能建议一些方法来量化它们的使用吗?我有一个标量函数它依赖于它在空间f(x,y,z)中的位置,我计算这个函数围绕N个物体。目前,我为每个主体运行一个循环,但是否有一种方法可以修改代码,使反斜杠操作符一次性跨所有它们工作,返回N组系数?

约翰D 'Errico

这个应该可以了…

看看下面这两行:

[Q,R,E] = qr(M,0);

polymodel.Coefficients(E) = R\(Q'*depvar);

假设权重的列向量为W,每个数据点对应一个。把上面的行改为:

[Q,R,E] = qr(bsxfun(@times,W(:),M),0);

polymodel.Coefficients (E) = R \ (Q * (W(:)。* depvar));

这将为您提供参数的加权估计。我认为参数的统计也会起作用,还有R^2等等。当然,你也需要传入权重向量作为参数。

上述更改假设您使用的版本具有可用的bsxfun。它已经存在一段时间了,所以应该不是问题。即使没有,调用repmat也会替换它。

伟大的功能!谢谢你!有没有简单的方法把它转换成加权最小二乘。我想尝试自己做这件事,也许你有一些见解或警告?

约翰D 'Errico

线性回归(polyfitn)的核心使用反斜杠。例如,假设有人希望使用polyfit来估计一个变量中的二次模型。这个模型是

Y = a1*x^2 + a2*x + a3(+噪声)

Polyfit能够通过最小化残差平方和的方式求解线性(但超定)方程组来估计三个未知系数[a1,a2,a3]。同样,反斜杠是下面的引擎。

现在,假设有人希望估计一个更复杂的模型,有两个自变量z(x,y)。我把它写成一个完整的双变量二次模型,所以有6个未知系数。

Z = a1*x^2 + a2*x*y + a3*y^2 + a4*x + a5*y + a6 +噪声

只要你有足够的数据来估计六个未知系数,问题本质上是一样的。它们线性地进入问题,就像我们之前看到的单变量问题一样。反斜杠可以解决系数估计的超定问题。

事实上,polyfit和polyfitn都使用设计矩阵的QR分解来实现它们的目标。这是因为它们都被设计为对参数产生协方差矩阵估计。QR因式分解是同时满足这两个目的的最简单方法,但它只是线性代数。

D先生

谢谢,你的函数解决了我的问题,但说实话,我对二维曲面的线性回归概念很好奇。我们必须对每一行进行多拟合(一阶)来得到曲面的最佳拟合平面吗?

谢谢。

MATLAB版本兼容性
使用R2007b创建
与任何版本兼容
平台的兼容性
窗户 macOS Linux

社区寻宝

在MATLAB Central中找到宝藏,并发现社区如何帮助您!

开始狩猎!

PolyfitnTools /演示/ html /