vb_gaussian_estimate

 Parameter estimation of Gaussian distribution
  [xmax,s,C,err] = vb_gaussian_estimate(Nhist,level,jmax)

 Nhist : histogram
 level : level of histogram
 jmax  : Max index of histogram to estimation

 xmax: Peak point of histogram : center of Gaussian
 s   : Variance 
 C   : normalization constant
 err : Estimated error

 pdf(x) = C * exp( - ((x - xmax)/s)^2/2 )
 
 P(x) = exp( - (x/s)^2/2 )
      = exp( - A * x^2 )
  s^2 = 1/(2*A)

 dP/dA = - x^2 * P(x)
  E    = <( y - P(x) )^2>
    dE/dA  = < (y - P) * x^2 * P>
 d^2E/dA^2 = < (x^2 * P)^2 > + ...
 a_new = a + da
 dE/da(a_new) = dE/da + da * d^2E/da^2 = 0

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [xmax,s,C,err] = vb_gaussian_estimate(Nhist,level,jmax)
0002 % Parameter estimation of Gaussian distribution
0003 %  [xmax,s,C,err] = vb_gaussian_estimate(Nhist,level,jmax)
0004 %
0005 % Nhist : histogram
0006 % level : level of histogram
0007 % jmax  : Max index of histogram to estimation
0008 %
0009 % xmax: Peak point of histogram : center of Gaussian
0010 % s   : Variance
0011 % C   : normalization constant
0012 % err : Estimated error
0013 %
0014 % pdf(x) = C * exp( - ((x - xmax)/s)^2/2 )
0015 %
0016 % P(x) = exp( - (x/s)^2/2 )
0017 %      = exp( - A * x^2 )
0018 %  s^2 = 1/(2*A)
0019 %
0020 % dP/dA = - x^2 * P(x)
0021 %  E    = <( y - P(x) )^2>
0022 %    dE/dA  = < (y - P) * x^2 * P>
0023 % d^2E/dA^2 = < (x^2 * P)^2 > + ...
0024 % a_new = a + da
0025 % dE/da(a_new) = dE/da + da * d^2E/da^2 = 0
0026 %
0027 % Copyright (C) 2011, ATR All Rights Reserved.
0028 % License : New BSD License(see VBMEG_LICENSE.txt)
0029 
0030 %clear all
0031 %%% DEBUG
0032 %j0 = 20;
0033 %jmax  = 100;
0034 %level = (0:jmax)/jmax;
0035 %Nhist = 2* exp( - 3 * (level - level(j0)).^2 );
0036 
0037 dmin  = 1e-10;
0038 Niter = 1000;
0039 
0040 if nargin < 3, jmax = length(Nhist); end;
0041 
0042 Nhist = Nhist(1:jmax);
0043 level = level(1:jmax);
0044 
0045 [ymax ,imax]= max(Nhist);
0046 xmax = level(imax);
0047 
0048 % Parameter estimation of Gaussian distribution
0049 x = level(imax:jmax) - xmax;
0050 y = Nhist(imax:jmax)/ymax;
0051 
0052 % Initial estimation by using log(y)
0053 % ( log(y) = A * x^2 )
0054 A = sum(log(y).* x.^2)/sum(x.^4);
0055 dmin = abs(A) * dmin;
0056 
0057 % Quasi-Newton iteration
0058 for n=1:Niter
0059     P   = exp( - A * x.^2 );
0060     dE  = sum((y - P) .* x.^2 .* P);
0061     ddE = sum((x.^2 .* P).^2)  ;
0062     dA  = - dE/ddE;
0063     
0064     if abs(dA) < dmin, break; end;
0065     
0066     A   = A + dA;
0067 end
0068 
0069 % y = C * P(x)
0070 P  = exp( - A * x.^2 );
0071 C  = sum(y.* P)/sum(P.^2);
0072 
0073 % s : Variance of Gaussian distribution
0074 %     s^2 = 1/(2*A)
0075 s = 1/sqrt(abs(2*A));
0076 
0077 % Normalization constant
0078 C = C * ymax ;
0079 
0080 ix = find( level >= (xmax + 3*s) );
0081 
0082 if isempty(ix)
0083     kmax = jmax;
0084 else
0085     kmax = min(ix(1) , jmax);
0086 end
0087 
0088 % pdf  : Gaussian histogram
0089 pdf = C * exp(-0.5*((level(1:kmax) - xmax)/s).^2);
0090 
0091 err = sum(abs(Nhist(1:kmax) - pdf))/sum(abs(Nhist(1:kmax)));
0092 
0093 return
0094 
0095 plot(level,Nhist)
0096 hold on
0097 plot(level,pdf,'-r')
0098