%% Construct bootstrap samples for a nonlinear regression and estimate
%% bootstrap confidence intervals. 

%% 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.
betaGuess = [1, 1, 1]; % Initial guess for parameter values
[betaHat, residuals, Jac, Cov, mse] = nlinfit(dose, yExpt, @Hill, betaGuess);

%% Resample the residuals using the bootstrp function
nboot = 200; % Number of bootstrap replicates
[bootresiduals, bootIndices] = bootstrp(nboot, [], residuals);
bootIndices(:,1:3) % Print the first 3 columns of bootIndices
bootResiduals = residuals(bootIndices);
bootResiduals(:, 1:3) % Print the first 3 columns of bootResiduals

%% Generate bootstrap datasets by adding the resampled residuals 
%%to the best fit curve
yMod = Hill(betaHat, dose);
yBoot = repmat(yMod, 1, nboot) + bootResiduals;

%% Fit each of the bootstrap datasets to the Hill function to 
%% build up the bootstrap distributions of parameter estimates
betaBoot = zeros(nboot, 3);
for i=1:nboot
	betaBoot(i,:) = nlinfit(dose, yBoot(:,i), @Hill, betaGuess);
end

% Estimate 95% CIs
bootCI = prctile(betaBoot,[2.5 97.5]);
betaHat
bootCI

%% Plot histograms of bootstrap distributions
for i = 1:3
	subplot(3,1,i)
	hist(betaBoot(:,i), 30)
	h = findobj(gca,'Type','patch');
	set(h, 'EdgeColor', 'w')
end

axes(subplot(311))
xlabel('y_{max}')
axes(subplot(312))
xlabel('log_{10}(EC_{50})')
axes(subplot(313))
xlabel('Hill coefficient, n')