%% clear workspace
close all
clear all
clc

%% parameters
n_boot = 1e3;
n_data = 20;
sigma = 0.5;
p0 = 5;
p1 = 2;

%% generate data
t = (1:n_data)';
x_true = p0 + p1 * t;

%% generate noise
x_noise =  sigma * randn(n_data,1);
x = x_true + x_noise;

%% fit and draw bootstrap samples#
p_est = zeros(2, n_boot);
for i = 1:n_boot
    p_est(:, i) = polyfit( t, x, 1); % estimated parameters
    x_best_fit = polyval( p_est(:, 1), t ); % fitted line
    x_residuals = x_best_fit - x; % residuals of i-th fit
    x_res_boot = x_residuals( randi(n_data, n_data, 1) ); % resample residuals with replacement
    x = x_best_fit + x_res_boot; % add resampled residuals to i-th fitted line
end

%% plot
figure;
subplot(2,2,1:2)
hold on
for i = 1:n_boot
    plot( t, polyval( p_est(:, i), t ), 'r-')
end
plot(t, x, 'k.')
hold off
title('original noisy data and bootstrapped fit')
subplot(2,2,3)
hist( p_est(2, :) )
title('p0')
subplot(2,2,4)
hist( p_est(1, :) )
title('p1')