function p = mejor(x, y)

n = length(x(:,1))
k = length(x(1,:)) + 1

fcal = finv(0.95, 1, n-2)
fmax = fcal;

mod = [];
for l = 1:k-1
   for m = 1:k-1
      m 
      mod
      if(isempty(mod) | prod(m ~= mod) ==1) 
   	   [SSE, SSR, p] = LineaRecta(x, y, [mod, m]);
         f = SSR/(SSE/(n-2))
         if(f > fmax)
            fmax = f;
            kmax = m
         end;
      end;
    end;
    mod = [mod, kmax]
end;

function [SSE, SSR, p] = LineaRecta(x, y, mod)

K = length(x(1, :));

M = ones(length(x(:,1)),1);

for m = 1: K
    if(m == mod) 
       M = [M, x(:,m)];
    end;
end;


p = inv(M'*M)*M'*y;

SSE = (M*p-y)'*(M*p-y);

yp = mean(y);

SSR = (M*p-yp)'*(M*p-yp);