ATLAB 脚本进行最小二乘法拟合

[复制链接]

字体大小: 正常放大

2744557306

1188 主题	4 听众	2931 积分

该用户从未签到

电梯直达

1^#

发表于 2023-12-31 16:12 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

这是一个 MATLAB 脚本，用于进行最小二乘法拟合。脚本首先要求用户输入已知点的 x 和 y 坐标，然后输入拟合的多项式次数 n。脚本使用最小二乘法拟合数据，并绘制了原始数据点和拟合曲线的图表。以下是对代码的主要部分的解释：
function fp = fitpt()
% 最小二乘
% 基取 {1, x, ...}
% fitpt.m

% 默认算例为课本：P65，例3.2
% x = [0,1,2,3,4,5,6,7]
% y = [3.95,6.82,9.78,12.91,15.74,19.26,21.73,24.07]
% 结果：P(x) = 4.005 + 2.936x  平方误差=0.6162

% MatLab函数：polyfit(x, y, n)

s = input('<最小二乘>\n输入已知点的x坐标：（回车表示[0,1,2,3,4,5,6,7]）\n', 's');
if isempty(s)
      s = '[0,1,2,3,4,5,6,7]';
else
      if (s(1) ~= '[')
         s = strcat('[', s);
         s = strcat(s, ']');
      end
end
x = sym(s);

s = input('输入已知点的y坐标：（回车表示[3.95,6.82,9.78,12.91,15.74,19.26,21.73,24.07]）\n', 's');
if isempty(s)
      s = '[3.95,6.82,9.78,12.91,15.74,19.26,21.73,24.07]';
else
      if (s(1) ~= '[')
         s = strcat('[', s);
         s = strcat(s, ']');
      end
end
y = sym(s);
sz = size(x);
sz = sz(2);
n = input('输入多项式次数n：');
if (n + 1 > sz)
      n = input('多项式次数需要小于已知点个数，请重新输入n：');
end
if (n + 1 > sz)
      error('多项式次数不能小于已知点个数！');
end
fp = s_fitpt_p(x, y, n);

% 绘制原始数据点和拟合曲线
plot(double(x), double(y), 'r*')
hold on
a = double(x(1));
b = double(x(sz));
x = a:abs(b - 1)/100:b;
y = subs(fp, x);
plot(x, y)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = s_fitpt_p(x, y, n)
% 用 n 次多项式实现的最小二乘法

sz = size(x);
sz = sz(2);
A = zeros(sz, n + 1);
v = vh(n);
for i = 1:sz
      A(i,

= subs(v, double(x(i)));
end
f = linsolve(A' * A, A' * y');
f = vpa(f, 4);
f = v * f;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function v = vh(n)
% Create vector in horizontal style, such as
% v = [1, x, x^2, ..., x^n]

if (n < 0 || n > 9)
      error('Make sure ''n'' is in range of [0, 9]')
end
s = '';
for i = 0:n
      s = strcat(s, ',x^');
      s = strcat(s, num2str(i));
end
s(1) = '[';
sz = size(s);
s(sz(2) + 1) = ']';

v = simplify(sym(s));
end

这个脚本首先获取用户输入的已知点的 x 和 y 坐标，然后使用最小二乘法进行拟合。最后，脚本绘制了原始数据点和拟合曲线的图表。

zan