%% Nonlinear regression. Generate data for a dose response curve by adding
%% some normally distributed noise to a Hill function. Use nlinfit to
%% perform nonlinear regression

%% Specify model parameters. The Hill function is given by 
%% y(x) = ymax/(1 + (EC50/dose)^n)
ymax = 1;
Lec50 = 0;
n = 2;
modelpar = [ymax Lec50 n];

%% Generate some data
dose = (logspace(-1, 1, 15))';
yTrue = Hill(modelpar, dose);

%% Add normal error to simulate experimental noise
NoiseStd = 0.1; % The standard deviation of the noise
err = NoiseStd*randn(length(yTrue), 1);
yExpt = yTrue + err;

%% Call nlinfit to perform nonlinear regression. Add an option to see the
%% SSR value after each iteration. Als
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)
hold on
doseFit = (logspace(-1, 1, 201))';
yFit = Hill(betaHat, doseFit);
plot(doseFit, yFit, 'b-', 'LineWidth', 2)
plot(doseFit, Hill(modelpar, doseFit), 'b--', 'LineWidth', 1)
xlim([0.09, 11])
ylim([-0.1, 1.15])
xlabel('dose')
ylabel('response')
grid on
legend('Data', 'Nonlinear fit', 'True model', 'Location', 'NW')