vb_gamma_dist_param

 Estimate threshold by approximating histgram by gamma distribution
    [yth] = vb_gamma_dist_param(y, p_val)
    [yth, A ] = vb_gamma_dist_param(y, p_val)
    [yth, A ,hy ,ylist] = vb_gamma_dist_param(y, p_val)
    
 y    : signal
 p_val: P-value corresponding to the threshold

 --- Output
 yth   : threshold value
 A     : gamma distribution parameter
         log(p(y)) = A(1)*y + A(2)*log(y) + A(3)
 hy : histgram of y
 ylist : list of y corresponding to hy

 2009-09-05 Masa-aki Sato

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

0001 function    [yth, A ,hy ,ylist] = vb_gamma_dist_param(y, p_val)
0002 % Estimate threshold by approximating histgram by gamma distribution
0003 %    [yth] = vb_gamma_dist_param(y, p_val)
0004 %    [yth, A ] = vb_gamma_dist_param(y, p_val)
0005 %    [yth, A ,hy ,ylist] = vb_gamma_dist_param(y, p_val)
0006 %
0007 % y    : signal
0008 % p_val: P-value corresponding to the threshold
0009 %
0010 % --- Output
0011 % yth   : threshold value
0012 % A     : gamma distribution parameter
0013 %         log(p(y)) = A(1)*y + A(2)*log(y) + A(3)
0014 % hy : histgram of y
0015 % ylist : list of y corresponding to hy
0016 %
0017 % 2009-09-05 Masa-aki Sato
0018 %
0019 % Copyright (C) 2011, ATR All Rights Reserved.
0020 % License : New BSD License(see VBMEG_LICENSE.txt)
0021 
0022 debug=0;
0023 
0024 % y0 : upper limit to estimate mean(y) : cumsum(p(y)) = p_0
0025 p_0  = 2/3; %p_0  = 3/4;
0026 rate = 0.1;
0027 
0028 T = length(y);
0029 y = abs(y(:));
0030 
0031 Nbin = min(5000,fix(T/10));
0032 
0033 % histgram
0034 [hy ,ylist] = hist(y,Nbin);
0035 
0036 % --- DEBUG
0037 if debug==1
0038     ymin = ylist(1)
0039 end
0040 
0041 % cummulative probality
0042 py = cumsum(hy)/T;
0043 
0044 % estimate mean & variance excluding outlier
0045 ix = find( py > p_0 );
0046 y0 = ylist(ix(1));
0047 jj = find(y < y0);
0048 ym = mean(y(jj));
0049 yd = sqrt(mean((y(jj)-ym).^2));
0050 
0051 % upper & lower limit of y when estimating gamma distribution
0052 ylow = max(ym-yd, eps);
0053 yup  = ym + 2*yd;
0054 %yup  = y0;
0055 
0056 if debug==1
0057     fprintf('ym=%g, yup=%g, ylow=%g\n',ym, y0,ylow)
0058 end
0059 
0060 ii = find((ylist <= yup));
0061 %ii = find( (ylist >= ylow) & (ylist <= yup));
0062 %ii = find( (ylist >= ylow) );
0063 %ii = find( (ylist >= ylow) & (ylist <= y0));
0064 %ii = find( (ylist >= eps) & (ylist <= yup));
0065 
0066 z = hy(ii);
0067 x = ylist(ii);
0068 
0069 % exclude low value of histgram
0070 ix = find( z > max(z)*rate );
0071 
0072 if debug==1
0073     Nhist = length(ix)
0074 end
0075 
0076 x = x(ix);
0077 z = z(ix);
0078 x = x(:);
0079 z = z(:);
0080 N = length(x);
0081 
0082 % parameter estimation for Gamma PDF
0083 Z  = log(z);
0084 X  = [ x , log(x), ones(N,1)];
0085 XX = X' * X;
0086 XZ = X' * Z;
0087 A  = XX \ XZ;
0088 
0089 %py = exp( A(1)*x + A(2)*log(x) + A(3));
0090 %errA = mean(abs(z - py));
0091 errA = mean(abs(Z - X * A));
0092 
0093 % parameter estimation for exponential PDF
0094 X = [ x , ones(N,1)];
0095 XX = X' * X;
0096 XZ = X' * Z;
0097 B = XX \ XZ;
0098 
0099 %py = exp( B(1)*x + B(2));
0100 %errB = mean(abs(z - py));
0101 errB = mean(abs(Z - X * B));
0102 
0103 if A(1) > 0 || A(2) < 0 || errA > errB
0104     fprintf('ym=%g, ylow=%g\n',ym/max(y), ylow/max(y))
0105     fprintf('Max Y for PDF = %g\n', yup/max(y))
0106     fprintf('Exponential PDF seems to be better than Gamma PDF\n')
0107     fprintf('Estimated Gamma PDF = exp( %g *y + %g *log(y) )\n',A(1:2))
0108     fprintf('Estimated Exp   PDF = exp( %g *y )\n',B(1))
0109     A = [B(1) 0 B(2)];
0110 end
0111 
0112 if A(1) > 0
0113     fprintf('gamma pdf estimation error\n')
0114 end
0115 
0116 % estimated gamma distribution
0117 if A(2)==0
0118     ii = find( ylist >= 0);
0119     
0120     yl = ylist(ii);
0121     p  = exp( A(1)*yl + A(3));
0122 else
0123     ii = find( ylist > 0);
0124     
0125     yl = ylist(ii);
0126     p  = exp( A(1)*yl + A(2)*log(yl) + A(3));
0127 end
0128 
0129 psum = cumsum(p)/sum(p);
0130 
0131 % threshold corresponding to p_val
0132 ii = find(psum > 1 - p_val);
0133 yth = yl(ii(1));
0134 
0135 % x = gaminv(p,a,b)
0136 % �K���}�z���A�m�p �^�l x �v
0137 %
0138 % gamma(x,a,b) �exp(-x/b) x^(a-1) = exp( -x/b + (a-1)*log(x) )
0139 % pdf  = exp( A(1)*x + A(2)*log(x) );
0140 % p = int_0^x dt gamma(t,a,b)
0141 % a = A(2) + 1  : A(2) = a - 1
0142 % b = - 1/A(1)  : A(1) = -1/b
0143 %
0144 
0145 if debug == 2
0146     ix = find( (y >= ylow) & (y <= yup));
0147     B = gamfit(y(ix));
0148     A = [ - 1/B(2) , B(1) - 1, 1];
0149     yth0 = yth;
0150     yth  = gaminv(1-p_val, A(2)+1, -1/A(1))/max(y);
0151     dif_yth  = (yth - yth0)
0152 end
0153 
0154 return
0155 
0156 
0157 % y = C * py
0158 
0159 C = (y'*py) / sum(py.^2);
0160 
0161 py = C * py;