function [u,sig2,result] = fit_ML_normal( x,hAx )
% fit_ML_normal - Maximum Likelihood fit of the normal distribution of i.i.d. samples!.
%                  Given the samples of a normal distribution, the PDF parameter is found
%
%    fits data to the probability of the form: 
%        p(r) = sqrt(1/2/pi/sig^2)*exp(-((r-u)^2)/(2*sig^2))
%    with parameters: u,sig^2
%
% format:   result = fit_ML_normal( x,hAx )
%
% input:    x   - vector, samples with normal distribution to be parameterized
%           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
%                      sig^2,u          - fitted parameters
%                      CRB_sig2,CRB_u   - Cram?r-Rao Bound for the estimator value
%                      RMS              - RMS error of the estimation 
%                      type             - 'ML'
%
% example:  fit_ML_normal( randn(1,10000)*3 - 1 )
%               or
%           figure; 
%           fit_ML_normal( randn(1,10000)*3 - 1,gca );
%

%
% Algorithm
% ===========
%
% We use the ML algorithm to estimate the PDF from the samples.
% The normal destribution is given by:
%
%    p(x;u,sig^2) = sqrt(1/2/pi/sig^2)*exp(-((x-u).^2)/(2*sig^2))
%
%    where x are the samples which distribute by the function p(x;u,sig^2)
%            and are assumed to be i.i.d !!!
%
% The ML estimator is given by:
%
%    a         = parameters vector = [u,sig^2]
%    f(Xn,a)   = sqrt(1/2/pi/sig^2)*exp(-((Xn-u).^2)/(2*sig^2))
%    L(a)      = f(X,a) = product_by_n( f(Xn,a) )
%              = (2*pi*sig^2)^(-N/2)*exp(-sum((Xn-u)^2)/(2*sig^2))
%    log(L(a)) = -N/2*log(2*pi) - N/2*log(sig^2) - sum((Xn-u)^2)/(2*sig^2)
%
%    The maximum likelihood point is found by the derivative of log(L(a)) with respect to "a":
%
%    diff(log(L(a)),u) = -2*sum(Xn-u)/(2*sig^2) = N/sig^2 * ( u - sum(Xn)/N )
%                    = J(u) * (u_estimation - u)  
%    diff(log(L(a)),sig^2) = -N/(2*sig^2) + sum((Xn-u)^2)/(2*sig^4)
%                          = N/(2*sig^4) * ( sum((Xn-u)^2)/N - sig^2 )
%                          = J(sig^2) * (sig^2_estimator - sig^2)
%
%    Therefore, the (efficient) estimators are given by:
%
%               u     = sum(Xn)/N
%               sig^2 = sum((Xn-u)^2)/N
%
%    The Cram?r-Rao Bounds for these estimator are:
%
%               VAR( u )     = 1/J(u)     = (sig^2) / N
%               VAR( sig^2 ) = 1/J(sig^2) = (2*sig^4) / N
%
%    NOTE: the ML estimator does not detect a deviation from the model.
%          therefore, check the RMS value !
%

if (nargin<1)
    error( 'fit_ML_normal - insufficient input arguments' );
end

% Estimation
% =============
x       = x(:);                 % should be column vectors !
N       = length(x);
u       = sum(x)/N;       
sig2    = (x-u)'*(x-u)/N;
CRB_u   = sig2 / N;
CRB_sig2= (2*sig2^2) / N;
[n,x_c] = hist( x,100 );
n       = n / sum(n*abs(x_c(2)-x_c(1)));
y       = sqrt(1/2/pi/sig2)*exp(-((x_c-u).^2)/(2*sig2));
RMS     = sqrt( (y-n)*((y-n)')/ (x_c(2)-x_c(1))^2 / (length(x_c)-1) );

% finish summarizing results
% ==================================
 result = struct( 'u',u,'sig2',sig2,'CRB_u',CRB_u,'CRB_sig2',CRB_sig2,'RMS',RMS,'type','ML' );

% plot distribution if asked for
% ==================================
if (nargin>1)
    xspan = linspace(min(x),max(x),100);
    if ishandle( hAx )
        plot_normal( xspan,result,hAx,1,14 );
    else
        figure;
        plot_normal( xspan,result,gca,1,14 );
    end
end