% eda12_12
% illustration of Bootstrap statistics of slope in least
% squares fit of a straight line to a small fragment of
% Black Rock Forest temperature data

% load data and set up time axis
dorig = load('brf_fragment.txt');
N = length(dorig);
Dt=1.0;
torig = 10*Dt+Dt*[0:N-1]';  % hours from start of orignal record

% least-squares slope and its variance of original data
M=2;
G=zeros(N,M);
G(:,1)=1;
G(:,2)=torig;
mest=(G'*G)\(G'*dorig);
slopeorig=mest(2);
    
% usual error estimate based on propagagtion
% of posterior estimate of sigma2 of data into
% a sigma2 of slope
e = dorig-G*mest;
E = e'*e;
sigma2dorig = E/N;
Cmorig = sigma2dorig * inv(G'*G);
sigma2morig = Cmorig(2,2);

% number of realizations;
Nr = 10000;
slope = zeros(Nr,1);

% loop over realizations
for p = [1:Nr]
    
    % resample
    rowindex=unidrnd(N,N,1);
    t = torig(rowindex);
    d = dorig(rowindex);

    % straight line fit
    M=2;
    G=zeros(N,M);
    G(:,1)=1;
    G(:,2)=t;
    mest=(G'*G)\(G'*d);
    slope(p)=mest(2);
        
end

% compute histogram of slopes and normalize
% normalize to a probability density distribution
% with Dt*sum(pbootstrap)=1;
Nbins = 100;
[shist, bins]=hist(slope, Nbins);
Db = bins(2)-bins(1);
pbootstrap = shist / (Db*sum(shist));

% compute Normal distribution of slopes
% using parameters estimated from original data
pnormal = (1/sqrt(2*pi*sigma2morig)) * exp(-((bins-slopeorig).^2)/(2*sigma2morig));

figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [bins(1), bins(Nbins), 0, 1.2*max(pbootstrap)]); 
plot(bins,pbootstrap,'k-','LineWidth',2);
plot(bins,pnormal,'k-','LineWidth',2);
xlabel('slope, b');
ylabel('p(b)');

% mean (expectation) and variance of bootstrap distibution
mb = Db*sum(bins.*pbootstrap);
vb = Db*sum(((bins-mb).^2).*pbootstrap);

% 5% and 95% confidence of bootstrap distibution
Pbootstrap = Db*cumsum(pbootstrap);
ilo = find(Pbootstrap>=0.025,1);
ihi = find(Pbootstrap>=0.975,1);
blo = bins(ilo);
bhi = bins(ihi);

disp(sprintf('standard method:   slope %f +/- %f (2 sigma)', slopeorig, 2*sqrt(sigma2morig)));
disp(sprintf('bootstrap method:  slope %f +/- %f (2 sigma)', mb, 2*sqrt(vb)));
disp(sprintf('bootstrap method:  slope between %f and %f (95 percent confidence)', blo, bhi));


