matlab/profiles.m

   1 clear;
   2
   3 % ---------- Parameters for the plots ----------
   4
   5 rows    = 3;                            % Number of rows in each figure
   6 columns = 2;                            % Number of columns
   7 polynomial = 3;                         % Degree of polynomial for
   8                                         % fitting Y = f(ln(X))
   9 prediction_steps = 100;                 % The number of prediction
  10                                         % steps in the plotted interval
  11 extrapolate = 0.1;                      % Extrapolate this much
  12                                         % outside the available
  13                                         % data points
  14 text_offset = 0.2;                      % The text below each graph
  15                                         % will be placed this far
  16                                         % below the X-axis.
  17
  18 % ---------- Autogenerated code below ----------
  19
  20 str_format = '%30s';
  21
  22 col_names = [ sprintf(str_format,'par1');
  23               sprintf(str_format,'par2');
  24               sprintf(str_format,'par3');
  25               sprintf(str_format,'par4');
  26               sprintf(str_format,'par5');
  27               sprintf(str_format,'par6');
  28               sprintf(str_format,'par7');
  29               sprintf(str_format,'par8');
  30               sprintf(str_format,'par9');
  31               sprintf(str_format,'par10');
  32               sprintf(str_format,'par11');
  33               sprintf(str_format,'par12');
  34               sprintf(str_format,'par13')];
  35
  36 goal = 3.84;
  37
  38 filename = 'llp_rr.csv';
  39
  40 % ---------- End autogenerated code ----------
  41
  42 % The rest of the code does not change between runs
  43
  44 profiling_log = load(filename);
  45
  46 [a, b]           = size( profiling_log );
  47
  48 orig_est         = profiling_log(1,1:2:b);
  49 orig_low         = profiling_log(2,1:2:b);
  50 orig_upp         = profiling_log(3,1:2:b);
  51
  52 [junk, n_par]    = size( orig_est );
  53
  54 delta_ofv        = profiling_log(:,2:2:b);
  55 fixed_values     = profiling_log(:,1:2:b);
  56 fixed_lower      = profiling_log(2:2:a,1:2:b);
  57 fixed_upper      = profiling_log(3:2:a,1:2:b);
  58 log_fixed_values = log( fixed_values );
  59
  60
  61 % Loop over all parameters, creating plots
  62 row = 1;
  63 col = 1;
  64 for i = 1:n_par,
  65   % Plot and subplot handling
  66   if col == (columns+1)
  67     col = 1;
  68     row = row + 1;
  69   end
  70   if row == (rows+1) || i == 1
  71     h  = figure;
  72     col = 1;
  73     row = 1;
  74   end
  75
  76   % prepare fit
  77   if min(fixed_values(:,i)) <= 0 & max(fixed_values(:,i)) >= 0
  78     X = fixed_values(:,i);
  79   else
  80     X = log_fixed_values(:,i);
  81   end
  82   x_nolog = fixed_values(:,i);
  83   Y = delta_ofv(:,i);
  84   idx = isnan(X) | isnan(Y);            % remove NaN's
  85   X = X(~idx);
  86   x_nolog = x_nolog(~idx);
  87   Y = Y(~idx);
  88
  89   % Fit the polynomial
  90   P = polyfit( X, Y, polynomial);
  91
  92   % Prepare the prediction
  93   x_max = max( x_nolog );
  94   x_min = min( x_nolog );
  95   extra = ( x_max - x_min ) * extrapolate;
  96   X_pred_max = x_max + extra;
  97   X_pred_min = x_min - extra;
  98   X_pred = X_pred_min:(X_pred_max-X_pred_min)/prediction_steps:X_pred_max;
  99
 100   % Predict profile
 101   if min(fixed_values(:,i)) <= 0 & max(fixed_values(:,i)) >= 0
 102     pred_Y = polyval( P, X_pred );
 103   else
 104     pred_Y = polyval( P, log(X_pred) );
 105   end
 106
 107   % Prepare plot
 108   nplot = col+(row-1)*columns;
 109   subplot(rows, columns, nplot);
 110
 111   % Plot the gathered profile
 112   plot( x_nolog, Y, 'o' );
 113   hold on;
 114
 115   % Plot the predicted profile
 116   plot( X_pred, pred_Y );
 117
 118   ax = axis;
 119   % Add line for maximum likelihood estimate
 120   h = line( [orig_est(i), orig_est(i)], ax(3:4) );
 121
 122   % Add lines for the lower and upper confidence limits based on
 123   % the standard errors
 124   h = line( [orig_low(i), orig_low(i)], ax(3:4) );
 125   set(h, 'LineStyle', '--' );
 126   h = line( [orig_upp(i), orig_upp(i)], ax(3:4) );
 127   set(h, 'LineStyle', '--' );
 128
 129   % Add lines for the lower and upper confidence limits as the are
 130   % determined by the log-likelihood profiling tool
 131   lower = fixed_lower(:,i);
 132   idx = isnan( lower );
 133   lower = lower(~idx);
 134   lower = lower(length(lower));
 135   h = line( [lower, lower], ax(3:4) );
 136   upper = fixed_upper(:,i);
 137   idx = isnan( upper );
 138   upper = upper(~idx);
 139   upper = upper(length(upper));
 140   h =line([upper, upper], ax(3:4));
 141
 142   % Add a reference line for the ofv increase
 143   h = line( ax(1:2),[goal, goal]);
 144
 145   % ----------- Add some text and figures ----------
 146
 147   % Maximum likelihood estimate
 148   y_text = ax(3) - (ax(4)-ax(3))*text_offset;
 149   text( orig_est(i), y_text, num2str(orig_est(i)) );
 150
 151   % upper SE confidence limit
 152   y_text = ax(3) - (ax(4)-ax(3))*text_offset*2/3;
 153   text( orig_upp(i), y_text, num2str(orig_upp(i)) );
 154
 155   % lower SE confidence limit
 156   text( orig_low(i), y_text, num2str(orig_low(i)) );
 157
 158   % Upper likelihood-profiling confidence limit
 159   y_text = ax(3) - (ax(4)-ax(3))*text_offset;
 160   text( upper, y_text, num2str(upper) );
 161
 162   % Lower likelihood-profiling confidence limit
 163   text( lower, y_text, num2str(lower) );
 164
 165   hold off;
 166
 167   % Add title to each graph
 168   title( col_names( i, : ) );
 169
 170   col = col + 1;
 171 end
 172
 173
 174