statanalysis_yp

 This function calculates the mean value, the standard deviation, and the
 relative accuracy of the mean value. The function should exhibit a
 stationnary region (plateau)

   INPUT :

       - y: noisy signal [1,m]
       - yt: minimum relative accuracy of the mean value [1,1]

   OUTPUT :

       - iy: valid index [1,p]
       - y: noisy signal [1,p]
       - ym: mean value [1,1]
       - ys: standard deviation [1,1]
       - yr: relative accuracy of the mean value [1,1]
       - iyr: rejected index [1,m-p]

 by Y. Peysson DRFC-CEA <yves.peysson@cea.fr> and J. Decker DRFC-CEA <joan.decker@cea.fr>

0001 function [iy,y,ym,ys,yr,iyr] = statanalysis_yp(yin,yt,display_mode)
0002 %
0003 % This function calculates the mean value, the standard deviation, and the
0004 % relative accuracy of the mean value. The function should exhibit a
0005 % stationnary region (plateau)
0006 %
0007 %   INPUT :
0008 %
0009 %       - y: noisy signal [1,m]
0010 %       - yt: minimum relative accuracy of the mean value [1,1]
0011 %
0012 %   OUTPUT :
0013 %
0014 %       - iy: valid index [1,p]
0015 %       - y: noisy signal [1,p]
0016 %       - ym: mean value [1,1]
0017 %       - ys: standard deviation [1,1]
0018 %       - yr: relative accuracy of the mean value [1,1]
0019 %       - iyr: rejected index [1,m-p]
0020 %
0021 % by Y. Peysson DRFC-CEA <yves.peysson@cea.fr> and J. Decker DRFC-CEA <joan.decker@cea.fr>
0022 %
0023 if nargin < 3,
0024     display_mode = 0;% default value, no display
0025 end
0026 %
0027 if nargin < 2,
0028     yt = NaN;% default value, two times the standard deviation
0029 end
0030 %
0031 nyin = length(yin);
0032 y = yin;
0033 ny = nyin;
0034 %
0035 ymask = true(1,ny);% mask of valid indices
0036 %
0037 ym = mean(y);
0038 ys = std(y);
0039 yd = ny*ys/ym;
0040 yr = Inf;
0041 %
0042 while ny > 1,%abs(yr) > abs(yt),%
0043     %
0044     iyc = find(abs(y - ym) == max(abs(y - ym)),1,'first');% find value furthest from average
0045     %
0046     y(iyc) = [];
0047     %
0048     ymc = mean(y);
0049     ysc = std(y);
0050     ydc = ny*ysc/ymc;
0051     %
0052     yr = abs(yd - ydc);
0053     %
0054     if isnan(yt),
0055         yrmax = 3*ys/ym;
0056     else
0057         yrmax = yt;
0058     end
0059     %
0060     if yr <= yrmax,
0061         break
0062     end
0063     %
0064     % take away rejected index
0065     %
0066     yind = find(ymask);
0067     ymask(yind(iyc)) = false;
0068     ny = ny - 1;
0069     %
0070     ym = ymc;
0071     ys = ysc;
0072     yd = ny*ys/ym;
0073     %
0074 end
0075 %
0076 iy = find(ymask);
0077 iyr = find(~ymask);
0078 y = yin(ymask);
0079 %
0080 if display_mode >= 2,
0081     %
0082     figure(1),clf
0083     %
0084     graph1D_jd(1:nyin,yin/ym,0,0,'n','y/ym','',NaN,[0,nyin+1],NaN,'-','+','r',2,20+14i,gca,0.9,0.7,0.7);
0085     graph1D_jd(iy,y/ym,0,0,'','','',NaN,[0,nyin+1],NaN,'none','o','b',2,20+14i);
0086     %
0087     graph1D_jd([0,nyin+1],[1+yrmax,1+yrmax],0,0,'','','',NaN,[0,nyin+1],NaN,'--','none','k',2,20+14i);
0088     graph1D_jd([0,nyin+1],[1-yrmax,1-yrmax],0,0,'','','',NaN,[0,nyin+1],NaN,'--','none','k',2,20+14i);
0089     %
0090     graph1D_jd([0,nyin+1],[1+ys/ym,1+ys/ym],0,0,'','','',NaN,[0,nyin+1],NaN,'--','none','g',2,20+14i);
0091     graph1D_jd([0,nyin+1],[1-ys/ym,1-ys/ym],0,0,'','','',NaN,[0,nyin+1],NaN,'--','none','g',2,20+14i);
0092     %
0093 end
0094     
0095