%% clear workspace
close all
clear all
clc

%% parameters
n_boot = 1e2;
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 original data
p_est = zeros(2, n_boot);
p_est(:, 1) = polyfit( t, x, 1);
x_residuals = polyval( p_est(:, 1), t ) - x;

%% draw bootstrap samples
x_res_boot = x_residuals( randi( n_data, n_data, n_boot ) );
x_boot = bsxfun( @plus, x_res_boot, x_true );

%% fit bootstrap samples
for i = 2:n_boot
    p_est(:, i) = polyfit( t, x_boot(:, i), 1 );
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')