发布网友 发布时间:2022-04-28 17:55
共1个回答
热心网友 时间:2022-06-22 17:56
使用最小二乘法拟合比较简单:
x_r=[abscissa ones(size(abscissa))]\ordinates;
求出来即为题中的x和γ。如果不限方法,也可以使用多项式拟合:
p = polyfit(abscissa, ordinates,1);
得到的结果是一致的(但二者分别是列向量和行向量)。
使用绝对值最小的拟合方法稍微复杂一些:
e = ones(size(abscissa));
f = [0; 0; e];
A1 = [-abscissa -e -eye(length(abscissa))];
b1 = -ordinates;
A2 = [abscissa e -eye(length(abscissa))];
b2 = ordinates;
c = linprog(f,[A1;A2],[b1;b2]);
求出来的即为所需的系数。
完整的代码如下(含绘图):
% Linear fit with least square and least absolute error
slope=3; intercept=-2;
abscissa = (-5:5)'; m = length(abscissa);
WhiteNoise = 5*randn(m,1);
ordinates = slope*abscissa + intercept + WhiteNoise;
GrossError=80;
ordinates(6)=ordinates(6)+GrossError;
ordinates(10)=ordinates(10)-GrossError;
x_r=[abscissa ones(size(abscissa))]\ordinates;
y2=[abscissa ones(size(abscissa))]*x_r;
e = ones(size(abscissa));
f = [0; 0; e];
A1 = [-abscissa -e -eye(length(abscissa))];
b1 = -ordinates;
A2 = [abscissa e -eye(length(abscissa))];
b2 = ordinates;
c = linprog(f,[A1;A2],[b1;b2]);
y1=[abscissa ones(size(abscissa))]*c(1:2);
plot(abscissa, ordinates, 'ro',abscissa, y2,'b-',abscissa, y1,'g--')
legend('Original data', 'Least L^2 error', 'Least L^1 error')
某次运行的结果如下(因数据为随机生成,每次运行结果不同):
对于L1最小拟合,也可以使用优化工具箱的函数fminunc来做:
objfun = @(x)sum( abs(abscissa*x(1)+x(2)-ordinates) );
c2=fminunc(objfun,[2 -2])
y3=[abscissa ones(size(abscissa))]*c2.';
plot(abscissa, ordinates, 'ro',abscissa, y2,'b-',abscissa, y1,'g--',abscissa, y3,'m:')
legend('Original data', 'Least L^2 error', 'Least L^1 error', 'fminunc')
拟合的结果可能与前述线性规划的结果不同(但目标函数基本都能达到最小),这也说明这个问题本质上是多极值的,存在多个局部最优点。