%
% Reset to a known state.
%
clear
rand('state',0);
randn('state',0);
%
% Make X, Y, and SIGMA global.
%
global X;
global Y;
global SIGMA;
%
% Generate the data set.
%
X=(1.0:0.25:7.0)';
ptrue=[1.0; -0.5; 1.0; -0.75];
ytrue=rawfunc(X,ptrue);
SIGMA=0.01*ones(size(ytrue));
Y=ytrue+0.01*randn(size(ytrue));
%
% Output the data set.
%
disp('Data set for Example 9.2:')
pause
[X,Y]
%
% Now, go ahead and solve the problem with many different random starting 
% points.
%
disp('Now we solve the problem with many different random starting points:')
pause
N=20;
models=zeros(4,N);
chisqvals=zeros(N,1);
normg=zeros(N,1);
itercnts=zeros(N,1);
for i=1:N
  p0=2*(rand(4,1)-0.5*ones(4,1));
  [pest,itercnts(i)]=lm('fun','jac',p0,1.0e-14,500);
  models(:,i)=pest;
  chisqvals(i)=chi2(pest,X,Y,SIGMA);
  normg(i)=norm(jac(pest)'*fun(pest));
end
%
% With the randomly generated starting points used here, it happens
% that the best of the locally optimal solutions is found from the 
% second starting point, the next best is found from the fifth 
% starting point, and the third best solution is found from
% starting point one.
%
% Make a table of the locally optimal solutions. (2, 5, 1).  
%
disp('Generate a table of locally optimal solutions from 20 starting points:')
pause
table=[models(1,2) models(2,2) models(3,2) models(4,2) chisqvals(2);
       models(1,5) models(2,5) models(3,5) models(4,5) chisqvals(5);
       models(1,1) models(2,1) models(3,1) models(4,1) chisqvals(1)]
%
% Find the best parameters, and covariance matrix.
%
pbest=models(:,2)
Jbest=jac(pbest);
covbest=inv(Jbest'*Jbest);
corbest=zeros(4,4);
for i=1:4
  for j=1:4
    corbest(i,j)=covbest(i,j)/(sqrt(covbest(i,i)*covbest(j,j))); 
  end 
end
fprintf('The correlation matrix is \n');
corbest 
%
% Plot the fitted model and data points.
%
disp('The fitted model:')
pause
yfit=rawfunc(X,pbest);
figure(1);
clf;
bookfonts;
plot(X,yfit,'k');
hold on;
%
% one sigma error bars.
%
errorbar(X,Y,SIGMA,'ko');
xlabel('x');
ylabel('y');
%print -deps nlregfit.eps
figure(2);
clf;
bookfonts
plot(X,fun(pbest),'ko');
axis([0 7 -2 2]);
xlabel('x');
ylabel('normalized residual');
%print -deps residuals.eps



