Parameter estimation of Rayleigh distribution
[s,pdf,err] = vb_rayleigh_estimate(Nhist,level,jmax)
Nhist : histogram
level : level of histogram
jmax : Max index of histogram to estimation
s : Variance of Rayleigh distribution
pdf : Estimated Rayleigh histogram
err : Estimated error
背景ノイズが従うRayleigh分布の分散を推定
Rayleigh分布
x >= 0
P(x) = a * x * exp( - a * x^2/2 )
Cummulative probability
Q(x) = ∫ P(x) = 1 - exp( - a * x^2/2 )
= 1 - exp( - (x/s)^2/2 )
P(s) : Max of P(x)
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [s,pdf,err] = vb_rayleigh_estimate(Nhist,level,jmax) 0002 % Parameter estimation of Rayleigh distribution 0003 % [s,pdf,err] = vb_rayleigh_estimate(Nhist,level,jmax) 0004 % 0005 % Nhist : histogram 0006 % level : level of histogram 0007 % jmax : Max index of histogram to estimation 0008 % 0009 % s : Variance of Rayleigh distribution 0010 % pdf : Estimated Rayleigh histogram 0011 % err : Estimated error 0012 % 0013 % 背景ノイズが従うRayleigh分布の分散を推定 0014 % Rayleigh分布 0015 % x >= 0 0016 % P(x) = a * x * exp( - a * x^2/2 ) 0017 % Cummulative probability 0018 % Q(x) = ∫ P(x) = 1 - exp( - a * x^2/2 ) 0019 % = 1 - exp( - (x/s)^2/2 ) 0020 % P(s) : Max of P(x) 0021 % 0022 % Copyright (C) 2011, ATR All Rights Reserved. 0023 % License : New BSD License(see VBMEG_LICENSE.txt) 0024 0025 if nargin < 3, jmax = length(Nhist); end; 0026 0027 % Parameter estimation of Rayleigh distribution 0028 x = level(1:jmax); 0029 y = Nhist(1:jmax); 0030 0031 [Hmax ,imax]= max(y); 0032 0033 % If max point is at imax = 1, this distribution is not Rayleigh 0034 if imax == 1, imax = 3; end; 0035 0036 xmax = level(imax); 0037 0038 % Cumulative histogram 0039 ysum = cumsum(y); 0040 fsum = 1-exp(-0.5*(x/xmax).^2); 0041 0042 % proportinal constant estimation ( ysum = A * fsum ) 0043 A = sum(ysum.*fsum)/sum(fsum.^2); 0044 % 1 - (Cumulative histogram) 0045 % = exp( - (x/s)^2/2 ) 0046 ysum = 1 - ysum/A; 0047 0048 % s : Variance of Rayleigh distribution 0049 % ysum = exp( - (x/s)^2/2 ) 0050 s = 1/sqrt(abs(2*(log(ysum(imax))/x(imax)^2))); 0051 0052 % Check Rayleigh distribution 0053 % pdf : estimated Rayleigh histogram 0054 0055 xx = [0 x]; 0056 % 1 - (Cumulative histogram) 0057 fsum = exp(-0.5*(xx/s).^2); 0058 % Estimated Rayleigh distribution 0059 pdf = A*[fsum(1:end-1) - fsum(2:end)]; 0060 0061 ix = find( x >= 4*s ); 0062 0063 if isempty(ix) 0064 kmax = jmax; 0065 else 0066 kmax = ix(1); 0067 end 0068 0069 %kmax = min(4*imax, jmax); 0070 err = sum(abs(y(1:kmax) - pdf(1:kmax)))/sum(abs(y(1:kmax)));