% Script that compares standard deviation of slope and intercept % for a first-order least-squares fit computed by Monte Carlo % simulation of 1000 repeats to predictions made by closed-form % algebraic equations. Change the values in lines 5-10 to suit. % Version 2, April 2015, changed STD function to be compatible with Octave. % Noise=5; % Standard deviation of the noise in the data TrueSlope=1; % Actual underlying slope of the data TrueIntercept=0; % Actual underlying y-intercept of the data minx=1; % minimum value of x maxx=100; % maximum value of x NumPoints=100; % Number of points in the data % % Monte Carlo simulation FullResults=zeros(10000,2); for repeat=1:1000 x=minx:(maxx-minx)/NumPoints:maxx; NoiseArray=Noise*(randn(size(x))); % Constant noise added to each data point y=TrueIntercept+x*TrueSlope+NoiseArray; PredSet=polyfit(x,y,1); if repeat==1, FullResults=PredSet; else FullResults=[FullResults; PredSet]; end end % repeat FullMean=mean(FullResults); if IsOctave, FullSTD=[stdev(FullResults(:,1)) stdev(FullResults(:,2))]; else FullSTD=std(FullResults); end % % Algebraic least-squares computation. 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]; % % Results comparison disp(' ') if IsOctave, disp(['NumPoints = ', num2str(NumPoints) ' SD Noise = ', num2str(stdev(NoiseArray)) ' x-range = ' num2str(maxx-minx)]) else disp(['NumPoints = ', num2str(NumPoints) ' SD Noise = ', num2str(std(NoiseArray)) ' x-range = ' num2str(maxx-minx)]) end disp(' MC Simulation Algebraic Prediction') disp(' slope intercept slope intercept') disp(['Mean: ', num2str([FullMean PredMean]) ]) disp(['Standard deviation: ', num2str([FullSTD PredSTD]) ]) disp('--------------------------------------------------------------') plot(x,y,'o',linspace(min(x),max(x)),TrueSlope*linspace(min(x),max(x))+TrueIntercept) xlabel('x') ylabel('y')