% ch13_bootstrap_fit.m
% Companion code for Chapter 13: Optimization and Data Fitting.
%
% Estimate straight-line parameter uncertainties by pairs bootstrap.
% Resample complete (x,y) pairs with replacement and refit each sample.
% Change the data and number of bootstrap samples below.
% This example treats observed pairs as independent sampling units. For
% fixed experimental x values with a specified noise model, consider a
% residual or parametric bootstrap instead. Correlated measurements require
% a resampling scheme that preserves their dependence (e.g., blocks).
% Bootstrap standard errors are not confidence intervals and do not account
% for systematic errors or an incorrect fitting model. No toolbox is needed.

clear; close all; clc;

%% Example data and settings
seed = 1;
rng(seed, 'twister');
% Alternative: comment out the line above and uncomment the following.
% rng('shuffle');

x = linspace(0,1,25)';
y = 2.0 + 0.5*x + 0.15*randn(size(x));
B = 2000;                 % Number of resampled data sets
make_plot = true;

% For measured data, replace the example above, for instance:
% data = readmatrix('measurements.txt'); x = data(:,1); y = data(:,2);
% Keep paired observations together. This program uses unweighted fits.

%% Original fit and bootstrap
[estimate, samples, info] = bootstrap_line(x,y,B);
% Parameter order throughout: [intercept; slope].
fprintf('\nPairs bootstrap for y = a0 + a1*x\n');
fprintf('Measurements:      %d\n',numel(x));
fprintf('Bootstrap samples: %d\n',B);
fprintf('%12s %16s %16s %16s\n','Parameter','Original fit','Bootstrap SE','Estimated bias');
fprintf('%12s %16.8g %16.8g %16.8g\n','Intercept', ...
    estimate(1),info.standard_errors(1),info.bias(1));
fprintf('%12s %16.8g %16.8g %16.8g\n','Slope', ...
    estimate(2),info.standard_errors(2),info.bias(2));
fprintf('SE is the sample standard deviation of bootstrap estimates.\n');
fprintf('Bias is the bootstrap mean minus the original estimate.\n');
fprintf('Increase B to check stability of the reported uncertainties.\n');
% Full bootstrap covariance matrix is available as info.covariance.
% samples(:,b) contains the parameters of bootstrap fit b.

%% Optional parameter histograms
if make_plot
    figure('Color','w');
    labels = {'Intercept a_0','Slope a_1'};
    for j=1:2
        subplot(2,1,j);
        hh = histogram(samples(j,:),30,'FaceColor',[0.85 0.85 0.85], ...
            'EdgeColor','k');
        hold on;
        limits = ylim;
        hl = plot([estimate(j) estimate(j)],limits,'k-','LineWidth',1.5);
        xlabel(labels{j},'FontSize',14); ylabel('Count','FontSize',14);
        title('Pairs-bootstrap parameter distribution', ...
            'FontSize',14,'FontWeight','normal');
        legend([hh,hl],{'Bootstrap estimates','Original estimate'},'Location','best');
        set(gca,'FontSize',14); grid on; box on;
    end
end

%% Reusable numerical routine
function [estimate,samples,info] = bootstrap_line(x,y,B)
% Return the original [intercept; slope] and a 2-by-B array of estimates.
% Copy this routine and its helper to bootstrap_line.m for reuse.
    validateattributes(x,{'double'},{'real','vector','finite','nonempty'});
    validateattributes(y,{'double'},{'real','vector','finite','nonempty'});
    validateattributes(B,{'double'},{'real','scalar','integer','finite','>=',2});
    x=x(:); y=y(:); n=numel(x);
    if numel(y)~=n || n<3
        error('Supply matching x and y vectors with at least three points.');
    end
    estimate=line_coefficients(x,y);
    samples=zeros(2,B);
    for b=1:B
        idx=randi(n,n,1);
        % The same indices select both coordinates of each observed pair.
        samples(:,b)=line_coefficients(x(idx),y(idx));
    end
    info=struct('standard_errors',std(samples,0,2), ...
        'bias',mean(samples,2)-estimate,'covariance',cov(samples.'));
end

function coefficients=line_coefficients(x,y)
    xbar=mean(x); ybar=mean(y);
    dx=x-xbar; Sxx=sum(dx.^2);
    if Sxx<=0 || ~isfinite(Sxx)
        error(['A data set has no usable spread in x, so its slope is undefined. ' ...
            'Pairs bootstrap may be unsuitable for this small or repeated-x data set.']);
    end
    slope=sum(dx.*(y-ybar))/Sxx;
    coefficients=[ybar-slope*xbar; slope];
end