%% Nonlinear regression using nlinfit (Statistics toolbox)

%% Fit simulated dose-response data to a Hill function, given by 
%% y(x) = ymax/(1 + (EC50/dose)^n)

%% Load data
DoseResponse = load('DoseResponseData.dat', '-ascii');
dose = DoseResponse(:,1);
yExpt = DoseResponse(:,2);

%% Call nlinfit to perform nonlinear regression. Add an option to view 
%% SSR after each iteration
options = statset('Display', 'iter');
betaGuess = [1, 1, 1] % Initial guess for parameter values
betaHat = nlinfit(dose, yExpt, @Hill, betaGuess, options)

%% Not all initial guesses work. This one fails to find the minimum. 
%betaGuess = [0.1, 1, 4]
%nlinfit(dose, yExpt, @Hill, betaGuess, options)

%% Plot data and fit
semilogx(dose, yExpt, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r')
hold on
doseFit = (logspace(-1, 1, 201))';
yFit = Hill(betaHat, doseFit);
plot(doseFit, yFit, 'b-', 'LineWidth', 2)
modelpar = [1, 0, 2]; % Make sure that these match the values in 'NonlinearRegressionMakeDat.m'
plot(doseFit, Hill(modelpar, doseFit), 'b--', 'LineWidth', 1)
xlim([0.09, 11])
ylim([-0.05, 1.15])
xlabel('dose')
ylabel('response')
grid on
legend('Data', 'Nonlinear fit', 'True model', 'Location', 'NW')

%% Compute and plot residuals for the best fit
yMod = Hill(betaHat, dose);
residuals = yExpt - yMod;

figure()
semilogx(dose, residuals, 'bo--', 'MarkerSize', 7, 'MarkerFaceColor', 'b')
xlabel('dose')
ylabel('residuals')
grid on