matlab/histograms.m

   1 clear;
   2
   3 % ---------- Parameters for the plots ----------
   4
   5 rows    = 2;                            % Number of rows in each figure
   6 columns = 2;                            % Number of columns
   7 bars    = 20;                           % Number of histogram bars
   8
   9 only_print = NaN;                       % Use to make the script print
  10                                         % only one parameter
  11
  12 text_size = 10;                         % Text size in legend and
  13                                         % on axes
  14
  15 create_emf = 0;                         % Create .emf files. Only
  16                                         % possible on windows.
  17
  18 factor_minimization_successful     = 1; % Condition the graphs on
  19                                         % the minimization outcome
  20 factor_covariance_step_successful  = 0; % Condition the graphs on
  21                                         % the covariance step outcome
  22 factor_no_covariance_step_warnings = 0; % Condition the graphs on
  23                                         % the presence of warning
  24                                         % messages in the
  25                                         % covariance step
  26 factor_no_estimate_near_boundary   = 1; % Condition the graphs on
  27                                         % the presence of parameter
  28                                         % estimates near boundaries
  29
  30 % The conditioning alternatives above create times-two extra graphs
  31 % per option set to 1. In total, 16 multiples of histograms may be
  32 % created. In reality, some combinations of errors an warnings are
  33 % impossible or very unlikely which reduce the dimensions
  34 % substancially.
  35
  36 confidence_limits = [2.5, 97.5 ];       % The percentiles to get
  37                                         % bootstrap estimates for
  38
  39 percentile_bca_limit = 100;
  40 % If the number of available bootstrap samples, given the
  41 % conditions specified above, is lower than this value, a normal
  42 % percentile approach is used even if a BCa was requested.
  43
  44 % ---------- Autogenerated code below ----------
  45
  46 use_bca = 1;                            % Was a BCa-type of
  47                                         % bootstrap run?
  48 bs_samples = 2000;                      % Number of bootstrap samples
  49 jk_samples = 36;                        % Number of (BCa) jackknife samples
  50
  51 col_names = { 'Significant Digits',
  52               'OFV',
  53               'par1',
  54               'par2',
  55               'par3',
  56               'par4',
  57               'par5',
  58               'par6',
  59               'par7',
  60               'par8',
  61               'par9' };
  62
  63 fixed_columns = [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1 ];
  64 same_columns = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ];
  65
  66 minimization_successful_col    = 2;     % Column number for the
  67                                         % minimization sucessful flag
  68 covariance_step_successful_col = 3;     % As above for cov-step warnings
  69 covariance_step_warnings_col   = 4;     % etc
  70 estimate_near_boundary_col     = 5;     % etc
  71
  72 not_data_cols = 12;                     % Number of columns in the
  73                                         % beginning that are not
  74                                         % parameter estimates.
  75
  76 filename = 'bootstrap_rr.csv';
  77
  78 % ---------- End autogenerated code ----------
  79
  80 % The rest of the code does not change between runs
  81
  82 bootstrap_log = csvread(filename);
  83
  84 [junk, n_columns] = size(fixed_columns);
  85
  86 use_rows     = ones(bs_samples,1);
  87 combinations = [NaN];
  88 comb_names   = {NaN};
  89 short_names   = {NaN};
  90 tmp = [];
  91 tmp_comb  = [];
  92 tmp_short  = {};
  93 tmp_names = {};
  94
  95 % Create the conditioned (sub)data sets from the bootstrap data
  96 if( factor_minimization_successful )
  97   tmp(:,1) = bootstrap_log(2:bs_samples+1,minimization_successful_col) == 1;
  98   tmp(:,2) = bootstrap_log(2:bs_samples+1,minimization_successful_col) == 0;
  99   use_rows = tmp;
 100   combinations = [1,0];
 101   comb_names   = {'minimization_successful=1','minimization_successful=0'};
 102   short_names   = {'m1','m0'};
 103 end
 104 tmp_comb = [];
 105 tmp_short  = {};
 106 tmp_names = {};
 107 if( factor_covariance_step_successful )
 108   [a,b] = size( use_rows );
 109   for i = 1:b,
 110     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_successful_col) == 1);
 111     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_successful_col) == 0);
 112     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);1];
 113     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);0];
 114     n1 = {'covariance_step_successful=1'};
 115     n2 = {'covariance_step_successful=0'};
 116     tmp_names(:,(i-1)*2+1) = [comb_names(:,i);n1];
 117     tmp_names(:,(i-1)*2+2) = [comb_names(:,i);n2];
 118     tmp_short(:,(i-1)*2+1) = [short_names(:,i);{'c1'}];
 119     tmp_short(:,(i-1)*2+2) = [short_names(:,i);{'c0'}];
 120   end
 121   use_rows     = tmp;
 122   combinations = tmp_comb;
 123   comb_names   = tmp_names;
 124   short_names  = tmp_short;
 125 end
 126 tmp_comb = [];
 127 tmp_short  = {};
 128 tmp_names = {};
 129 if( factor_no_covariance_step_warnings )
 130   [a,b] = size( use_rows );
 131   for i = 1:b,
 132     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_warnings_col) == 0);
 133     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_warnings_col) == 1);
 134     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);0];
 135     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);1];
 136     n1 = {'covariance_step_warnings=0'};
 137     n2 = {'covariance_step_warnings=1'};
 138     tmp_names(:,(i-1)*2+1) = [comb_names(:,i);n1];
 139     tmp_names(:,(i-1)*2+2) = [comb_names(:,i);n2];
 140     tmp_short(:,(i-1)*2+1) = [short_names(:,i);{'w0'}];
 141     tmp_short(:,(i-1)*2+2) = [short_names(:,i);{'w1'}];
 142   end
 143   use_rows = tmp;
 144   combinations = tmp_comb;
 145   comb_names   = tmp_names;
 146   short_names  = tmp_short;
 147 end
 148 tmp_comb = [];
 149 tmp_short  = {};
 150 tmp_names = {};
 151 if( factor_no_estimate_near_boundary )
 152   [a,b] = size( use_rows );
 153   for i = 1:b,
 154     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,estimate_near_boundary_col) == 0);
 155     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,estimate_near_boundary_col) == 1);
 156     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);0];
 157     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);1];
 158     n1 = {'estimate_near_boundary=0'};
 159     n2 = {'estimate_near_boundary=1'};
 160     tmp_names(:,(i-1)*2+1) = [comb_names(:,i);n1];
 161     tmp_names(:,(i-1)*2+2) = [comb_names(:,i);n2];
 162     tmp_short(:,(i-1)*2+1) = [short_names(:,i);{'b0'}];
 163     tmp_short(:,(i-1)*2+2) = [short_names(:,i);{'b1'}];
 164   end
 165   use_rows = tmp;
 166   combinations = tmp_comb;
 167   comb_names   = tmp_names;
 168   short_names  = tmp_short;
 169 end
 170
 171 [a,b] = size( use_rows );
 172 n_runs = sum( use_rows );
 173 orig_est = bootstrap_log(1,not_data_cols+1:not_data_cols+n_columns);
 174 all_data = bootstrap_log(2:bs_samples+1,not_data_cols+1: ...
 175                          not_data_cols+n_columns);
 176 if use_bca
 177         jk_data  = bootstrap_log(bs_samples+3:bs_samples+jk_samples+2, ...
 178                          not_data_cols+1: not_data_cols+n_columns);
 179 end
 180
 181 % Create directory for figures
 182 mkdir('matlab_plots');
 183 chdir('matlab_plots');
 184
 185 lower = [];
 186 upper = [];
 187 ci = [];
 188 % loop the subsets
 189 for j = 1:b
 190   % Only those subsets that have some runs are of interest (obviously)
 191   if n_runs( j ) > 0
 192     row = 1;
 193     col = 1;
 194     [R,C] = find( use_rows(:,j) );
 195
 196     % Create subset
 197     data = all_data(R,:);
 198
 199     % Calculate the confidence limits according to each method
 200     if ( use_bca == 1 ) & ( n_runs(j) > percentile_bca_limit )
 201       [lower(j,:), upper(j,:)] = bca( orig_est, data, jk_data, ...
 202                                       confidence_limits(1),  confidence_limits(2) );
 203       ci = [ ci, [[combinations(:,j),combinations(:,j)];[ lower(j,:)', upper(j,:)']] ];
 204     else
 205       lower(j,:) = prctile(data,confidence_limits(1),1);
 206       upper(j,:) = prctile(data,confidence_limits(2),1);
 207       ci         = [ ci, [[combinations(:,j),combinations(:,j)];[ lower(j,:)', upper(j,:)' ]] ];
 208     end
 209   end
 210 end
 211
 212 fignum = 1;
 213 % Loop over all parameters
 214 for i = 1:n_columns,
 215     if fixed_columns( i ) == 1 || ...
 216         same_columns( i ) == 1 || ...
 217       ( ~isnan(only_print) && i ~= only_print )
 218         continue;
 219     end
 220     % Loop over the possible subsets
 221     row = 1;
 222     col = 1;
 223     for j = 1:b
 224         % Only those subsets that have some runs are of interest (obviously)
 225         if n_runs( j ) > 0
 226              [R,C] = find( use_rows(:,j) );
 227
 228 %             % Create subset
 229              data = all_data(R,:);
 230              data = data(:,i);
 231
 232 %             lower = [];
 233 %             upper = [];
 234 %             % Calculate the confidence limits according to each method
 235 %             if ( use_bca == 1 ) & ( n_runs(j) > percentile_bca_limit )
 236 %                 [lower, upper] = bca( orig_est(i), data, jk_data(:,i), ...
 237 %                     confidence_limits(1),  confidence_limits(2) );
 238 %                 ci = [ ci, [[combinations(:,j),combinations(:,j)];[ lower', upper']] ];
 239 %             else
 240 %                 lower = prctile(data,confidence_limits(1),1);
 241 %                 upper = prctile(data,confidence_limits(2),1);
 242 %                 ci = [ ci, [[combinations(:,j),combinations(:,j)];[ lower', upper' ]] ];
 243 %             end
 244
 245             if col == (columns+1)
 246                 col = 1;
 247                 row = row + 1;
 248             end
 249             if row == (rows+1) || j == 1
 250                 if row == (rows+1)
 251                     file = strrep( col_names( i ), ' ', '_' );
 252                     file = strrep( file, '/', '_' );
 253                     file = strrep( file, '\', '_' );
 254                     if create_emf == 1
 255                         saveas( gcf, cell2mat(file), 'emf' );
 256                     end
 257                     saveas( gcf, cell2mat(file), 'epsc2' );
 258                 end
 259                 h_base  = figure;
 260                 set( h_base, 'Visible', 'off' );
 261                 set( h_base, 'PaperType', 'A4');
 262                 set( h_base, 'PaperOrientation', 'landscape');
 263                 set( h_base, 'PaperPosition', [0.634517 0.634517 28.4084 19.715] );
 264
 265                 col = 1;
 266                 row = 1;
 267             end
 268
 269             % Prepare plot
 270             nplot = col+(row-1)*columns;
 271             subplot(rows, columns, nplot);
 272
 273             % Create a histogram
 274             hist(data,bars);
 275             h = get(gca,'Children');
 276             set(h,'EdgeColor','w');
 277             hold on;
 278
 279             ax = axis;
 280             ax_s = ax(2:2:4)-ax(1:2:3);
 281             ax_n(1:2:3) = ax(1:2:3)-0.01*ax_s;
 282             ax_n(2:2:4) = ax(2:2:4)+0.01*ax_s;
 283             axis(ax_n);
 284
 285             % Add lines for the lower and upper confidence limits as they are
 286             % determined by the bootstrap
 287             h = line( [lower(j,i), lower(j,i)], ax(3:4) );
 288             h = line( [upper(j,i), upper(j,i)], ax(3:4) );
 289
 290             if adjust_axes( i ) == 2 && not ( lower(j,i) == upper(j,i) )
 291               axis([(lower(j,i)-0.5*abs(upper(j,i)-lower(j,i))), ...
 292                     (upper(j,i)+0.5*abs(upper(j,i)-lower(j,i))), ...
 293                     ax(3:4)]);
 294             end
 295
 296             % ----------- Add some text and figures ----------
 297
 298             [LEGH,OBJH,OUTH,OUTM] = legend( sprintf('%14s%5.3g','lower limit = ',lower(j,i)), ...
 299                 sprintf('%14s%5.3g','upper limit = ',upper(j,i)), ...
 300                 sprintf('%14s%5.0f','n = ', n_runs( j ) ), ...
 301                 'Location', 'NorthEast' );
 302             legend boxoff;
 303             set( gca, 'FontSize', text_size);
 304             %      set( OBJH(1), 'HorizontalAlignment', 'right' );
 305             %      set( OBJH(2), 'HorizontalAlignment', 'right' );
 306             %      set( OBJH(3), 'HorizontalAlignment', 'right' );
 307             set( OBJH(4), 'Visible', 'off');
 308             set( OBJH(5), 'Visible', 'off');
 309             set( OBJH(7), 'Visible', 'off');
 310
 311             % Add a title to each plot
 312             [c1,c2] = size(short_names);
 313             name = [];
 314             %      name = [col_names( i ),'-'];
 315             for k = 1:c1
 316                 name = [ name, short_names(k,j)];
 317             end
 318             if( col == 1 && row == 1 )
 319                 title( cell2mat([col_names( i ),{' - '},name]), 'EdgeColor', 'none' );
 320             else
 321                 title( cell2mat(name), 'EdgeColor', 'none' );
 322             end
 323
 324             hold off;
 325             col = col + 1;
 326         end
 327         %   close;
 328     end
 329
 330     if adjust_axes( i ) == 1
 331         plot_children = get( gcf, 'Children' );
 332         plot_axes = [];
 333         max_axis = [];
 334         min_axis = [];
 335         for j = 1 : length(plot_children)
 336             if strcmp( get( plot_children(j), 'Type'), 'axes') && ...
 337                   ~strcmp( get( plot_children(j), 'Tag'), 'legend')
 338                 plot_axes = [plot_axes, plot_children(j)];
 339                 max_axis = max([max_axis;axis(plot_children(j))],[],1);
 340                 min_axis = min([min_axis;axis(plot_children(j))],[],1);
 341             end
 342         end
 343         axis(plot_axes,[min_axis(1),max_axis(2),min_axis(3),max_axis(4)]);
 344     end
 345
 346     file = strrep( col_names( i ), ' ', '_' );
 347     file = strrep( file, '/', '_' );
 348     file = strrep( file, '\', '_' );
 349     if create_emf == 1
 350         saveas( gcf, cell2mat(file), 'emf' );
 351     end
 352     saveas( gcf, cell2mat(file), 'epsc2' );
 353 end
 354 close all;
 355 chdir('..');
 356 save 'matlab_ci.txt' ci -ascii;
 357
 358 % % The rest of the code does not change between runs
 359
 360 % bootstrap_log = csvread(filename);
 361
 362 % [junk, n_columns] = size(fixed_columns);
 363
 364 % use_rows     = ones(bs_samples,1);
 365 % combinations = [NaN];
 366 % comb_names   = {NaN};
 367 % tmp = [];
 368 % tmp_comb  = [];
 369 % tmp_names = {};
 370
 371 % % Create the conditioned (sub)data sets from the bootstrap data
 372 % if( factor_minimization_successful )
 373 %   tmp(:,1) = bootstrap_log(2:bs_samples+1,minimization_successful_col) == 1;
 374 %   tmp(:,2) = bootstrap_log(2:bs_samples+1,minimization_successful_col) == 0;
 375 %   use_rows = tmp;
 376 %   combinations = [1,0];
 377 %   comb_names   = {'minimization_successful=1','minimization_successful=0'};
 378 % end
 379 % tmp_comb = [];
 380 % tmp_names = {};
 381 % if( factor_covariance_step_successful )
 382 %   [a,b] = size( use_rows );
 383 %   for i = 1:b,
 384 %     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_successful_col) == 1);
 385 %     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_successful_col) == 0);
 386 %     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);1];
 387 %     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);0];
 388 %     n1 = {'covariance_step_successful=1'};
 389 %     n2 = {'covariance_step_successful=0'};
 390 %     tmp_names(:,(i-1)*2+1) = [comb_names(:,i);n1];
 391 %     tmp_names(:,(i-1)*2+2) = [comb_names(:,i);n2];
 392 %   end
 393 %   use_rows     = tmp;
 394 %   combinations = tmp_comb;
 395 %   comb_names   = tmp_names;
 396 % end
 397 % tmp_comb = [];
 398 % tmp_names = {};
 399 % if( factor_no_covariance_step_warnings )
 400 %   [a,b] = size( use_rows );
 401 %   for i = 1:b,
 402 %     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_warnings_col) == 0);
 403 %     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,covariance_step_warnings_col) == 1);
 404 %     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);0];
 405 %     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);1];
 406 %     n1 = {'covariance_step_warnings=0'};
 407 %     n2 = {'covariance_step_warnings=1'};
 408 %     tmp_names(:,(i-1)*2+1) = [comb_names(:,i);n1];
 409 %     tmp_names(:,(i-1)*2+2) = [comb_names(:,i);n2];
 410 %   end
 411 %   use_rows = tmp;
 412 %   combinations = tmp_comb;
 413 %   comb_names   = tmp_names;
 414 % end
 415 % tmp_comb = [];
 416 % tmp_names = {};
 417 % if( factor_no_estimate_near_boundary )
 418 %   [a,b] = size( use_rows );
 419 %   for i = 1:b,
 420 %     tmp(:,(i-1)*2+1) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,estimate_near_boundary_col) == 0);
 421 %     tmp(:,(i-1)*2+2) = use_rows(:,i).*(bootstrap_log(2:bs_samples+1,estimate_near_boundary_col) == 1);
 422 %     tmp_comb(:,(i-1)*2+1) = [combinations(:,i);0];
 423 %     tmp_comb(:,(i-1)*2+2) = [combinations(:,i);1];
 424 %     n1 = {'estimate_near_boundary=0'};
 425 %     n2 = {'estimate_near_boundary=1'};
 426 %     tmp_names(:,(i-1)*2+1) = [comb_names(:,i);n1];
 427 %     tmp_names(:,(i-1)*2+2) = [comb_names(:,i);n2];
 428 %   end
 429 %   use_rows = tmp;
 430 %   combinations = tmp_comb;
 431 %   comb_names   = tmp_names;
 432 % end
 433
 434 % [a,b] = size( use_rows );
 435 % n_runs = sum( use_rows )
 436 % orig_est = bootstrap_log(1,not_data_cols+1:not_data_cols+n_columns);
 437 % all_data = bootstrap_log(2:bs_samples+1,not_data_cols+1: ...
 438 %                        not_data_cols+n_columns);
 439 % if use_bca
 440 %       jk_data  = bootstrap_log(bs_samples+2:bs_samples+jk_samples+1, ...
 441 %                        not_data_cols+1: not_data_cols+n_columns);
 442 % end
 443
 444 % % Loop over the possible subsets
 445 % ci = [];
 446 % fignum = 1;
 447 % for j = 1:b
 448 %   % Only those subsets that have some runs are of interest (obviously)
 449 %   if n_runs( j ) > 0
 450 %     row = 1;
 451 %     col = 1;
 452 %     [R,C] = find( use_rows(:,j) );
 453
 454 %     % Create subset
 455 %     data = all_data(R,:);
 456
 457 %     lower = [];
 458 %     upper = [];
 459 %     % Calculate the confidence limits according to each method
 460 %     if ( use_bca == 1 ) & ( n_runs(j) > percentile_bca_limit )
 461 %       [lower, upper] = bca( orig_est, data, jk_data, ...
 462 %                           confidence_limits(1),  confidence_limits(2) );
 463 %       ci = [ ci, [[combinations(:,j),combinations(:,j)];[ lower', upper']] ];
 464 %     else
 465 %       lower = prctile(data,confidence_limits(1),1);
 466 %       upper = prctile(data,confidence_limits(2),1);
 467 %       ci = [ ci, [[combinations(:,j),combinations(:,j)];[ lower', upper' ]] ];
 468 %     end
 469
 470 %     % Loop over all parameters
 471 %     for i = 1:n_columns,
 472 %       if fixed_columns( i ) == 1 || same_columns( i ) == 1
 473 %       continue;
 474 %       end
 475 %       if col == (columns+1)
 476 %       col = 1;
 477 %       row = row + 1;
 478 %       end
 479 %       if row == (rows+1) || i == 1
 480 %       if i ~= 1
 481 %         saveas( gcf, num2str(fignum), 'epsc2' );
 482 %           fignum = fignum + 1;
 483 %           close;
 484 %         end
 485 %       h  = figure;
 486 %       col = 1;
 487 %       row = 1;
 488 %       end
 489
 490 %       % Prepare plot
 491 %       nplot = col+(row-1)*columns;
 492 %       subplot(rows, columns, nplot);
 493
 494 %       % Create a histogram
 495 %       hist(data(:,i),bars);
 496 %       hold on;
 497
 498 %       ax = axis;
 499
 500 %       if col == 1 & row == 1
 501 %       x_text = ax(1) - (ax(2)-ax(1))*0.7;
 502 %       [c1,c2] = size(comb_names);
 503 %       for k = 1:c1
 504 %         y_text = ax(4) + (ax(4)-ax(3))*(0.1*k);
 505 %         text( x_text, y_text, comb_names(k,j) );
 506 %       end
 507 %       end
 508
 509 %       % Add lines for the lower and upper confidence limits as they are
 510 %       % determined by the bootstrap
 511 %       h = line( [lower(i), lower(i)], ax(3:4) );
 512 %       h = line( [upper(i), upper(i)], ax(3:4) );
 513
 514 %       % ----------- Add some text and figures ----------
 515
 516 %       % Upper likelihood-profiling confidence limit
 517 %       y_text = ax(3) - (ax(4)-ax(3))*text_offset;
 518 %       text( upper(i), y_text, num2str(upper(i)) );
 519
 520 %       % Lower likelihood-profiling confidence limit
 521 %       text( lower(i), y_text, num2str(lower(i)) );
 522
 523 %       % Add a title to each plot
 524 %       title( col_names( i, : ) );
 525
 526 %       hold off;
 527 %       col = col + 1;
 528 %     end
 529 %     saveas( gcf, num2str(fignum), 'epsc2' );
 530 %     fignum = fignum + 1;
 531 %     close;
 532 %   end
 533 % end
 534
 535 % save 'matlab_ci.prn' ci -ascii;
 536
 537
 538
 539
 540
 541