matlab中的多项式拟合谷点

Question

我正在尝试找到最适合输入点集的多项式。

这是我目前的代码：

    x=(1:length(meanValues));
    y=meanValues(:);        

    A=fliplr(vander(x));
    v=A \ y;
    P(1: length(x))=0;
    for i=1: length(x)
        for j=1: length(v)
            P(i)=P(i)+v(j)*x(i).^(j-1);
        end
    end    
    plot(x,y,'r*');
    hold on;
    plot(x, P);

• meanValues 是 [1x127] 向量，填充了 (0.0000-5.0000] 之间的双精度值

下面是绘制的平均值：

结果：

有人知道错误在哪里吗？

编辑 1：

所以这一次，我遍历了所有多项式阶并找到了最合适的阶。这是否更好？我可以优化这段代码吗？计算大约需要 1 秒，因此总共需要大约 30 秒。

    tic
    x=(1:length(meanValues));
    y=meanValues(:)'; 
    for i=1:length(meanValues)-1
        [p,s,mu] = polyfit(x,y, i);
        [f,delta] = polyval(p,x,s,mu);
        if i==1
            minf=f;
            minmse = mean(delta.^2);
            minp=p;
        elseif minmse>mean(delta.^2)
            minf=f;
            minmse = mean(delta.^2);
            minp=p;
        end
    end
    toc
    plot(x,y,'r*',x,minf,'-');
    axis([0 length(meanValues) 0 max(meanValues)]);

Answer 1

您需要开发一些方法来迭代您想要考虑拟合的所有多项式阶数。在整个迭代过程中，您将计算模型 P 和数据 y 之间的误差。我建议，均方误差是衡量相似性的常用方法。

目前您无法更改模型顺序，实际上它非常高 (127)，因此您的最终结果不稳定。

在这个修改后的代码中，我生成了我自己的嘈杂平均值，这将是使用二阶拟合的最佳拟合。然而阶值设置为 4，所以你会发现 v 中的三阶和四阶系数与 0 阶、1 阶和 2 阶系数相比非常小。

至少对于我生成的数据，您应该能够验证二阶拟合的 y 和 P 之间的 MSE 低于四阶拟合。您的数据似乎没有太大的趋势，因此您最好测试几个不同的订单并选择 MSE 最低的订单。这并不是说它正确地模拟了生成数据的系统，所以要小心。

clear all;
meanValues = (1:127)/25;
meanValues(:) = meanValues(:).^2;
for i = 1:length(meanValues)
    meanValues(i) = meanValues(i) + rand(1,1)*4;
end

x=(1:length(meanValues));
y=meanValues(:);

Order = 4;
A(:,1) = ones(127,1);
for j = 1:Order
    A(:,j+1) = (x'.^j);
end
%   A=fliplr(vander(x));
v=A \ y;
P(1: length(x))=0;
for i=1: length(x)
    for j=1: length(v)
        P(i)=P(i)+v(j)*x(i).^(j-1);
    end
end
plot(x,y,'r*');
hold on;
plot(x, P);

编辑： 此版本计算 MSE 并找到最小订单。只需 0.324198 秒即可检查最多 100 次拟合。也许使用 polyfit 有一些优势……我不确定。

clear all;
meanValues = (1:127)/25;
meanValues(:) = meanValues(:).^2;
for i = 1:length(meanValues)
    meanValues(i) = meanValues(i) + rand(1,1)*4;
end

x=(1:length(meanValues));
y=meanValues(:);
tic
minMSE = Inf;
nOrder = 100;
for Order = 1:nOrder

    A(:,1) = ones(127,1);
    for j = 1:Order
        A(:,j+1) = (x'.^j);
    end
    %   A=fliplr(vander(x));
    v=A \ y;
    P = zeros(1,length(x));
    for i=1: length(x)
        for j=1: length(v)
            P(i)=P(i)+v(j)*x(i).^(j-1);
        end
    end
   P = P';
    newMSE = norm(P-y);
    if (newMSE < minMSE)
        minMSE = newMSE;
        minOrder = Order;
        minP = P';
    end
end

toc
plot(x,y,'r*');
hold on;
plot(x, minP);
minMSE
minOrder

Answer 2

这段代码可以正常工作：

% Data and regression
y = cumsum(randn(100, 1));
x=(1:length(y));
x = x(:);
A=fliplr(vander(x));
A = A(:, 1:7);
v=A \ y;

% Calculate P your way
P(length(x))=0;
for i=1: length(x)
  for j=1: length(v)
    P(i)=P(i)+v(j)*x(i).^(j-1);
  end
end

% Calculate P by vectorization
Q = A * v;

% P and Q should be the same - they are!
tmp = P - Q';
plot(tmp, '.')

% Plot data and fitted data
figure
plot(x,y,'r*');
hold on;
plot(x, P, '-b');
plot(x, Q, '-g');

这个回归等价于

p = fliplr(polyfit(x,y,6))';

返回警告

Warning: Polynomial is badly conditioned. Add points with distinct X
         values, reduce the degree of the polynomial, or try centering
         and scaling as described in HELP POLYFIT. 

如果你试试这个

A=fliplr(vander(x));
A = A(:, 1:8);
v=A \ y;

它返回一个警告：

Warning: Rank deficient, rank = 7,  tol =   5.9491e+000. 

因为A(end) 是1.0000e+014。

所以，你看，多项式回归是讨厌的方法。你必须另寻出路。