0

I'm trying to map a set of data points that seems to form a kind of a flat top gaussian in one dimension. I'm used to fitting data to a gaussian using non-linear least squares. I'm trying to understand how to do the same, for a flat top gaussian?

dev_nut
  • 161
  • 3
  • 11

2 Answers2

2

The approach is the same. A Gaussian distribution has two parameters, the mean and standard deviation. You add up the square of the error of each point to get an error function, which depends on those two parameters, then minimize the error function to find the best fit set of parameters. For your case, you need to define the "flat topped Gaussian" as a function of some number of parameters. You might just chop off the head of the Gaussian at $\pm n \sigma$, leading to something like $$M(\mu , \sigma , n,x)=\begin {cases} \frac K{\sqrt {2 \pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma ^2}}&|x| \gt n \sigma \\ \frac K{\sqrt {2 \pi \sigma^2}} e^{-\frac{(n \sigma-\mu)^2}{2\sigma ^2}}&|x| \le n \sigma\end {cases}$$ where you choose $K$ as a function of $n$ so the integral over the distribution remains $1$ Now compute the error of each point from $M$, square it, add them all up, and that is your error function to minimize over the parameters.

Ross Millikan
  • 374,822
  • is it possible for your to elaborate why the choice of $n\sigma$, acts as a sort of step function in this case to flatten out the top? – dev_nut Nov 05 '15 at 18:46
  • In the normal Gaussian, when $x = \mu$ the exponential factor is $1$. I have forced the exponential factor to be less than $1$ when $x$ is near the mean. This form will have corners at $\pm n \sigma$, which you may not like. I suggest you plot this function for $\mu=0, \sigma=1,$ and $n=0.2, 0.5, 1, 2$ (don't worry about $K$-you can set it to $1$ for these plots) and see if you like the result – Ross Millikan Nov 05 '15 at 18:50
  • The choice of a fitting function can be art, not science, unless you have a theoretical basis for the function. You want relatively few parameters. The corners I mentioned before may look strange and tempt you to add another parameter to be the radius of a circle that rounds them off. I don't think it will matter much, but with more parameters you can get funny fits. – Ross Millikan Nov 05 '15 at 18:53
  • I like your formulation. I think it's simple and great. I just don't have the mathematical intuitiveness to understand how it works. I will graph it out and find out graphically. – dev_nut Nov 05 '15 at 19:55
  • Should the first term contain a negative sign in the exponential, to form a gaussian like? – dev_nut Nov 05 '15 at 20:08
  • The exponential of both terms should be negated. Then, it forms the gaussian, otherwise an inverted gaussian. – dev_nut Nov 05 '15 at 20:42
  • You are correct, it should be negated. I have fixed it. – Ross Millikan Nov 05 '15 at 21:20
  • The $n$ in the equation, does not mean the number of points of the data set right? It's just another variable to be determined. I'm asking just to clear out any ambiguity. – dev_nut Nov 05 '15 at 23:25
  • Yes, it is the number of $\sigma$ that your flat top extends each side of the mean. If $n=0$ you get a standard Gaussian. If $n=1$ the flat top covers $\mu \pm 1\sigma$ – Ross Millikan Nov 06 '15 at 00:04
-1

This is my implementation in Matlab. I use FWHM (full width half maximum) which is the location where the Gaussian is equal to 0.5. FunctionResult

function[y] = flatTopGaussian(sizex,FWHM,plot_result)
% Produces a flat top Gaussian.
%   Equal to 1 between -FWHM/2,FWHM/2
%   Gaussian function for the other values
% Inputs
%   sizex:  length of y
%   FWHM:   Full width half maximum
%   plot_result: should we plot the result

%%
% sizex=200;
% FWHM=2*20;
x_vec=-sizex/2+1:sizex/2;
sigma=FWHM/2.355;
[y] = 2*createGaussianPeak_1D(sizex, sigma);
y(abs(x_vec)<FWHM/2)=1;

%% Plotting
if plot_result
    flat_win=zeros(size(x_vec)); flat_win(abs(x_vec)<FWHM/2)=1;
    figure(1); plot(x_vec,flat_win,x_vec,y);
end
end

function [y] = createGaussianPeak_1D(sizex,output_sigma)
sxh = ceil((sizex) / 2);
mult = -0.5 / (output_sigma * output_sigma);
y=zeros(1,sizex);
for j=1:sizex
    jh = j - sxh;
    y(j) = exp(mult * (jh * jh));
end
end
lahmania
  • 101