0
# Confidence Intervals
1
2
LMFIT provides comprehensive tools for confidence interval estimation and error analysis, going beyond simple covariance-based uncertainties to explore parameter space and provide robust confidence bounds even for challenging nonlinear problems.
3
4
## Capabilities
5
6
### Confidence Interval Functions
7
8
Statistical analysis tools for determining parameter confidence bounds through parameter space exploration.
9
10
```python { .api }
11
def conf_interval(minimizer, result, p_names=None, sigmas=(0.674, 0.95, 0.997),
12
trace=False, maxiter=200, prob_func=None, verbose=False):
13
"""
14
Calculate confidence intervals for parameters using profile likelihood.
15
16
Args:
17
minimizer (Minimizer): Minimizer instance used for original fit
18
result (MinimizerResult): Results from minimization
19
p_names (list): Parameter names to analyze (all varied parameters if None)
20
sigmas (tuple): Confidence levels as probabilities (default: 1σ, 2σ, 3σ)
21
trace (bool): Return full trace of parameter exploration
22
maxiter (int): Maximum iterations for each confidence bound
23
prob_func (callable): Function to calculate probability from chi-squared
24
verbose (bool): Print progress information
25
26
Returns:
27
dict: Parameter names mapped to confidence interval dictionaries
28
Each parameter dict contains confidence bounds at each sigma level
29
"""
30
31
def conf_interval2d(minimizer, result, x_name, y_name, nx=10, ny=10, limits=None):
32
"""
33
Calculate 2D confidence regions for parameter pairs.
34
35
Args:
36
minimizer (Minimizer): Minimizer instance from original fit
37
result (MinimizerResult): Results from minimization
38
x_name (str): Name of first parameter
39
y_name (str): Name of second parameter
40
nx (int): Number of grid points in x direction
41
ny (int): Number of grid points in y direction
42
limits (tuple): ((x_min, x_max), (y_min, y_max)) for grid bounds
43
44
Returns:
45
tuple: (x_values, y_values, probability_grid)
46
x_values: Array of x parameter values
47
y_values: Array of y parameter values
48
probability_grid: 2D array of probability values
49
"""
50
51
def f_compare(best_fit, new_fit):
52
"""
53
F-test comparison of two fits.
54
55
Args:
56
best_fit (MinimizerResult): Results from better fit
57
new_fit (MinimizerResult): Results from comparison fit
58
59
Returns:
60
tuple: (f_statistic, p_value)
61
"""
62
```
63
64
### Utility Functions
65
66
Helper functions for confidence interval analysis and parameter management.
67
68
```python { .api }
69
def copy_vals(params):
70
"""
71
Save current parameter values and uncertainties.
72
73
Args:
74
params (Parameters): Parameters to save
75
76
Returns:
77
dict: Saved parameter state
78
"""
79
80
def restore_vals(tmp_params, params):
81
"""
82
Restore parameter values from saved state.
83
84
Args:
85
tmp_params (dict): Saved parameter state
86
params (Parameters): Parameters to restore
87
"""
88
89
def map_trace_to_names(trace, params):
90
"""
91
Map trace array to parameter names.
92
93
Args:
94
trace (array): Parameter trace array
95
params (Parameters): Parameter object
96
97
Returns:
98
dict: Parameter names mapped to trace values
99
"""
100
```
101
102
## Usage Examples
103
104
### Basic Confidence Intervals
105
106
```python
107
import numpy as np
108
from lmfit import minimize, Parameters, conf_interval, fit_report
109
110
def residual(params, x, data):
111
"""Exponential decay model"""
112
a = params['amplitude']
113
t = params['decay']
114
c = params['offset']
115
model = a * np.exp(-t * x) + c
116
return model - data
117
118
# Generate sample data
119
x = np.linspace(0, 10, 101)
120
data = 5.0 * np.exp(-0.8 * x) + 1.0 + np.random.normal(size=101, scale=0.1)
121
122
# Set up parameters and fit
123
params = Parameters()
124
params.add('amplitude', value=10, min=0)
125
params.add('decay', value=1, min=0)
126
params.add('offset', value=1)
127
128
# Perform fit
129
result = minimize(residual, params, args=(x, data))
130
131
# Calculate confidence intervals
132
ci = conf_interval(result.minimizer, result)
133
134
# Print results with confidence intervals
135
print(fit_report(result, show_correl=False))
136
print("\nConfidence Intervals:")
137
for param_name, bounds in ci.items():
138
print(f"{param_name}:")
139
for sigma, (lower, upper) in bounds.items():
140
percent = sigma * 100
141
print(f" {percent:.1f}%: [{lower:.4f}, {upper:.4f}]")
142
```
143
144
### Custom Confidence Levels
145
146
```python
147
# Calculate confidence intervals at custom probability levels
148
custom_sigmas = (0.68, 0.90, 0.95, 0.99) # 68%, 90%, 95%, 99%
149
ci = conf_interval(result.minimizer, result, sigmas=custom_sigmas, verbose=True)
150
151
# Process results
152
for param_name, bounds in ci.items():
153
print(f"\n{param_name} confidence intervals:")
154
for prob, (lower, upper) in bounds.items():
155
print(f" {prob*100:.0f}%: [{lower:.6f}, {upper:.6f}]")
156
```
157
158
### Confidence Intervals for Specific Parameters
159
160
```python
161
# Calculate confidence intervals only for specific parameters
162
specific_params = ['amplitude', 'decay']
163
ci_specific = conf_interval(result.minimizer, result,
164
p_names=specific_params,
165
maxiter=500)
166
167
# Get trace information for detailed analysis
168
ci_with_trace = conf_interval(result.minimizer, result,
169
p_names=['amplitude'],
170
trace=True, verbose=True)
171
172
# Access trace data
173
amplitude_trace = ci_with_trace['amplitude']['trace']
174
print(f"Trace points explored: {len(amplitude_trace)}")
175
```
176
177
### 2D Confidence Regions
178
179
```python
180
# Calculate 2D confidence region for parameter pair
181
x_vals, y_vals, prob_grid = conf_interval2d(result.minimizer, result,
182
'amplitude', 'decay',
183
nx=20, ny=20)
184
185
# Plot 2D confidence region
186
import matplotlib.pyplot as plt
187
188
plt.figure(figsize=(8, 6))
189
contour = plt.contour(x_vals, y_vals, prob_grid,
190
levels=[0.68, 0.95, 0.997],
191
colors=['blue', 'green', 'red'])
192
plt.clabel(contour, inline=True, fontsize=10,
193
fmt={0.68: '68%', 0.95: '95%', 0.997: '99.7%'})
194
195
# Mark best-fit point
196
best_amp = result.params['amplitude'].value
197
best_decay = result.params['decay'].value
198
plt.plot(best_amp, best_decay, 'ko', markersize=8, label='Best fit')
199
200
plt.xlabel('Amplitude')
201
plt.ylabel('Decay Rate')
202
plt.title('2D Confidence Regions')
203
plt.legend()
204
plt.grid(True, alpha=0.3)
205
plt.show()
206
```
207
208
### Advanced Confidence Analysis
209
210
```python
211
from lmfit import Minimizer
212
213
# Custom probability function for confidence intervals
214
def custom_prob_func(nvarys, ndata, chimin, chisqr):
215
"""Custom probability calculation"""
216
dof = ndata - nvarys
217
delta_chi = chisqr - chimin
218
# Custom logic for probability calculation
219
return np.exp(-0.5 * delta_chi)
220
221
# Create minimizer instance for detailed control
222
minimizer = Minimizer(residual, params, fcn_args=(x, data))
223
result = minimizer.minimize()
224
225
# Calculate confidence intervals with custom probability function
226
ci_custom = conf_interval(minimizer, result,
227
prob_func=custom_prob_func,
228
maxiter=1000, verbose=True)
229
```
230
231
### Confidence Intervals with Model Interface
232
233
```python
234
from lmfit import Model
235
from lmfit.models import ExponentialModel
236
237
# Using Model interface
238
exp_model = ExponentialModel()
239
params = exp_model.guess(data, x=x)
240
241
# Fit model
242
result = exp_model.fit(data, params, x=x)
243
244
# Calculate confidence intervals for model parameters
245
ci = conf_interval(result.minimizer, result)
246
247
# Report with confidence intervals included
248
from lmfit import ci_report
249
print(result.fit_report())
250
print("\n" + ci_report(ci))
251
```
252
253
### Comparing Fits with F-test
254
255
```python
256
# Compare nested models using F-test
257
from lmfit.models import LinearModel, QuadraticModel
258
259
# Fit simple linear model
260
linear_model = LinearModel()
261
linear_params = linear_model.guess(data, x=x)
262
linear_result = linear_model.fit(data, linear_params, x=x)
263
264
# Fit more complex quadratic model
265
quad_model = QuadraticModel()
266
quad_params = quad_model.guess(data, x=x)
267
quad_result = quad_model.fit(data, quad_params, x=x)
268
269
# F-test to compare models
270
f_stat, p_value = f_compare(linear_result, quad_result)
271
print(f"F-statistic: {f_stat:.4f}")
272
print(f"P-value: {p_value:.6f}")
273
274
if p_value < 0.05:
275
print("Quadratic model is significantly better")
276
else:
277
print("Linear model is adequate")
278
```
279
280
### Error Propagation with Uncertainties
281
282
```python
283
# Create uncertainties variables for error propagation
284
uvars = result.params.create_uvars()
285
286
# Use uncertainties for calculations with error propagation
287
amplitude_u = uvars['amplitude']
288
decay_u = uvars['decay']
289
offset_u = uvars['offset']
290
291
# Calculate derived quantities with propagated uncertainties
292
half_life = np.log(2) / decay_u
293
initial_value = amplitude_u + offset_u
294
295
print(f"Half-life: {half_life}")
296
print(f"Initial value: {initial_value}")
297
```
298
299
### Saving and Loading Confidence Results
300
301
```python
302
import json
303
304
# Convert confidence intervals to JSON-serializable format
305
ci_data = {}
306
for param_name, bounds in ci.items():
307
ci_data[param_name] = {}
308
for sigma, (lower, upper) in bounds.items():
309
ci_data[param_name][str(sigma)] = [lower, upper]
310
311
# Save to file
312
with open('confidence_intervals.json', 'w') as f:
313
json.dump(ci_data, f, indent=2)
314
315
# Load from file
316
with open('confidence_intervals.json', 'r') as f:
317
loaded_ci = json.load(f)
318
```