% Script that compares standard deviation of slope and intercept for a
% first-order least-squares fit computed by random-number simulation to
% predictions made by closed-form algebraic equations and to a bootstrap.
NumPoints=100; % Number of data points in signal
Noise=10; % Standard deviation of the random noise
TrueSlope=2;
TrueIntercept=10;
minx=1;
maxx=50;
SmoothWidth=0;

FullResults=zeros(100,2);
% Prediction using empirical equation derived from Monte Carlo simulations
% disp(['3*log(pi)*Noise./((maxx-minx)*(sqrt(NumPoints))) = '
% num2str(3*log(pi)*Noise./((maxx-minx)*(sqrt(NumPoints))))])
for repeat=1:100
    x=minx:(maxx-minx)/NumPoints:maxx;
    NoiseArray=Noise*(randn(size(x))); % normal white constant noise
    % NoiseArray=bimodal(x,Noise,30,-30); % bimodal white constant noise
    % NoiseArray=(TrueIntercept+x*TrueSlope).*(randn(size(x)))./10; % normal white proportional noise
    y=TrueIntercept+x*TrueSlope+NoiseArray;
    if SmoothWidth,y=fastsmooth(y,SmoothWidth,3,1);end
    PredSet=polyfit(x,y,1);
    if repeat==1,
        FullResults=PredSet;
    else
        FullResults=[FullResults; PredSet];
    end
end % repeat
FullMean=mean(FullResults);
FullSTD=std(FullResults);
% Least-squares computation without polyfit function.
Sxx=sum((x-mean(x)).^2);
Syy=sum((y-mean(y)).^2);
Sxy=sum((x-mean(x)).*(y-mean(y)));
CalcSlope=Sxy./Sxx;
CalcIntercept=mean(y)-CalcSlope*mean(x);
Sy=sqrt((Syy-CalcSlope^2*Sxx)/(NumPoints-2));
SDSlope=Sy/sqrt(Sxx);
SDIntercept=Sy*sqrt(1./(NumPoints-(sum(x).^2)./sum(x.^2)));
PredMean=[CalcSlope CalcIntercept];
PredSTD=[SDSlope SDIntercept];
BootstrapResultsMatrix=zeros(100,2);
clear BootstrapResultsMatrix xx yy
cutoff=0.5;
for trial=1:500,
    n=1;
    xx=x;
    yy=y;
    while n<length(x)-1,
        if rand>cutoff,
            xx(n)=x(n+1);
            yy(n)=y(n+1);
        end
        n=n+1;
    end
    PolyFitResults=polyfit(xx,yy,1);
    BootstrapResultsMatrix(trial,:)=[PolyFitResults(1) PolyFitResults(2)]';
end
%subplot(2,2,1)
SubsetMean=mean(real(BootstrapResultsMatrix));
SubsetSTD=sqrt(2)*std(BootstrapResultsMatrix);
if SmoothWidth,disp(['Smooth Width = ', num2str(SmoothWidth)]);end
disp(['NumPoints = ', num2str(NumPoints) '     SD Noise = ', num2str(std(NoiseArray)) '    x-range = ' num2str(maxx-minx)])
disp('    Simulation        Algebraic equation    Bootstrap method')
disp('    SDslope  SDint      SDslope  SDint      SDslope  SDint')
disp([FullSTD PredSTD SubsetSTD]);
disp('--------------------------------------------------------------')
figure(1)
plot(x,y,'o',linspace(min(x),max(x)),TrueSlope*linspace(min(x),max(x))+TrueIntercept)
figure(2)
plot(xx,yy,'o',linspace(min(xx),max(xx)),TrueSlope*linspace(min(xx),max(xx))+TrueIntercept)