%  Figure 14.1 - Plots actual data, a logistic regression, and a
%  probit regression, in order to show that (a) the fits are good, and
%  (b) the logistic and probit models are very similar.
%
%  Figure caption: Two curves fitted to the data in Figure 8.9.  The
%  fitted curve from probit regression (dashed line) is shown together
%  with the fitted curve from logistic regression.  The fits are very
%  close to eachother.

%  This data comes from the "HSP" data set in R's MPDiR library.
%  These are the observations for subject ("Obs") "SS" and run
%  ("Run") "R1".
lightIntensity = [1.382017 1.575188 1.767898 1.959041 2.151982 2.344981];
percentDetected = [0.00 0.04 0.18 0.54 0.94 1.00];
nObservations = numel(percentDetected);

%  glmfit requres a row of ones.
covariates = [percentDetected', ones(nObservations)];

%  Fit models.
logisticModel = glmfit(lightIntensity, covariates, 'binomial', 'link', 'logit');
probitModel = glmfit(lightIntensity, covariates, 'binomial', 'link', 'probit');



%  Plot observed data.
h(1) = plot(lightIntensity, percentDetected, '.k', 'Markersize', 16);

%  Set axes and labels.
set(gca,'Box','off', 'FontSize', 18, 'TickDir', 'out', ...
        'XLim', [1, 2.5], 'YLim', [0, 1], ...
        'XTick', 1:0.5:2.5, 'YTick', 0:0.2:1)

xlabel('Log Intensity', 'FontSize', 18)
ylabel('Proportion Perceived', 'FontSize', 18)
hold on;

%  Get models' predictions and plot them.
intensityValues = 1:0.01:2.5;
logisticPredictions <- ...
    1 ./ (1 + exp( - (logisticModel(1) + intensityValues * logisticModel(2))))
probitPredictions   <- ...
    normcdf(probitModel(1) + intensityValues * probitModel(2), 0, 1)

h(2) = plot(intensityValues, logisticPredictions, '-k', 'LineWidth', 3);
h(3) = plot(intensityValues, probitPredictions, '--k', 'LineWidth', 3);

lh = legend(h(2:3), 'Logit', 'Probit', 'Location', 'Northwest');
set(lh, 'Box', 'off')