tab_3_3_delta_0_for_chi_2

% Table 3.3 to 3.5 page 68
 determine delta_0 for chi^2 distribution

 Wolfgang Förstner 2015
 last changes: Susanne Wenzel 09/16
 wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de

0001 %% Table 3.3 to 3.5 page 68
0002 % determine delta_0 for chi^2 distribution
0003 %
0004 % Wolfgang Förstner 2015
0005 % last changes: Susanne Wenzel 09/16
0006 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0007 
0008 
0009 disp('---- determine lower bounds delta for non-central chi-distribution ----')
0010 
0011 %% set parameters
0012 
0013 % fix alpha's
0014 alphas = [0.003:-0.00005:0.002];
0015 Nalpha = length(alphas);
0016 
0017 % set beta's
0018 betas = [0.8] %,0.9,0.95];
0019 
0020 % fix degrees of freedom
0021 dof =[2]; %,3,4,8,10,20,50];
0022 Ndof = length(dof);
0023 
0024 
0025 %% complete table including left and top margin
0026 table_alpha_d = zeros(1+Nalpha,1+Ndof);
0027 table_alpha_d(2:Nalpha+1,1) = alphas;
0028 table_alpha_d(1,2:1+Ndof) = dof;
0029 
0030 % for all beta's
0031 for ib = 1:length(betas);
0032 
0033     beta_0 = betas(ib);
0034     
0035     % for all alpha's
0036     for ia = 1:Nalpha
0037         
0038         alpha = alphas(ia);
0039         
0040         % for all dof's
0041         for idf = 1:length(dof)
0042         
0043             df = dof(idf);
0044             
0045             % critical value for chi-square test
0046             c = chi2inv(1-alpha,df);
0047             
0048             % set boundaries for searching for delta_0
0049             mindelta = 0.1;
0050             maxdelta = 6*df;
0051             meandelta = (maxdelta+mindelta)/2;
0052             
0053             % use noncentral chi-square cumulative distribution function to
0054             % determine boundaries
0055             minpower = 1 - ncx2cdf(c,df,mindelta^2);
0056             maxpower = 1 - ncx2cdf(c,df,maxdelta^2);
0057             
0058             % search binary, half intervals until right and left interval
0059             % ends have the same power
0060             while maxpower-minpower > 10^-10
0061                 
0062                 d_power=maxpower-minpower;
0063                 meanpower= 1-ncx2cdf(c,df,meandelta^2);
0064                 
0065                 if meanpower > beta_0
0066                     maxdelta  =  meandelta;
0067                     maxpower = 1-ncx2cdf(c,df,maxdelta^2);
0068                     meandelta = (maxdelta+mindelta)/2;
0069                 else
0070                     mindelta = meandelta;
0071                     minpower = 1-ncx2cdf(c,df,mindelta^2);
0072                     meandelta = (maxdelta+mindelta)/2;
0073                 end                
0074                 
0075             end
0076             
0077             table_alpha_d(ia+1,idf+1) = meandelta;
0078             
0079         end
0080     end
0081     beta = beta_0 %#ok<*NOPTS>
0082     t_alpha_0_d = table_alpha_d
0083 
0084 end