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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
"""
Fit using the Model interface
=============================
This notebook shows a simple example of using the ``lmfit.Model`` class. For
more information please refer to:
https://lmfit.github.io/lmfit-py/model.html#the-model-class.
"""
import numpy as np
from pandas import Series
from lmfit import Model, Parameter, report_fit
###############################################################################
# The ``Model`` class is a flexible, concise curve fitter. I will illustrate
# fitting example data to an exponential decay.
def decay(t, N, tau):
return N*np.exp(-t/tau)
###############################################################################
# The parameters are in no particular order. We'll need some example data. I
# will use ``N=7`` and ``tau=3``, and add a little noise.
t = np.linspace(0, 5, num=1000)
data = decay(t, 7, 3) + np.random.randn(t.size)
###############################################################################
# **Simplest Usage**
model = Model(decay, independent_vars=['t'])
result = model.fit(data, t=t, N=10, tau=1)
###############################################################################
# The Model infers the parameter names by inspecting the arguments of the
# function, ``decay``. Then I passed the independent variable, ``t``, and
# initial guesses for each parameter. A residual function is automatically
# defined, and a least-squared regression is performed.
#
# We can immediately see the best-fit values:
print(result.values)
###############################################################################
# and use these best-fit parameters for plotting with the ``plot`` function:
result.plot()
###############################################################################
# We can review the best-fit `Parameters` by accessing `result.params`:
result.params.pretty_print()
###############################################################################
# More information about the fit is stored in the result, which is an
# ``lmfit.MimimizerResult`` object (see:
# https://lmfit.github.io/lmfit-py/fitting.html#lmfit.minimizer.MinimizerResult)
###############################################################################
# **Specifying Bounds and Holding Parameters Constant**
#
# Above, the ``Model`` class implicitly builds ``Parameter`` objects from
# keyword arguments of ``fit`` that match the arguments of ``decay``. You can
# build the ``Parameter`` objects explicitly; the following is equivalent.
result = model.fit(data, t=t,
N=Parameter('N', value=10),
tau=Parameter('tau', value=1))
report_fit(result.params)
###############################################################################
# By building ``Parameter`` objects explicitly, you can specify bounds
# (``min``, ``max``) and set parameters constant (``vary=False``).
result = model.fit(data, t=t,
N=Parameter('N', value=7, vary=False),
tau=Parameter('tau', value=1, min=0))
report_fit(result.params)
###############################################################################
# **Defining Parameters in Advance**
#
# Passing parameters to ``fit`` can become unwieldy. As an alternative, you
# can extract the parameters from ``model`` like so, set them individually,
# and pass them to ``fit``.
params = model.make_params()
params['N'].value = 10
params['tau'].value = 1
params['tau'].min = 0
result = model.fit(data, params, t=t)
report_fit(result.params)
##############################################################################
# Keyword arguments override ``params``, resetting ``value`` and all other
# properties (``min``, ``max``, ``vary``).
result = model.fit(data, params, t=t, tau=1)
report_fit(result.params)
###############################################################################
# The input parameters are not modified by ``fit``. They can be reused,
# retaining the same initial value. If you want to use the result of one fit
# as the initial guess for the next, simply pass ``params=result.params``.
###############################################################################
# #TODO/FIXME: not sure if there ever way a "helpful exception", but currently
# #it raises a ``ValueError: The input contains nan values``.
#
# #*A Helpful Exception*
#
# #All this implicit magic makes it very easy for the user to neglect to set a
# #parameter. The ``fit`` function checks for this and raises a helpful exception.
# #result = model.fit(data, t=t, tau=1) # N unspecified
###############################################################################
# #An *extra* parameter that cannot be matched to the model function will
# #throw a ``UserWarning``, but it will not raise, leaving open the possibility
# #of unforeseen extensions calling for some parameters.
###############################################################################
# *Weighted Fits*
#
# Use the ``sigma`` argument to perform a weighted fit. If you prefer to think
# of the fit in term of ``weights``, ``sigma=1/weights``.
weights = np.arange(len(data))
result = model.fit(data, params, t=t, weights=weights)
report_fit(result.params)
###############################################################################
# *Handling Missing Data*
#
# By default, attempting to fit data that includes a ``NaN``, which
# conventionally indicates a "missing" observation, raises a lengthy exception.
# You can choose to ``omit`` (i.e., skip over) missing values instead.
data_with_holes = data.copy()
data_with_holes[[5, 500, 700]] = np.nan # Replace arbitrary values with NaN.
model = Model(decay, independent_vars=['t'], nan_policy='omit')
result = model.fit(data_with_holes, params, t=t)
report_fit(result.params)
###############################################################################
# If you don't want to ignore missing values, you can set the model to raise
# proactively, checking for missing values before attempting the fit.
#
# Uncomment to see the error
# #model = Model(decay, independent_vars=['t'], nan_policy='raise')
# #result = model.fit(data_with_holes, params, t=t)
#
# The default setting is ``nan_policy='raise'``, which does check for NaNs and
# raises an exception when present.
#
# Null-chekcing relies on ``pandas.isnull`` if it is available. If pandas
# cannot be imported, it silently falls back on ``numpy.isnan``.
###############################################################################
# *Data Alignment*
#
# Imagine a collection of time series data with different lengths. It would be
# convenient to define one sufficiently long array ``t`` and use it for each
# time series, regardless of length. ``pandas``
# (https://pandas.pydata.org/pandas-docs/stable/) provides tools for aligning
# indexed data. And, unlike most wrappers to ``scipy.leastsq``, ``Model`` can
# handle pandas objects out of the box, using its data alignment features.
#
# Here I take just a slice of the ``data`` and fit it to the full ``t``. It is
# automatically aligned to the correct section of ``t`` using Series' index.
model = Model(decay, independent_vars=['t'])
truncated_data = Series(data)[200:800] # data points 200-800
t = Series(t) # all 1000 points
result = model.fit(truncated_data, params, t=t)
report_fit(result.params)
###############################################################################
# Data with missing entries and an unequal length still aligns properly.
model = Model(decay, independent_vars=['t'], nan_policy='omit')
truncated_data_with_holes = Series(data_with_holes)[200:800]
result = model.fit(truncated_data_with_holes, params, t=t)
report_fit(result.params)
|
