%% Generate the data
x = [    87    92    92    93    93    95    98    99    99   105   105  101  80    89    91    92    93    94    95    99    99   104   105   110    86   90    91    92    93    95    97    99    99   105   105    94]';
ndat = size(x,1);
y = ones(ndat,1) + 0.02* x -1.3*sin(x.^2/30) + exp(-0.2*x.^3/40);
y = y.*(1 + 0.2*rand(ndat,1));
%y = [   76    71    40    43    69    84    38    41    30    22    38    37    90    42    39    43    51    79    33    33    50    30    32    24    59  31    70    67    59    76    30    27    23    41    35    53]';
[x2,ix] = sort(x);
y2 = y(ix);
xs = [min(x):max(x)];

%% 
figure(1);clf;
hold on;
plot(x2,y2,'+');

%% now fit the data with models from polynomials
for nfit = 6:2:6
    
    J = zeros(ndat,nfit);
    for k = 1:nfit
        J(:,k) = x2.^(k-1);
    end
    afit = J\y2;
    yfit = afit(1);
    for j = 2:nfit
        yfit = yfit + afit(j) * xs.^(j-1);
    end
    figure(1); hold on;
    plot(x2,y2,'b+');
    plot(xs,yfit,'r-');
    pause(1);
end

%% "exact fit with N = M
nfit = ndat;
   J = zeros(ndat,nfit);

    for k = 1:nfit
        J(:,k) = x2.^(k-1);
    end
 %      alpha = 1e-12*trace(J); % regularisation !
       alpha = 0;

    afit = [J;alpha*speye(nfit)]\[y2;zeros(nfit,1)];
    yfit = afit(1);
    for j = 2:nfit
        yfit = yfit + afit(j) * xs.^(j-1);
    end
    
    % plot this on log scale because of instability
    figure(2);clf; 
    semilogy(xs,yfit,'r-');
    hold on;
    semilogy(x2,y2,'b+');