% 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=30; 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))))]) disp(['NumPoints = ', num2str(NumPoints) ' SD Noise = ', num2str(Noise) ' x-range = ' num2str(maxx-minx)]) disp(' Simulation Algebraic equation Bootstrap method') disp(' SDslope SDint SDslope SDint SDslope SDint') for test=1:10, for repeat=1:100 x=minx:(maxx-minx)/NumPoints:maxx; % Choose one of the noise models below NoiseArray=Noise*(randn(size(x))); % normal white constant noise % NoiseArray=Noise*bluenoise(x); % Blue noise % NoiseArray=Noise*pinknoise(length(x)); % Pink 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 ncutoff, 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([FullSTD PredSTD SubsetSTD]); 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) end disp('--------------------------------------------------------------')