matlab/histograms.m

   1 clear;
   2
   3 % ---------- Parameters for the plots ----------
   4
   5 rows    = 3;                            % Number of rows in each figure
   6 columns = 4;                            % Number of columns
   7 bars    = 20;                           % Number of histogram bars
   8 text_offset = 0.2;                      % The text below each graph
   9                                         % will be placed this far
  10                                         % below the X-axis.
  11 factor_minimization_successful     = 1; % Condition the graphs on
  12                                         % the minimization outcome
  13 factor_covariance_step_successful  = 1; % Condition the graphs on
  14                                         % the covariance step outcome
  15 factor_no_covariance_step_warnings = 0; % Condition the graphs on
  16                                         % the presence of warning
  17                                         % messages in the
  18                                         % covariance step
  19 factor_no_estimate_near_boundary   = 1; % Condition the graphs on
  20                                         % the presence of parameter
  21                                         % estimates near boundaries
  22
  23 % The conditioning alternatives above create times-two extra graphs
  24 % per option set to 1. In total, 16 multiples of histograms may be
  25 % created. In reality, some combinations of errors an warnings are
  26 % impossible or very unlikely which reduce the dimensions
  27 % substancially.
  28
  29 confidence_limits = [2.5, 97.5 ];       % The percentiles to get
  30                                         % bootstrap estimates for
  31
  32 percentile_bca_limit = 1000;
  33 % If the number of available bootstrap samples, given the
  34 % conditions specified above, is lower than this value, a normal
  35 % percentile approach is used even if a BCa was requested.
  36
  37 % ---------- Autogenerated code below ----------
  38
  39 use_bca = 1;                            % Was a BCa-type of
  40                                         % bootstrap run?
  41 bs_samples = 2000;                      % Number of bootstrap samples
  42 jk_samples = 36;                        % Number of (BCa) jackknife samples
  43
  44 str_format = '%30s';
  45
  46 col_names = [ sprintf(str_format,'Significant Digits');
  47               sprintf(str_format,'OFV');
  48               sprintf(str_format,'par1');
  49               sprintf(str_format,'par2');
  50               sprintf(str_format,'par3');
  51               sprintf(str_format,'par4');
  52               sprintf(str_format,'par5');
  53               sprintf(str_format,'par6');
  54               sprintf(str_format,'par7');
  55               sprintf(str_format,'par8');
  56               sprintf(str_format,'par9')];
  57
  58 fixed_columns = [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1 ];
  59 same_columns = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ];
  60
  61 minimization_successful_col    = 2;     % Column number for the
  62                                         % minimization sucessful flag
  63 covariance_step_successful_col = 3;     % As above for cov-step warnings
  64 covariance_step_warnings_col   = 4;     % etc
  65 estimate_near_boundary_col     = 5;     % etc
  66
  67 not_data_cols = 5;                      % Number of columns in the
  68                                         % beginning that are not
  69                                         % parameter estimates.
  70
  71 filename = 'bootstrap_rr.csv';
  72
  73 % ---------- End autogenerated code ----------
  74
  75 % The rest of the code does not change between runs
  76
  77 bootstrap_log = load(filename);
  78
  79 [junk, n_columns] = size(fixed_columns);
  80
  81 use_rows     = ones(bs_samples,1);
  82 combinations = [NaN];
  83 tmp = [];
  84 tmp_comb = [];
  85
  86 % Create the conditioned (sub)data sets from the bootstrap data
  87 if( factor_minimization_successful )
  88   tmp(:,1) = bootstrap_log(2:bs_samples+1,minimization_successful_col) == 1;
  89   tmp(:,2) = bootstrap_log(2:bs_samples+1,minimization_successful_col) == 0;
  90   use_rows = tmp;
  91   combinations = [1,0];
  92 end
  93 tmp_comb = [];
  94 if( factor_covariance_step_successful )
  95   [a,b] = size( use_rows );
  96   for i = 1:b,
  97     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_successful_col) == 1);
  98     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_successful_col) == 0);
  99     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);1];
 100     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);0];
 101   end
 102   use_rows     = tmp;
 103   combinations = tmp_comb;
 104 end
 105 tmp_comb = [];
 106 if( factor_no_covariance_step_warnings )
 107   [a,b] = size( use_rows );
 108   for i = 1:b,
 109     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_warnings_col) == 0);
 110     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_warnings_col) == 1);
 111     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);0];
 112     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);1];
 113   end
 114   use_rows = tmp;
 115   combinations = tmp_comb;
 116 end
 117 tmp_comb = [];
 118 if( factor_no_estimate_near_boundary )
 119   [a,b] = size( use_rows );
 120   for i = 1:b,
 121     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,estimate_near_boundary_col) == 0);
 122     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,estimate_near_boundary_col) == 1);
 123     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);0];
 124     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);1];
 125   end
 126   use_rows = tmp;
 127   combinations = tmp_comb;
 128 end
 129
 130 [a,b] = size( use_rows );
 131 n_runs = sum( use_rows )
 132 orig_est = bootstrap_log(1,not_data_cols+1:n_columns);
 133 all_data = bootstrap_log(2:bs_samples+1,not_data_cols+1: ...
 134                          n_columns);
 135 if use_bca
 136         jk_data  = bootstrap_log(bs_samples+2:bs_samples+jk_samples+1, ...
 137                          not_data_cols+1: n_columns);
 138 end
 139
 140 % Loop over the possible subsets
 141 ci = [];
 142 for j = 1:b
 143   % Only those subsets that have some runs are of interest (obviously)
 144   if n_runs( j ) > 0
 145     row = 1;
 146     col = 1;
 147     [R,C] = find( use_rows(:,j) );
 148
 149     % Create subset
 150     data = all_data(R,:);
 151
 152     lower = [];
 153     upper = [];
 154     % Calculate the confidence limits according to each method
 155     if ( use_bca == 1 ) & ( n_runs(j) > percentile_bca_limit )
 156       [lower, upper] = bca( orig_est, data, jk_data, ...
 157                             confidence_limits(1),  confidence_limits(2) );
 158       ci = [ ci, [[combinations(:,j),combinations(:,j)];[ lower', upper']] ];
 159     else
 160       lower = prctile(data,confidence_limits(1));
 161       upper = prctile(data,confidence_limits(2));
 162       ci = [ ci, [[combinations(:,j),combinations(:,j)];[ lower', upper' ]] ];
 163     end
 164
 165     % Loop over all parameters
 166     for i = 1:n_columns,
 167       if fixed_columns( i ) == 1 || same_columns( i ) == 1
 168         continue;
 169       end
 170       if col == (columns+1)
 171         col = 1;
 172         row = row + 1;
 173       end
 174       if row == (rows+1) || i == 1
 175         h  = figure;
 176         col = 1;
 177         row = 1;
 178       end
 179
 180       % Prepare plot
 181       nplot = col+(row-1)*columns;
 182       subplot(rows, columns, nplot);
 183
 184       % Create a histogram
 185       hist(data(:,i),bars);
 186       hold on;
 187
 188       ax = axis;
 189       % Add lines for the lower and upper confidence limits as they are
 190       % determined by the bootstrap
 191       h = line( [lower(i), lower(i)], ax(3:4) );
 192       h = line( [upper(i), upper(i)], ax(3:4) );
 193
 194       % ----------- Add some text and figures ----------
 195
 196       % Upper likelihood-profiling confidence limit
 197       y_text = ax(3) - (ax(4)-ax(3))*text_offset;
 198       text( upper(i), y_text, num2str(upper(i)) );
 199
 200       % Lower likelihood-profiling confidence limit
 201       text( lower(i), y_text, num2str(lower(i)) );
 202
 203       % Add a title to each plot
 204       title( col_names( i, : ) );
 205
 206       hold off;
 207       col = col + 1;
 208     end
 209
 210   end
 211 end
 212
 213 save 'matlab_ci.prn' ci -ascii;
 214
 215
 216
 217
 218
 219