%--------------------------------------------------------------------------
% Poly_Reg is a function developped by J.C. COMTE comtejc@gmail.com
% CNC UMR CNRS 5229.
% Last Release July. 2007
% This work is protected by the CeCILL-C Licence (see below).
%
%
% This function performs nonlinear or polynomial regression.
% x, and y are the two input vectors, ordr is the maximum degrees
% of the polynomial fit.
%
% The output vector P, contains the polynomial coefficents from
% higher to the lower power. Ex: P(1)*X.^3+P(2)*X.^2+P(3)*X+P(4).
%
%
% This software is governed by the [CeCILL|CeCILL-B|CeCILL-C] license under French law and
% abiding by the rules of distribution of free software. You can use,
% modify and/ or redistribute the software under the terms of the [CeCILL|CeCILL-B|CeCILL-C]
% license as circulated by CEA, CNRS and INRIA at the following URL
% "http://www.cecill.info".
%
% As a counterpart to the access to the source code and rights to copy,
% modify and redistribute granted by the license, users are provided only
% with a limited warranty and the software's author, the holder of the
% economic rights, and the successive licensors have only limited
% liability.
%
% In this respect, the user's attention is drawn to the risks associated
% with loading, using, modifying and/or developing or reproducing the
% software by the user in light of its specific status of free software,
% that may mean that it is complicated to manipulate, and that also
% therefore means that it is reserved for developers and experienced
% professionals having in-depth computer knowledge. Users are therefore
% encouraged to load and test the software's suitability as regards their
% requirements in conditions enabling the security of their systems and/or
% data to be ensured and, more generally, to use and operate it in the
% same conditions as regards security.
%
% The fact that you are presently reading this means that you have had
% knowledge of the [CeCILL|CeCILL-B|CeCILL-C] license and that you accept its terms.
%
%--------------------------------------------------------------------------
function [P,x,S]=Poly_Reg(X,Y,ordr)
[lx,cx]=size(X);
[ly,cy]=size(Y);
if (cx>lx) X=X'; end;
if (cy>ly) Y=Y'; end;
N=2*ordr;
x(:,1)=X;
for i=1:N,
x(:,i)=X.^i;
end;
x=mean(x);
for i=1:ordr+1,
for j=1:ordr+1,
if ((i==ordr+1)&&(j==ordr+1))
xm(j,i)=1;
else
xm(j,i)=x(N+2-i-j);
end;
end;
end;
for i=0:ordr,
yxm(i+1,1)=mean(Y.*(X.^(ordr-i)));
end;
P=inv(xm)*yxm;
S=0;
m=min(X);M=max(X);
stp=abs(M-m)/10;
stp=stp/100;
x=m-stp:stp:M;
for i=0:ordr,
S=S+P(i+1)*x.^(ordr-i);
end;