How to use the lmfit.models.GaussianModel function in lmfit

'g2_center': values2['cen'],
                       'g2_amplitude': values2['amp'],
                       'g2_sigma': values2['sig']}
        pars = model.make_params()
        pars['g1_sigma'].set(2)
        pars['g1_center'].set(1)
        pars['g1_amplitude'].set(3)
        pars['g2_sigma'].set(1)
        pars['g2_center'].set(2.4)
        pars['g2_amplitude'].set(1)

        result = model.fit(data, pars, x=self.x)
        assert_results_close(result.values, true_values, rtol=0.01, atol=0.01)

        # without suffix, the names collide and Model should raise
        model1 = models.GaussianModel()
        model2 = models.GaussianModel()
        f = lambda: model1 + model2
        self.assertRaises(NameError, f)

def test_composite_has_bestvalues(self):
        # test that a composite model has non-empty best_values
        model1 = models.GaussianModel(prefix='g1_')
        model2 = models.GaussianModel(prefix='g2_')

        mod = model1 + model2
        pars = mod.make_params()

        values1 = dict(amplitude=7.10, center=1.1, sigma=2.40)
        values2 = dict(amplitude=12.2, center=2.5, sigma=0.5)
        data = (gaussian(x=self.x, **values1) + gaussian(x=self.x, **values2)
                + 0.1*self.noise)

        pars['g1_sigma'].set(value=2)
        pars['g1_center'].set(value=1, max=1.5)
        pars['g1_amplitude'].set(value=3)
        pars['g2_sigma'].set(value=1)
        pars['g2_center'].set(value=2.6, min=2.0)
        pars['g2_amplitude'].set(value=1)

def test_composite_plotting(self):
        # test that a composite model has non-empty best_values
        pytest.importorskip("matplotlib")
        import matplotlib
        matplotlib.use('Agg')

        model1 = models.GaussianModel(prefix='g1_')
        model2 = models.GaussianModel(prefix='g2_')

        mod = model1 + model2
        pars = mod.make_params()

        values1 = dict(amplitude=7.10, center=1.1, sigma=2.40)
        values2 = dict(amplitude=12.2, center=2.5, sigma=0.5)
        data = (gaussian(x=self.x, **values1) + gaussian(x=self.x, **values2)
                + 0.1*self.noise)

        pars['g1_sigma'].set(2)
        pars['g1_center'].set(1, max=1.5)
        pars['g1_amplitude'].set(3)
        pars['g2_sigma'].set(1)
        pars['g2_center'].set(2.6, min=2.0)
        pars['g2_amplitude'].set(1)

def test_sum_composite_models(self):
        # test components of composite model created adding composite model
        model1 = models.GaussianModel(prefix='g1_')
        model2 = models.GaussianModel(prefix='g2_')
        model3 = models.GaussianModel(prefix='g3_')
        model4 = models.GaussianModel(prefix='g4_')

        model_total1 = (model1 + model2) + model3
        for mod in [model1, model2, model3]:
            self.assertTrue(mod in model_total1.components)

        model_total2 = model1 + (model2 + model3)
        for mod in [model1, model2, model3]:
            self.assertTrue(mod in model_total2.components)

        model_total3 = (model1 + model2) + (model3 + model4)
        for mod in [model1, model2, model3, model4]:
            self.assertTrue(mod in model_total3.components)

def create_model_params(x, y):
    """Create the model and parameters."""
    exp_mod = ExponentialModel(prefix='exp_')
    params = exp_mod.guess(y, x=x)

    gauss1 = GaussianModel(prefix='g1_')
    params.update(gauss1.make_params())

    gauss2 = GaussianModel(prefix='g2_')
    params.update(gauss2.make_params())

    params['g1_center'].set(value=105, min=75, max=125)
    params['g1_sigma'].set(value=15, min=3)
    params['g1_amplitude'].set(value=2000, min=10)

    params['g2_center'].set(value=155, min=125, max=175)
    params['g2_sigma'].set(value=15, min=3)
    params['g2_amplitude'].set(value=2000, min=10)

    model = gauss1 + gauss2 + exp_mod
    return model, params

import numpy as np

from lmfit.models import GaussianModel

data = np.loadtxt('model1d_gauss.dat')
x = data[:, 0]
y = data[:, 1]

y[44] = np.nan
y[65] = np.nan

# nan_policy = 'raise'
# nan_policy = 'propagate'
nan_policy = 'omit'

gmodel = GaussianModel()
result = gmodel.fit(y, x=x, amplitude=5, center=6, sigma=1,
                    nan_policy=nan_policy)

print(result.fit_report())

# make sure nans are removed for plotting:
x_ = x[np.where(np.isfinite(y))]
y_ = y[np.where(np.isfinite(y))]

plt.plot(x_, y_, 'bo')
plt.plot(x_, result.init_fit, 'k--', label='initial fit')
plt.plot(x_, result.best_fit, 'r-', label='best fit')
plt.legend(loc='best')
plt.show()
# 

dat = np.loadtxt('NIST_Gauss2.dat')
x = dat[:, 1]
y = dat[:, 0]

exp_mod = ExponentialModel(prefix='exp_')
pars = exp_mod.guess(y, x=x)

gauss1 = GaussianModel(prefix='g1_')
pars.update(gauss1.make_params())

pars['g1_center'].set(value=105, min=75, max=125)
pars['g1_sigma'].set(value=15, min=3)
pars['g1_amplitude'].set(value=2000, min=10)

gauss2 = GaussianModel(prefix='g2_')
pars.update(gauss2.make_params())

pars['g2_center'].set(value=155, min=125, max=175)
pars['g2_sigma'].set(value=15, min=3)
pars['g2_amplitude'].set(value=2000, min=10)

mod = gauss1 + gauss2 + exp_mod

init = mod.eval(pars, x=x)
out = mod.fit(y, pars, x=x)

print(out.fit_report(min_correl=0.5))

fig, axes = plt.subplots(1, 2, figsize=(12.8, 4.8))
axes[0].plot(x, y, 'b')
axes[0].plot(x, init, 'k--', label='initial fit')

pars['g0_height'].set(height, min = 0)

        for n in range(nPeaks)[1:]:

            if n == nPeaks - 1:
                cMax = bins.max()

            else:
                cMax = old_div((xMaxs[n] + xMaxs[n + 1]),2)

            cMin = old_div((xMaxs[n] + xMaxs[n - 1]),2)

            center = xMaxs[n]
            height = yMaxs[n]

            gModN = GaussianModel(prefix = 'g%s_' % n)
            parsN = gModN.guess(frequencies, x = bins)
            parsN['g%s_center' % n].set(center, min = cMin, max = cMax)
            parsN['g%s_height' % n].set(height, min = 0)

            gMod += gModN
            pars.update(parsN)

        out = gMod.fit(frequencies, pars, x = bins)
        fit = out.best_fit
        resonance = []
        stderr = []
        fwhm = []
        sigma = []

        for n in range(nPeaks):
            resonance.append(out.params['g%s_center' % n].value)

def histyFit(frequencies, bins, nPeaks = 1, xMaxs = [], yMaxs = []):

    if nPeaks == 1:
        gMod = GaussianModel()
        pars = gMod.guess(frequencies, x = bins)
        out = gMod.fit(frequencies, pars, x = bins)#Performs the fit, based on initial guesses
        resonance = out.params['center'].value
        stderr = out.params['center'].stderr
        fwhm = out.params['fwhm'].value
        sigma = out.params['sigma'].value
        fit = out.best_fit

        print('\t\tAverage peakpos: %s +/- %s nm' % (resonance, stderr))
        print('\t\tFWHM: %s nm\n' % fwhm)

    else:
        gMod = GaussianModel(prefix = 'g0_')
        pars = gMod.guess(frequencies, x = bins)

        center = xMaxs[0]