Least Trimmed Squares fitting

Is there any toolbox or written code for least trimmed square in Matlab?

Answer 1

[Matt J]

If you don't have too many unknown parameters to fit, you could implement it easily with fminsearch.

[NA]

How can I use the number of the trimming point in this function?

I have this data set.

x = (1:10)';
y = 10 - 2*x + randn(10,1);
y(10) = 6; y(6) = 15; y(4) = -10; y(1) = -10;
plot(x,y,'or')

Here, it is possible to remove 4 points. I do not know how to do this part.

About fminsearch, I use this

M = [ones(size(x)),x]; 
p0 = M\y;          
errfcn = @(p,y,M) sum((y-M*p).^2); 
p1 = fminsearch(@(p) errfcn(p,y,M),p0); 
hold on 
plot(x,y-M*p1,'.-')

I want to find the regression line, that removes this data y(10) = 6, y(6) = 15, y(4) = -10, y(1) = -10.

[Matt J]

x = (1:10)';
y = 10 - 2*x + randn(10,1);
y(10) = 6; y(6) = 15; y(4) = -10; y(1) = -10;

%%Initial guess
in=abs(y-movmedian(y,3))<=2;
p0=polyfit(x(in), y(in),1);

%%Optimize
k=4; %number to trim
p=fminsearch(@(p)trimlsq(p,x,y,k), p0);

%%Plot
xs=linspace(min(x),max(x),100);
plot(x,y,'o',xs, polyval(p,xs),'--')

function fval=trimlsq(p,x,y,k)

  r=abs(polyval(p,x)-y);
  rtrim=maxk(r,k);
  fval=norm(r)^2-norm(rtrim)^2;

end

[NA]

Thank you as always. For initial guess, I don't understand why you choose this way

%%Initial guess
in=abs(y-movmedian(y,3))<=2;
p0=polyfit(x(in), y(in),1);

[NA]

If x is matrix

x = [-122.9115,0.0000,-0.0000;-70.9142,0.0000,-0.0000;122.8232,-0.0000,0.0000;...
     0.0000,109.3731,-112.7629;0.0000,63.1032,-65.0589;0.0000,-116.0692,119.4565;...
     0,111.5914,0;0,0,109.0116]
y = [0.1415;0.1756;-0.9547;-4.5564;-0.6909;3.8466;5.6495;8.4522];

For initial guess

y = sparse(y);
x = sparse(x);
in = abs(y-movmedian(sparse(y),3))<=2;
p0 = polyfit(x(in), y(in),size(x,2)-1); % estimate 3 unknown variable
% or
% p0 = x\y; 
p = fminsearch(@(p)trimlsq(p,x,y,1), p0);

In this way, I could not get good regression

[Matt J]

Why would you choose k=1?

[NA]

Why would you choose k=1?

Removing more than 1 row from the x make the matrix rank deficient.

[Matt J]

I don't understand the meaning of x being a matrix when y is not.