function result = fit_rayleigh_pdf( x,y,W,hAx )
% fit_rayleigh_pdf - Non Linear Least Squares fit of the Rayleigh distribution.
%                    given the samples of the histogram of the samples, finds the 
%                    distribution parameter that fits the histogram samples.
%
%    fits data to the probability of the form: 
%        p(r)=r*exp(-r^2/(2*s))/s
%    with parameter: s
%
% format:   result = fit_rayleigh_pdf( x,y,W,hAx )
%
% input:    y   - vector, samples of the histogram to be fitted
%           x   - vector, position of the samples of the histogram (i.e. y = f(x,a))
%           W   - matrix or scalar, a square weighting matrix of the size NxN where
%                 N = length(y), or 0 to indicate no weighting is needed.
%           hAx - handle of an axis, on which the fitted distribution is plotted
%                 if h is given empty, a figure is created.
%
% output:   result  - structure with the fields
%                      s   - fitted parameter
%                      VAR - variance of the estimation
%                      type- weighted LS or not weighted LS
%                      iter- number of iteration for the solution
%

%
% Algorithm
% ===========
%
% We use the WLS algorithm to estimate the PDF from the samples.%
% The rayleigh distribution is given by:
%
%    p(x,s) = x*exp(-x^2/(2*s))/s
%
%    note that X is known and therefore, considered a constant vector
%
% The non liner WLS estimator is given by:
%
%    s(n+1) = s(n) + inv(H'*W*H)*(H') * (y-h) = s(n) + G * err
% 
%    where:   h = p(x,s)
%             H = diff( p(x,a) ) with respect to "s"
%             W = weighting matrix of size NxN  (N = length(y))
%             s = a single parameter to be estimated
%
% The error estimation is given by:
%
%    VAR( s ) = G * VAR( err ) * (G')
%
%       or when W=I and the noise is a gaussian noise 
%
%    VAR( s ) = inv( H' * H )
%

if (nargin<3)
    error( 'fit_rayleigh_pdf - insufficient input arguments' );
end

s       = x(find(y==max(y))).^2;        % initial guess
y       = y(:);                         % both should be column vectors !
x       = x(:);
x2      = x.^2;                         % save computation time
thresh  = 0.995;                        % convergence threshold for the loop
last_cnt= inf;
iter    = 0;

% check weight matrix input
if (size(W,1)==length(y)) & (size(W,2)==length(y))
    weights_flag    = 1;
    type            = 'WLS';
else
    weights_flag    = 0;
    type            = 'LS';
end


% Estimation
% =============
if (weights_flag)
    % loop for convergence (with weighting matrix)
    % =============================================
    while (1)
        iter    = iter + 1;
        h       = x.*exp( -x2/(2*s) )/s;
        H       = h.*(x2/(2*s^2) - 1/s);
        HTW     = H'*W;
        e       = inv( HTW * H ) * HTW * (y-h);
        s       = s + e;
        control = e*e;
        if ( control > (last_cnt * thresh) )
            break;
        else
            last_cnt = control;
        end
    end

    % summarize results
    h       = x.*exp( -x2/(2*s) )/s;
    H       = h.*(x2/(2*s^2) - 1/s);
    HTW     = H'*W;
    G       = inv( HTW * H ) * HTW;
    err     = (y-h);
    VAR     = G * var(err) * (G');
else

    % loop for convergence (without a weighting matrix)
    % ==================================================
    while (1)
        iter    = iter + 1;
        h       = x.*exp( -x2/(2*s) )/s;
        H       = h.*(x2/(2*s^2) - 1/s);
        HT      = H';
        control = inv( HT * H );
        s       = s + control * HT * (y-h);
        if ( control>(last_cnt * thresh) )
            break;
        else
            last_cnt = control;
        end
    end
    
	% summarize results
    h       = x.*exp( -x2/(2*s) )/s;
    H       = h.*(x2/(2*s^2) - 1/s);
    err     = (y-h);
    VAR     = inv( (H') * H );
end


% finish summarizing results
% ============================
result.s    = s;
result.VAR  = VAR;
result.RMS  = sqrt( (err')*err/ (x(2)-x(1))^2 / (length(err)-1) );
result.iter = iter;
result.type = type;

% plot distribution if asked for
% ===============================
if (nargin>3)
    if ishandle( hAx )
        plot_rayleigh( x,result,hAx,4 );
    else
        figure;
        plot_rayleigh( x,result1,gca,4 );
    end
end