function [MeasuredSTD PredictedSTD ResultsMatrix]=GaussFitMC(TrueHeight,TruePosition,TrueWidth,NumPoints,Noise,NumTrials)
% Monte Carlo simulation of the measurement of the peak height, position,
% and width of a noisy x,y Gaussian peak, by fitting a parabola 
% to a region of x vs ln(y) equal to one Gaussian half-width, 
% using the gaussfit function (included as a subfunction). Plots a 20-bin
% histogran of the NumTrials results of peak height, position,
% and width, plus a randomly selected single data set and Gaussian fit.
% The input arguments are:
% TrueHeight = True peak height of the Gaussian peak.
% TruePosition = True x-axis value at the peak maximum.
% TrueWidth =  True half-width (FWHM) of the Gaussian peak.
% NumPoints = Number of points taken for the least-squares fit
% Noise = standard deviaiton of (normally-distributed) random noise
% NumTrials = Number of trials used to compute standard deviation
% Optional Outputs:
% MeasuredSTD = computed standard deviation of [Height Position Width]
% PredictedSTD = Predicted standard deviation (SD) based on empirical
% equations based on a similar Monte Carlo simulation: 
%    SD of the peak height = 1.73*noise/sqrt(NumPoints),
%    SD of the peak position = noise/sqrt(NumPoints,) 
%    SD of the peak width = 3.62*noise/sqrt(NumPoints). 
% ResultsMatrix = Measured [Height Position Width] of each trial
% Example: GaussFitMC(100,100,100,50,1,1000)
% ans =
%     0.2454    0.1365    0.5182
startx=TruePosition-TrueWidth/2;
endx=TruePosition+TrueWidth/2;
deltax=(endx-startx)/(NumPoints-1);
x=startx:deltax:endx;
ResultsMatrix=zeros(3,NumTrials);
for trial=1:NumTrials,
    NoiseArray=Noise*randn(size(x));
    y=TrueHeight*gaussian([startx:deltax:endx],TruePosition,TrueWidth)+NoiseArray;
    [Height, Position, Width]=gaussfit(x,y);
    ResultsMatrix(:,trial)=[Height, Position, Width];
end
FullDataSetSTD=std(ResultsMatrix');
RootFactor=sqrt(length(x));
HeightFactor=1.73;
PositionFactor=.99;
WidthFactor=3.62;
% Display results
TypicalHeightPositionWidth=[Height, Position, Width];
MeasuredSTD=FullDataSetSTD;
PredictedSTD=[Noise.*HeightFactor./RootFactor Noise.*PositionFactor./RootFactor Noise.*WidthFactor./RootFactor];
subplot(2,2,1)
hist(ResultsMatrix(1,:),20)
title('Histogram of Peak Heights');
xlabel('Height')
ylabel('Probability')
subplot(2,2,2)
hist(ResultsMatrix(2,:),20)
title('Histogram of Peak Positions');
xlabel('Position')
ylabel('Probability')
subplot(2,2,3)
hist(ResultsMatrix(3,:),20)
title('Histogram of Peak Widths');
xlabel('Width')
ylabel('Probability')
subplot(2,2,4)
plot(x,y,'o',linspace(min(x),max(x)),Height*gaussian(linspace(min(x),max(x)),Position,Width))
title('Single data set');
xlabel('x')
ylabel('y')

function [Height, Position, Width]=gaussfit(x,y)
% Converts y-axis to a log scale, fits a parabola
% (quadratic) to the (x,ln(y)) data, then calculates
% the position, width, and height of the
% Gaussian from the three coefficients of the
% quadratic fit.  This works only if the data have
% no baseline offset (that is, trends to zero far off the
% peak) and if there are no zeros or negative values in y.
% Example: [Height, Position, Width]=gaussfit([1 2 3],[1 2 1]) 
% returns Height = 2, Position = 2, Width = 2
z=log(y);
coef=polyfit(x,z,2);
a=coef(3);
b=coef(2);
c=coef(1);
Height=exp(a-c*(b/(2*c))^2);
Position=-b/(2*c);
Width=2.35482/(sqrt(2)*sqrt(-c));

function g = gaussian(x,pos,wid)
%  gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  T. C. O'Haver, 1988
% Examples: gaussian([0 1 2],1,2) gives result [0.5000    1.0000    0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.6005615.*wid)) .^2);