gaussian_fit.py

from scipy import optimize
from scipy import stats


def gaussian(x, p):
    n = len(p) // 2
    return np.add.reduce([p[i] * stats.norm.pdf(x, scale=p[n + i]) for i in range(n)])


def fit_gauss_sym(x, y, n=2, origin=True, full=True, positive=True):
    """
    x, y:
        The right half (x>0) of an even function.
        x and y are monotonically increasing and decreasing respectively.
    """
    if origin:
        y1 = y[0] * np.linspace(1, 0, 2 * n + 1)[:-1]
    else:
        y1 = y[0] * np.linspace(1, 0, 2 * n + 2)[1:-1]
    x1 = np.interp(-y1, -y, x)

    if positive or n == 1:
        s0 = x1[n] / 1.18 * np.geomspace(0.8, 1.25, n + 2)[1:-1]
        a0 = y1[0] / stats.norm.pdf(0, scale=s0) / n
        p0 = [*a0, *s0]
    else:
        s0 = x1[n] / 1.18 * np.geomspace(0.8, 1.25, n + 2)[1:-1]
        a0 = y1[0] / stats.norm.pdf(0, scale=s0) / (n - 1)
        a0[0] *= 1.05
        a0[-1] *= -0.05
        p0 = [*a0, *s0]

    if positive:
        lnp0 = np.log(p0)
    else:
        lnp0 = *np.arcsinh(p0[:n]), *np.log(p0[n:])

    if full:
        args = (x, y)
    else:
        args = (x1, y1)

    def func_obj(lnp, x, y):
        if positive:
            p = np.exp(lnp)
        else:
            p = *np.sinh(lnp[:n]), *np.exp(lnp[n:])
        f = gaussian(x, p)
        return f - y

    with np.errstate(all='ignore'):
        res = optimize.root(func_obj, lnp0, args=args, method='lm', options={'ftol': 1e-4 * y1[0], 'maxiter': 2000})
    lnp = res.x

    if positive:
        p = np.exp(lnp)
    else:
        p = *np.sinh(lnp[:n]), *np.exp(lnp[n:])

    return p, p0, x1, y1, res
  
  
x = np.linspace(0, 20, 1000)
y = 0.2 * norm.pdf(x) + norm.pdf(x, scale=4)
# x = np.linspace(0, np.pi / 2, 100)
# y = np.cos(x)**2

n = 3
p, p0, x1, y1, res = fit_gauss_sym(x, y, n=n, origin=1, full=1, positive=0)
printarr(p0)
printarr(p)
printarr(res.success, (res.fun**2).mean()**0.5)

plt.plot(x, y)
plt.scatter(x1, y1)
ylim = plt.ylim()

plt.plot(x, func(x, p), ls='--')
plt.plot(x, func(x, p0), ls='-.')
for i in range(n):
    plt.plot(x, func(x, [p[i], p[n + i]]), ls=':')
plt.ylim(0, ylim[1])

gaussian_fit_simple.py

    x = np.hstack([-t[::-1], t])
    y = np.hstack([A[::-1], A])
    y = y - y[-1]

    f = CubicSpline(x, y, extrapolate=False)
    amp = f(0)

    sig1 = (-amp / f(0, nu=2))**0.5  # |d^2f/dt^2|=-A/sigma^2

    rt = f.derivative(2).roots(extrapolate=False)  # root of d^2f/dt^2 = 1 sigma
    sig2 = np.abs(rt[np.argmin(np.abs(rt))])  # root closest to 0, usually smaller than sig1?

    sig = (sig1 * sig2)**0.5  # average (purely empirical!)

    amp *= (2 * np.pi)**0.5 * sig