Fit Data to Gauß-Function with 2 peaks

Question

Currently, I am trying to fit data to a Gauß-Function with 2 peaks. Here is my current setup for that (Data at the End):

Plot of the data I am trying to fit:

Gauß Function:

def doublegaussian(x,x0,mu1,mu2,sigma,c,y0): 
    return y0 + 1/np.sqrt(2*np.pi*sigma**2)*np.exp(-((x+x0)-mu1)**2/(2*sigma**2)) +((1/np.sqrt(2*np.pi*sigma**2)*np.exp(-((x+x0)-mu2)**2/(2*sigma**2))+c))

Where x=data, x0 = x-offset, mu1/mu2=peak offsets, sigma=width of peak, c=offset cause the right peak is higher than the left, y0 = offset because the paks start at about 1.4.

def fitdoublegaussian(array,x0,mu1,mu2,sigma,c,y0,d1,d2):
    x = array[d1:d2,0]
    y = array[d1:d2,1]
    popt1,pcov1 = curve_fit(doublegaussian,x,y,p0=[x0,mu1,mu2,sigma,c,y0])
    print(popt1)
    plt.scatter(x,y,label="data",s=2)
    plt.plot(x,doublegaussian(x,*popt1),label="fit")
    plt.legend(loc=0)
    plt.savefig("doublegauß.png")

fitdoublegaussian(bf1,730,-45,40,60,0.1,1.4,40,70)

So these are the parameters I thought would work, but they don't. Thank you for your help!

Data:

400,000000  2,828000E-12
406,000200  2,860000E-12
412,000400  2,898000E-12
418,000600  2,960000E-12
424,000800  3,017000E-12
430,001000  3,104000E-12
436,001200  3,308000E-12
442,001399  3,283000E-12
448,001599  3,480000E-12
454,001799  3,450000E-12
460,001999  3,420000E-12
466,002199  3,478000E-12
472,002399  3,746000E-12
478,002599  3,739000E-12
484,002799  3,843000E-12
490,002999  3,978000E-12
496,003199  4,073000E-12
502,003399  4,213000E-12
508,003599  4,592000E-12
514,003798  4,511000E-12
520,003998  4,651000E-12
526,004198  4,827000E-12
532,004398  4,973000E-12
538,004598  5,152000E-12
544,004798  5,395000E-12
550,004998  5,592000E-12
556,005198  5,781000E-12
562,005398  5,991000E-12
568,005598  6,277000E-12
574,005798  6,537000E-12
580,005998  6,820000E-12
586,006198  7,165000E-12
592,006397  7,600000E-12
598,006597  8,103000E-12
604,006797  8,635000E-12
610,006997  9,473000E-12
616,007197  1,007200E-11
622,007397  1,077200E-11
628,007597  1,164200E-11
634,007797  1,268300E-11
640,007997  1,414100E-11
646,008197  1,502000E-11
652,008397  1,751500E-11
658,008597  1,832500E-11
664,008796  1,963600E-11
670,008996  2,125600E-11
676,009196  2,273500E-11
682,009396  2,316600E-11
688,009596  2,341400E-11
694,009796  2,308400E-11
700,009996  2,212600E-11
706,010196  2,071000E-11
712,010396  1,944100E-11
718,010596  1,827400E-11
724,010796  1,752200E-11
730,010996  1,708400E-11
736,011196  1,696500E-11
742,011395  1,757800E-11
748,011595  1,878600E-11
754,011795  2,060600E-11
760,011995  2,236800E-11
766,012195  2,431600E-11
772,012395  2,491800E-11
778,012595  2,426500E-11
784,012795  2,282800E-11
790,012995  2,046100E-11
796,013195  1,830100E-11
802,013395  1,628500E-11
808,013595  1,462500E-11
814,013794  1,324000E-11
820,013994  1,190800E-11
826,014194  1,087500E-11
832,014394  9,980000E-12
838,014594  9,130000E-12
844,014794  8,379000E-12
850,014994  7,726000E-12
856,015194  7,110000E-12
862,015394  6,756000E-12
868,015594  6,301000E-12
874,015794  5,896000E-12
880,015994  5,535000E-12
886,016194  5,263000E-12
892,016393  4,995000E-12
898,016593  4,793000E-12
904,016793  4,571000E-12
910,016993  4,382000E-12
916,017193  4,208000E-12
922,017393  4,109000E-12
928,017593  3,962000E-12
934,017793  4,111000E-12
940,017993  3,753000E-12
946,018193  3,617000E-12
952,018393  4,186000E-12
958,018593  3,429000E-12
964,018792  3,320000E-12
970,018992  3,232000E-12
976,019192  3,120000E-12
982,019392  3,070000E-12
988,019592  3,126000E-12
994,019792  2,945000E-12
1000,019992 3,705000E-12

I tried fitting it with the double-gauß function, but it doesn't work.

Answer 1

Inappropriate model. In

y0 + 1/np.sqrt(2*np.pi*sigma**2)*np.exp(-((x+x0)-mu1)**2/(2*sigma**2)) +((1/np.sqrt(2*np.pi*sigma**2)*np.exp(-((x+x0)-mu2)**2/(2*sigma**2))+c))

c is redundant. x0 is redundant. sigma needs to vary between the two peaks. Peaks need outer amplitudes.

The following illustrates a fit based on the section that you chose, 40-70:

from io import StringIO

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

with StringIO(content) as f:
    bf1 = np.loadtxt(f)


def double_gaussian(
    x: np.ndarray,
    a: float, mu1: float, sigma1: float,
    b: float, mu2: float, sigma2: float,
    y0: float,
):
    return (
        y0
        + a/np.sqrt(2*np.pi*sigma1**2) * np.exp(-(x-mu1)**2/2/sigma1**2)
        + b/np.sqrt(2*np.pi*sigma2**2) * np.exp(-(x-mu2)**2/2/sigma2**2)
    )


def fit_double_gaussian(array, p0):
    x, y = array.T
    popt, _ = curve_fit(double_gaussian, *array[40:70, :].T, p0)
    print(popt)
    plt.scatter(x, y, label="data", s=2)
    plt.plot(x, double_gaussian(x, *popt), label="fit")
    plt.legend(loc=0)
    plt.show()


fit_double_gaussian(bf1, (5e-9, 685, 5,
                          1e-9, 775, 4, 0))

Bad, but headed in the right direction. Without slicing your input data, if you care about the whole array then a supergaussian fit is better:

from io import StringIO

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

with StringIO(content) as f:
    bf1 = np.loadtxt(f)


def double_gaussian(
    x: np.ndarray,
    a: float, mu1: float, sigma1: float, p1: float,
    b: float, mu2: float, sigma2: float, p2: float,
    y0: float,
):
    return (
        y0
        + a * np.exp(-np.abs(x-mu1)**p1/sigma1)
        + b * np.exp(-np.abs(x-mu2)**p2/sigma2)
    )


def fit_double_gaussian(array, p0):
    x, y = array.T
    popt, _ = curve_fit(
        double_gaussian, *array.T,
        p0=p0,
        # bounds=(
        #     (0,    600,    1,   0,
        #      0,    600,    1,   0,   0),
        #     (1e-9, 800,  100,   3,
        #      1e-9, 800,  100,   3,   3e-12),
        # )
    )
    print(popt)
    plt.scatter(x, y, label="data", s=2)
    plt.plot(x, double_gaussian(x, *popt), label="fit")
    plt.legend(loc=0)
    plt.show()


fit_double_gaussian(bf1, (2e-11, 685, 70, 1,
                          2e-11, 775, 40, 1, 3e-12))