least absolute deviation when we have data set
8 views (last 30 days)
Show older comments
I have this data
x = (1:10)';
y = 10 - 2*x + randn(10,1);
y(10) = 0;
how can I use least absolute value regression?
0 Comments
Accepted Answer
Bjorn Gustavsson
on 22 Aug 2020
You can do it rather straight-forwardly with fminsearch (or other similar tools on the file exchange: fminsearchbnd, minimize etc):
M = [ones(size(x)),x]; % Matrix for linear LSQ-regression, we could do centering and scaling etc...
p0 = M\y; % straight least-square fit - to get ourselves a sensible start-guess (hopefully)
errfcn = @(p,y,M) sum(abs(y-M*p)); % L1 error-function
p1 = fminsearch(@(p) errfcn(p,y,M),p0); % L1-optimization
subplot(2,1,1)
plot(x,y,'.-')
hold on
plot(x,M*p0)
plot(x,M*p1)
subplot(2,1,2)
plot(x,y*0,'.-')
hold on
plot(x,y-M*p0,'.-')
plot(x,y-M*p1,'.-')
% For my test I got L1-error-function-value for the least-square-fit p0:
% errfcn(p0,y,M)
% ans =
% 22.058
% and for the L1-optimal parameters:
% >> errfcn(p1,y,M)
% ans =
% 20.067
This would generalize to more interesting problems too. Also have a look at Huber-norms, for an error-norm kind of intermediate between L1 and L2.
HTH
8 Comments
Bjorn Gustavsson
on 7 Sep 2020
Because they use different algorithms, and from the robust-fit documentation you can look up the weighting used for its different settings of wfun and tune. Do the regressions differ by much? How do they vary if you vary the different tuning-parameters? When using robust fitting you should always check the residuals and their relative contributions to the total error-function.
More Answers (1)
Bruno Luong
on 22 Aug 2020
Edited: Bruno Luong
on 22 Aug 2020
% Test data
x = (1:10)';
y = 10 - 2*x + randn(10,1);
y(10) = 0;
order = 1; % polynomial order
M = x(:).^(0:order);
m = size(M,2);
n = length(x);
Aeq = [M, speye(n,n), -speye(n,n)];
beq = y(:);
c = [zeros(1,m) ones(1,2*n)]';
%
LB = [-inf(1,m) zeros(1,2*n)]';
% no upper bounds at all.
UB = [];
sol = linprog(c, [], [], Aeq, beq, LB, UB);
Pest = sol(m:-1:1); % here is the polynomial
% Check
clf(figure(1));
plot(x, y, 'or', x, polyval(Pest,x), 'b');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/349955/image.png)
3 Comments
Bruno Luong
on 22 Dec 2020
Edited: Bruno Luong
on 22 Dec 2020
"2) I am totally new to the ways of linear programming, so I am wondering how come you have no inequality constraints? I am guessing you are saying the solution must adhere to the objective function, precisely? "
Because I don't need it. I formulate the problem as
M*P - u + v = y
where u and v a extra variables, they meant to be positive
v =( M*P - y) = u
so
argmin (u + v) is sum(abs( M*P - y)) is L1 norm of the fit.
I could formulate with inequality but they are equivalent. There is no unique way to formulate LP, as long as it does what we want.
And as comment; all LP can be showed to be equivalent to a "canonical form" where all the inequalities are replaced by only linear equalities + positive bounds
argmin f'*x
A*x = b
x >= 0
See Also
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!