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>
|
