%   Program to demonstrate the use of the Levenberg-Marquardt algorithm
%   It works well if the initial estimates of the parameters are 
%   reasonably close to the actual parameter values but can fail if the
%   initial estimates are extremely bad

%   A Gaussian-form function gf has three parameters
%       A - Gain Factor
%       m - Mean value
%       s - Standard deviation

%   Parameter vector    x = [A,m,s]'
%   Observations        y = gf(x,z)

%   Close all open figure windows and open a new one

    close all ; figure('Position',[20,20,800,600]) ;
    
%   Initialize the true values of the parameters in the parameter
%   column vector x0

	A = 10 ; m = 30 ; s = 20 ; x0 = [A,m,s]' ; 

%   Create a column vector of z values at which to compute 
%   the observations y 

	z = [-20:1:100]' ;
    
%   Create the true y column vector y0

	y0 = gf(x0,z) ;
	
%   Add noise to the true y values to create the actual noisy
%   observation column vector y0n

	N = length(z) ; y0n = y0 + randn(N,1)/5 ;

%   Create a column vector xe0 of initial guesses at the parameter
%   column vector x

	Ae = 15 ; me = 20 ; se = 50 ; xe0 = [Ae,me,se]' ;
	
%   Use the Levenberg-Marquardt algorithm to find the 
%   parameter column vector xe which minimizes the mean-squared
%   error between the observation column vector y0n and the
%   computed observation column vector from xe

	xe = levmar('gf',xe0,y0n,z) ;

%   Graph the noisy observation column vector y0h, the best estimate 
%   observation column vector ye and the observation column vector
%   computed from the initial guesses at the parameter column vector xe0
%   
    ye = gf(xe,z) ; ye0 = gf(xe0,z) ;
    plot(z,y0n,'k.',z,ye,'k',z,ye0,'k--') ;
    xlabel('z','FontName','Times','FontSize',18) ; 
    ylabel('gf(z)','FontName','Times','FontSize',18') ;
    set(gca,'FontName','Times','FontSize',14) ;
    legend('Noisy observations','Best estimate','Initial guess')