least absolute deviation matlab

Hi,

Is the a way to do parameter estimation using Least Absolute deviation with constraints in matlab.

If so, please let me know of a place to read about it.

thanks.

Answer 1

[Matt J]

If you are talking about the problem

 min_x  sum( abs(r(x)))

it is equivalent to

min_{x,d} sum(d)
s.t.
     r(x) <=d
    -r(x)<=d

The above formulation is differentiable if r(x) is differentiable, so FMINCON can handle it. You can obviously add any additional constraints you wish, so long as they too are differentiable.

[dav]

Hi,

Thanks for your answer.

I have a dataset yt (all positive - size (100,1)). I need to fit a straight line to this data with the constraints that intercept parameter is greater than zero and slope parameter is between 0 1nd 1.

Could you please help me with this code? It is a little difficult to figure it out using your answer.

Thanks much!

[Matt J]

So, the problem is as follows?

 min.

    f(m,b) = sum( abs(yi-m*xi-b))

 s.t. 
      b>=0
      0<=m<=1

[dav]

yes, but

b > 0 
0<= m < 1

Thanks.

[Matt J]

You cannot impose strict inequality constraints. The minimum can be ill-defined, for example in the minimization of f(a)=a^2 over a>0. If a=0 is not allowed, then what do you consider the minimizing point?

[dav]

Thanks. can I have something like a >0.00000001. Would you please help me with the code ?

[dav]

I can provide the code which generates the data. but I dont know how to estimate it?

[dav]

   The data set is yt. 
  clc;
  clear;

T = 300;

    a0 = 0.1; a1 = 0.4; 

    ra = zeros(T+2000,1);
    seed=123;
    rng(seed);

       ra = trnd(5,T+2000,1); 

         ytn = [];

    epsi=zeros(T+2000,1);
    simsig=zeros(T+2000,1);

    unvar = a0/(1-a1);

    for i = 1:T+2000
                     if (i==1) 
                         simsig(i) = unvar;
                         s=(simsig(i))^0.5;
                         epsi(i) = ra(i) * s;

                     else                     
                         simsig(i) = a0+ a1*(epsi(i-1))^2;
                         s=(simsig(i))^0.5;
                         epsi(i) = ra(i)* s;
                     end
    end

    epsi2 = epsi.^2;
    yt = epsi2(2001:T+2000);

[Matt J]

I don't understand your code, but here's an example for you to study

   %%Fake data

     xi=0:10;
     yi=2*xi+bt+randn(size(xi))/2;

    %%Linprog parameters 

     N=length(xi);
     e=ones(N,1);

     f=[0,0,e.'];
     A=[xi(:),e,-speye(N);-xi(:), -e, -speye(N)];
     b=[yi(:);-yi(:)];

     lb=zeros(N+2,1);
     ub=inf(N+2,1); ub(1)=1;

     %%Perform fit and test

     p=linprog(f,A,b,[],[],lb,ub);

     slope=p(1),
     intercept=p(2),

     yf=slope*xi+intercept;
     plot(xi,yi,xi,yf)

[dav]

Thank you very much

Could you please explain the following codes

f=[0,0,e.'];
A=[xi(:),e,-speye(N);-xi(:), -e, -speye(N)];
b=[yi(:);-yi(:)];

and also bt in yi=2*xi+bt+randn(size(xi))/2;

thanks

[Matt J]

Recall that the reformulation of the problem requires additional constraints

     r(x) <=d
    -r(x)<=d

where r(x) is your residual. The A,b pair that you've cited is the implementation of these additional constraints. bt is the true intercept of the simulated line data, yi. You can choose any value for it that you want, for simulation purposes, and see if the p(2) returned by linprog ends up being a good estimate of it.

[dav]

is it possible to add the constraint that the parameter estimates should be positive here?

[Matt J]

Not sure what you mean by "add". I already put positivity constraints in my example using the lb input argument to linprog. Similarly, I used ub to impose the upper bounds you originally mentioned.

[dav]

thank you.

I used your code to test the data set I have. BOTH parameter estimates should be 0.1. However, the estimates I get using your code are very different. is there a way to fix it. From a standard theory I know that when you regress the y vector I have mentioned on the x vector I should get parameters, both equal to 0.1

clc;
clear;

p=1;
T = 300;

a0 = 0.1; a1 = 0.1; 

seed=123;

ra = randn(T+2000,1); 

epsi=zeros(T+2000,1);
simsig=zeros(T+2000,1);

unvar = a0/(1-a1);

for i = 1:T+2000
                 if (i==1) 
                     simsig(i) = unvar;
                     s=(simsig(i))^0.5;
                     epsi(i) = ra(i) * s;

                   else                     
                       simsig(i) = a0+ a1*(epsi(i-1))^2;
                       s=(simsig(i))^0.5;
                       epsi(i) = ra(i)* s;
                   end
  end

epsi2 = epsi.^2;

y = epsi2(2001:T+2000); % THIS IS THE DATA SET I WANT TO TEST. 

len = length(y);

x = zeros(len,p);

for i = 1:p
       x(1+i:len,i) = y(1:len-i,1);  % THIS IS THE x VECTOR
end

       N=length(x);
       e=ones(N,1);

       f=[0,0,e.'];
       A=[x(:),e,-speye(N);-x(:), -e, -speye(N)];
       b=[y(:);-y(:)];

       lb=zeros(N+2,1);
       ub=inf(N+2,1); ub(1)=1;

       %%Perform fit and test

       p=linprog(f,A,b,[],[],lb,ub);

       slope=p(1),
       intercept=p(2),

[Matt J]

All I can say is that I tested the example code I gave you. It worked fine for me and produced accurate estimates.