%% 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;

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

%% fit bootstrap samples
x_all = [x, x_boot];
p_est = zeros(2, n_boot);
for i = 1:n_boot
    p_est(:, i) = polyfit( t, x_all(:, i), 1 ); % the actual objective function would be minimizing the
                        % sum of squares with x_true, but this is hard to incorporate
                        % in this example since my actual application isn't a fit per se,
                        % but a model that sort of de-noises x (or x_boot).
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')