File: doc_model_uncertainty2.py

package info (click to toggle)
lmfit-py 1.3.3-4
  • links: PTS, VCS
  • area: main
  • in suites: sid, trixie
  • size: 2,332 kB
  • sloc: python: 13,071; makefile: 130; sh: 30
file content (79 lines) | stat: -rw-r--r-- 3,148 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
# <examples/doc_model_uncertainty2.py>
import matplotlib.pyplot as plt
import numpy as np

from lmfit.models import ExponentialModel, GaussianModel

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

model = (GaussianModel(prefix='g1_') +
         GaussianModel(prefix='g2_') +
         ExponentialModel(prefix='bkg_'))

params = model.make_params(bkg_amplitude=100, bkg_decay=80,
                           g1_amplitude=3000,
                           g1_center=100,
                           g1_sigma=10,
                           g2_amplitude=3000,
                           g2_center=150,
                           g2_sigma=10)

result = model.fit(y, params, x=x)
print(result.fit_report(min_correl=0.5))

comps = result.eval_components(x=x)
dely = result.eval_uncertainty(sigma=3)

fig, axes = plt.subplots(2, 2, figsize=(12.8, 9.6))

axes[0][0].plot(x, y, 'o', color='#99002299', markersize=3, label='data')
axes[0][0].plot(x, result.best_fit, '-', label='best fit')
axes[0][0].plot(x, result.init_fit, '--', label='initial fit')
axes[0][0].set_title('data, initial fit, and best-fit')
axes[0][0].legend()

axes[0][1].plot(x, y, 'o', color='#99002299', markersize=3, label='data')
axes[0][1].plot(x, result.best_fit, '-', label='best fit')
axes[0][1].fill_between(x, result.best_fit-dely, result.best_fit+dely,
                        color="#8A8A8A", label=r'3-$\sigma$ band')
axes[0][1].set_title('data, best-fit, and uncertainty band')
axes[0][1].legend()

axes[1][0].plot(x, result.best_fit, '-', label=r'best fit, 3-$\sigma$ band')
axes[1][0].fill_between(x,
                        result.best_fit-result.dely,
                        result.best_fit+result.dely,
                        color="#8A8A8A")

axes[1][0].plot(x, comps['bkg_'], label=r'background, 3-$\sigma$ band')
axes[1][0].fill_between(x,
                        comps['bkg_']-result.dely_comps['bkg_'],
                        comps['bkg_']+result.dely_comps['bkg_'],
                        color="#8A8A8A")

axes[1][0].plot(x, comps['g1_'], label=r'Gaussian #1, 3-$\sigma$ band')
axes[1][0].fill_between(x,
                        comps['g1_']-result.dely_comps['g1_'],
                        comps['g1_']+result.dely_comps['g1_'],
                        color="#8A8A8A")

axes[1][0].plot(x, comps['g2_'], label=r'Gaussian #2, 3-$\sigma$ band')
axes[1][0].fill_between(x,
                        comps['g2_']-result.dely_comps['g2_'],
                        comps['g2_']+result.dely_comps['g2_'],
                        color="#8A8A8A")
axes[1][0].set_title('model components with uncertainty bands')
axes[1][0].legend()

axes[1][1].plot(x, result.best_fit, '-', label='best fit')
axes[1][1].plot(x, 10*result.dely, label=r'3-$\sigma$ total (x10)')
axes[1][1].plot(x, 10*result.dely_comps['bkg_'], label=r'3-$\sigma$ background (x10)')
axes[1][1].plot(x, 10*result.dely_comps['g1_'], label=r'3-$\sigma$ Gaussian #1 (x10)')
axes[1][1].plot(x, 10*result.dely_comps['g2_'], label=r'3-$\sigma$ Gaussian #2 (x10)')
axes[1][1].set_title('uncertainties for model components')
axes[1][1].legend()

plt.show()
# <end examples/doc_model_uncertainty2.py>