function [a b]=fit_HP(x,y)
[d,N]=size(x);
B=x*x';
a=sum(x');
BB=[B a'];
A=[BB;[a N]];
XX=[x;ones(1,N)];
AA=XX*XX';
c=[x*y';sum(y)];
g=inv(A)*c;
a=g(1:2);b=g(3);
err=sum((y-a'*x-b).^2);
a
b
fprintf(' mean square error of training = %f\n',err/length(y));
return

          
      
        