# What is the best function to fit onto a "flat top gaussian"?

I am studying some signal and I am trying to make an automated algorithm that can extract the parameters for my signal. First let me describe a bit, I have a light emiting object crossing a window of observation with a given speed. That object is an elipse, thus when I measure the signal I find something that looks like a flat top gaussian. (see scheme)

I have been using a Gaussian Mixture Model to describe it, but because the top is flatten, it finds many more distribution than what is actually there ! Thus I am looking for a mathematical function that could describe that kind of signal, so that I can implement a least square regression or something along those line. Idealy that function would have a "variance" and a "mean" parameter such that I can make analysis about those two later on. Thanks a lot !

• There is the "Frankenstein" distribution that is known as the Johnson S_U distribution. You can probably fit it so it looks like a a flat-top Gaussian but really... why? Nov 14 '15 at 7:15
• The anwser to why, is research XD I don't necessarily choose what my signal is going to look like ^^
– Xqua
Dec 10 '15 at 21:37
• This has a similar look. stats.stackexchange.com/questions/57069/…
– Sycorax
Sep 29 '16 at 22:01

You have probably solved this problem, since this post is from a year ago. But just in case: is there a reason why you can't directly fit a flat-top Gaussian functional representation to your distribution? The general Gauss function is $g(x) = \exp[-(x/a)^n]$, "$a$" being the half-width at $1/e$. $n=2$ is a standard Gaussian, $n>2$ a flat-top Gaussian, and $n<2$ a peaked Gaussian.

• This is what I ended up doing ! ;)
– Xqua
Oct 3 '16 at 19:36

Alright so I came up with a way to do it, it's kindof a hack: Take the left part of a gaussian, plus a plateau, then the right side of the gaussian.

def multinomial_normal_fusion(x, params):
# params here is a list of [mu_0, sigma_0, amplitude_0, spread_0, mu_1, ...]
par = np.reshape(params, (len(params) / 4, 4))
n = len(par)
curve = np.zeros(len(x))
for i in range(n):
mu, sigma, k, spread = par[i]
xi = np.searchsorted(x, mu)
curve[0:xi] += k * scipy.stats.norm.pdf(x[0:xi], loc=mu, scale=sigma)
xj = np.searchsorted(x, mu + spread)
curve[xi:xj] += k * scipy.stats.norm.pdf([mu] * len(curve[xi:xj]), loc=mu, scale=sigma)
curve[xj:] += k * scipy.stats.norm.pdf(x[xj:], loc=mu + spread, scale=sigma)
return curve


If anyone has a better idea ? Thanks

• Yes that's a bad idea. A much better idea would be what you would obtain if you convolved a Gaussian function with a rectangle, which would be something like (erf(x-p1) - erf(x-p2))/2, p1 and p2 being the positions of the edge of the 1D rectangle. Aug 17 '19 at 8:02