0
# Modeling and Fitting
1
2
Model fitting framework with astronomical models (Gaussian, Sérsic, etc.) and fitting algorithms for data analysis and parameter estimation.
3
4
## Capabilities
5
6
### Base Model Classes
7
8
Foundation classes for creating and manipulating models with parameters, constraints, and evaluation methods.
9
10
```python { .api }
11
class Model:
12
"""
13
Base class for all models.
14
15
Parameters:
16
- *args: positional arguments for model parameters
17
- **kwargs: keyword arguments for model parameters and metadata
18
"""
19
def __init__(self, *args, **kwargs): ...
20
21
def __call__(self, *args, **kwargs):
22
"""
23
Evaluate model at given inputs.
24
25
Parameters:
26
- *args: input coordinate arrays
27
- **kwargs: additional evaluation options
28
29
Returns:
30
ndarray: model values
31
"""
32
33
def evaluate(self, *args, **kwargs):
34
"""Evaluate model (implementation method)."""
35
36
def copy(self):
37
"""Create a copy of the model."""
38
39
def rename(self, name):
40
"""Rename the model."""
41
42
@property
43
def param_names(self):
44
"""List of parameter names."""
45
46
@property
47
def parameters(self):
48
"""Array of parameter values."""
49
50
@parameters.setter
51
def parameters(self, values):
52
"""Set parameter values."""
53
54
def __add__(self, other):
55
"""Add models together."""
56
57
def __sub__(self, other):
58
"""Subtract models."""
59
60
def __mul__(self, other):
61
"""Multiply models."""
62
63
def __truediv__(self, other):
64
"""Divide models."""
65
66
def __pow__(self, other):
67
"""Raise model to power."""
68
69
class Fittable1DModel(Model):
70
"""Base class for fittable 1D models."""
71
def __init__(self, *args, **kwargs): ...
72
73
class Fittable2DModel(Model):
74
"""Base class for fittable 2D models."""
75
def __init__(self, *args, **kwargs): ...
76
77
class CompoundModel(Model):
78
"""Model composed of multiple sub-models."""
79
def __init__(self, op, left, right, name=None): ...
80
```
81
82
### Model Parameters
83
84
Parameter class for managing model parameters with bounds, constraints, and units.
85
86
```python { .api }
87
class Parameter:
88
"""
89
Model parameter with value, bounds, and constraints.
90
91
Parameters:
92
- name: parameter name
93
- default: default value
94
- description: parameter description
95
- unit: parameter unit
96
- min: minimum allowed value
97
- max: maximum allowed value
98
- bounds: (min, max) bounds tuple
99
- tied: tie parameter to another parameter or function
100
- fixed: whether parameter is fixed during fitting
101
"""
102
def __init__(self, name='', default=0, description='', unit=None,
103
min=None, max=None, bounds=None, tied=False, fixed=False): ...
104
105
@property
106
def value(self):
107
"""Parameter value."""
108
109
@value.setter
110
def value(self, val):
111
"""Set parameter value."""
112
113
@property
114
def bounds(self):
115
"""Parameter bounds (min, max)."""
116
117
@property
118
def tied(self):
119
"""Parameter tied constraint."""
120
121
@property
122
def fixed(self):
123
"""Whether parameter is fixed."""
124
125
def copy(self):
126
"""Create copy of parameter."""
127
```
128
129
### 1D Astronomical Models
130
131
One-dimensional models commonly used in astronomical data analysis.
132
133
```python { .api }
134
class Gaussian1D(Fittable1DModel):
135
"""
136
1D Gaussian model.
137
138
Parameters:
139
- amplitude: Gaussian amplitude
140
- mean: Gaussian center
141
- stddev: Gaussian standard deviation
142
"""
143
def __init__(self, amplitude=1, mean=0, stddev=1, **kwargs): ...
144
145
class Lorentz1D(Fittable1DModel):
146
"""
147
1D Lorentzian model.
148
149
Parameters:
150
- amplitude: Lorentzian amplitude
151
- x_0: Lorentzian center
152
- fwhm: full width at half maximum
153
"""
154
def __init__(self, amplitude=1, x_0=0, fwhm=1, **kwargs): ...
155
156
class Voigt1D(Fittable1DModel):
157
"""
158
1D Voigt profile (convolution of Gaussian and Lorentzian).
159
160
Parameters:
161
- x_0: center position
162
- amplitude_L: Lorentzian amplitude
163
- fwhm_L: Lorentzian FWHM
164
- fwhm_G: Gaussian FWHM
165
"""
166
def __init__(self, x_0=0, amplitude_L=1, fwhm_L=1, fwhm_G=1, **kwargs): ...
167
168
class Moffat1D(Fittable1DModel):
169
"""
170
1D Moffat model for stellar profiles.
171
172
Parameters:
173
- amplitude: profile amplitude
174
- x_0: profile center
175
- gamma: profile width parameter
176
- alpha: profile power index
177
"""
178
def __init__(self, amplitude=1, x_0=0, gamma=1, alpha=1, **kwargs): ...
179
180
class Sersic1D(Fittable1DModel):
181
"""
182
1D Sérsic profile for galaxy surface brightness.
183
184
Parameters:
185
- amplitude: profile amplitude
186
- r_eff: effective radius
187
- n: Sérsic index
188
"""
189
def __init__(self, amplitude=1, r_eff=1, n=4, **kwargs): ...
190
191
class PowerLaw1D(Fittable1DModel):
192
"""
193
1D power law model.
194
195
Parameters:
196
- amplitude: amplitude at reference point
197
- x_0: reference point
198
- alpha: power law index
199
"""
200
def __init__(self, amplitude=1, x_0=1, alpha=1, **kwargs): ...
201
202
class BlackBody(Fittable1DModel):
203
"""
204
Blackbody spectrum model.
205
206
Parameters:
207
- temperature: blackbody temperature
208
- bolometric_flux: total bolometric flux
209
"""
210
def __init__(self, temperature=5000*u.K, bolometric_flux=1*u.erg/u.s/u.cm**2, **kwargs): ...
211
212
class Polynomial1D(Fittable1DModel):
213
"""
214
1D polynomial model.
215
216
Parameters:
217
- degree: polynomial degree
218
- c0, c1, c2, ...: polynomial coefficients
219
"""
220
def __init__(self, degree, **kwargs): ...
221
```
222
223
### 2D Astronomical Models
224
225
Two-dimensional models for image analysis and surface brightness modeling.
226
227
```python { .api }
228
class Gaussian2D(Fittable2DModel):
229
"""
230
2D Gaussian model.
231
232
Parameters:
233
- amplitude: Gaussian amplitude
234
- x_mean: x-coordinate of center
235
- y_mean: y-coordinate of center
236
- x_stddev: standard deviation in x
237
- y_stddev: standard deviation in y
238
- theta: rotation angle
239
"""
240
def __init__(self, amplitude=1, x_mean=0, y_mean=0, x_stddev=1, y_stddev=1, theta=0, **kwargs): ...
241
242
class Moffat2D(Fittable2DModel):
243
"""
244
2D Moffat model for stellar PSF.
245
246
Parameters:
247
- amplitude: profile amplitude
248
- x_0: x-coordinate of center
249
- y_0: y-coordinate of center
250
- gamma: profile width parameter
251
- alpha: power index
252
"""
253
def __init__(self, amplitude=1, x_0=0, y_0=0, gamma=1, alpha=1, **kwargs): ...
254
255
class Sersic2D(Fittable2DModel):
256
"""
257
2D Sérsic profile for galaxy modeling.
258
259
Parameters:
260
- amplitude: surface brightness at effective radius
261
- r_eff: effective radius
262
- n: Sérsic index
263
- x_0: x-coordinate of center
264
- y_0: y-coordinate of center
265
- ellip: ellipticity
266
- theta: position angle
267
"""
268
def __init__(self, amplitude=1, r_eff=1, n=4, x_0=0, y_0=0, ellip=0, theta=0, **kwargs): ...
269
270
class AiryDisk2D(Fittable2DModel):
271
"""
272
2D Airy disk for diffraction-limited PSF.
273
274
Parameters:
275
- amplitude: peak amplitude
276
- x_0: x-coordinate of center
277
- y_0: y-coordinate of center
278
- radius: Airy disk radius
279
"""
280
def __init__(self, amplitude=1, x_0=0, y_0=0, radius=1, **kwargs): ...
281
282
class Disk2D(Fittable2DModel):
283
"""
284
2D disk (top-hat) model.
285
286
Parameters:
287
- amplitude: disk amplitude
288
- x_0: x-coordinate of center
289
- y_0: y-coordinate of center
290
- R_0: disk radius
291
"""
292
def __init__(self, amplitude=1, x_0=0, y_0=0, R_0=1, **kwargs): ...
293
294
class Polynomial2D(Fittable2DModel):
295
"""
296
2D polynomial model.
297
298
Parameters:
299
- degree: polynomial degree
300
- c0_0, c1_0, c0_1, ...: polynomial coefficients
301
"""
302
def __init__(self, degree, **kwargs): ...
303
```
304
305
### Fitting Algorithms
306
307
Various fitting algorithms for parameter estimation and model optimization.
308
309
```python { .api }
310
class Fitter:
311
"""Base class for fitting algorithms."""
312
def __call__(self, model, x, y, weights=None, **kwargs):
313
"""
314
Fit model to data.
315
316
Parameters:
317
- model: model to fit
318
- x: independent variable data
319
- y: dependent variable data
320
- weights: data weights
321
- **kwargs: fitter-specific options
322
323
Returns:
324
Model: fitted model instance
325
"""
326
327
class LinearLSQFitter(Fitter):
328
"""Linear least squares fitter for linear models."""
329
def __init__(self, calc_uncertainties=False): ...
330
331
class LevMarLSQFitter(Fitter):
332
"""Levenberg-Marquardt least squares fitter."""
333
def __init__(self, calc_uncertainties=False): ...
334
335
class SLSQPLSQFitter(Fitter):
336
"""Sequential Least Squares Programming fitter."""
337
def __init__(self, calc_uncertainties=False): ...
338
339
class SimplexLSQFitter(Fitter):
340
"""Simplex algorithm fitter."""
341
def __init__(self, calc_uncertainties=False): ...
342
343
class DifferentialEvolutionLSQFitter(Fitter):
344
"""Differential evolution global optimizer."""
345
def __init__(self, calc_uncertainties=False): ...
346
347
class LMLSQFitter(Fitter):
348
"""Levenberg-Marquardt fitter using scipy.optimize.leastsq."""
349
def __init__(self, calc_uncertainties=False): ...
350
```
351
352
### Model Utilities
353
354
Utility functions for model evaluation, parameter handling, and model combinations.
355
356
```python { .api }
357
def custom_model(func, fit_deriv=None):
358
"""
359
Create custom model from function.
360
361
Parameters:
362
- func: function that defines the model
363
- fit_deriv: function that computes derivatives
364
365
Returns:
366
Model: custom model class
367
"""
368
369
def bind_bounding_box(model, bounding_box, order='C'):
370
"""Bind a bounding box to a model."""
371
372
def fix_inputs(model, fixed_inputs):
373
"""Fix certain inputs of a model."""
374
```
375
376
## Usage Examples
377
378
### Basic Model Fitting
379
380
```python
381
from astropy.modeling import models, fitting
382
import numpy as np
383
import matplotlib.pyplot as plt
384
385
# Generate synthetic data
386
np.random.seed(42)
387
x = np.linspace(-3, 3, 100)
388
y_true = 2 * np.exp(-0.5 * ((x - 0.5) / 0.8)**2) + 0.1
389
noise = np.random.normal(0, 0.1, len(x))
390
y = y_true + noise
391
392
# Create and fit Gaussian model
393
gaussian_model = models.Gaussian1D(amplitude=2, mean=0.5, stddev=0.8)
394
fitter = fitting.LevMarLSQFitter()
395
fitted_model = fitter(gaussian_model, x, y)
396
397
print(f"Fitted parameters:")
398
print(f" Amplitude: {fitted_model.amplitude.value:.3f}")
399
print(f" Mean: {fitted_model.mean.value:.3f}")
400
print(f" Stddev: {fitted_model.stddev.value:.3f}")
401
402
# Evaluate fitted model
403
y_fit = fitted_model(x)
404
```
405
406
### 2D Surface Brightness Modeling
407
408
```python
409
from astropy.modeling import models, fitting
410
import numpy as np
411
412
# Create 2D coordinate grids
413
y, x = np.mgrid[0:100, 0:100]
414
415
# Generate synthetic galaxy with Sérsic profile
416
true_model = models.Sersic2D(amplitude=1000, r_eff=15, n=2,
417
x_0=50, y_0=45, ellip=0.3, theta=0.5)
418
galaxy_data = true_model(x, y)
419
420
# Add noise
421
noise = np.random.poisson(galaxy_data + 100) - 100 # Poisson noise
422
galaxy_data_noisy = galaxy_data + noise
423
424
# Fit Sérsic model
425
initial_model = models.Sersic2D(amplitude=800, r_eff=10, n=1,
426
x_0=48, y_0=47, ellip=0.2, theta=0)
427
fitter = fitting.LevMarLSQFitter()
428
fitted_galaxy = fitter(initial_model, x, y, galaxy_data_noisy)
429
430
print("Galaxy fitting results:")
431
print(f" Effective radius: {fitted_galaxy.r_eff.value:.2f} pixels")
432
print(f" Sérsic index: {fitted_galaxy.n.value:.2f}")
433
print(f" Ellipticity: {fitted_galaxy.ellip.value:.3f}")
434
print(f" Position angle: {fitted_galaxy.theta.value:.3f} radians")
435
```
436
437
### Spectral Line Fitting
438
439
```python
440
# Fit multiple spectral lines
441
wavelength = np.linspace(6540, 6580, 200)
442
443
# Create composite model with multiple Gaussian lines
444
h_alpha = models.Gaussian1D(amplitude=100, mean=6562.8, stddev=0.5, name='H_alpha')
445
nii_6548 = models.Gaussian1D(amplitude=30, mean=6548.1, stddev=0.4, name='NII_6548')
446
nii_6583 = models.Gaussian1D(amplitude=90, mean=6583.5, stddev=0.4, name='NII_6583')
447
continuum = models.Polynomial1D(degree=1, c0=10, c1=0, name='continuum')
448
449
# Combine models
450
line_model = continuum + h_alpha + nii_6548 + nii_6583
451
452
# Generate synthetic spectrum
453
spectrum = line_model(wavelength)
454
spectrum += np.random.normal(0, 2, len(wavelength)) # Add noise
455
456
# Fit the composite model
457
fitter = fitting.LevMarLSQFitter()
458
fitted_lines = fitter(line_model, wavelength, spectrum)
459
460
# Extract individual line fluxes
461
h_alpha_flux = fitted_lines['H_alpha'].amplitude * fitted_lines['H_alpha'].stddev * np.sqrt(2*np.pi)
462
nii_ratio = fitted_lines['NII_6583'].amplitude / fitted_lines['H_alpha'].amplitude
463
464
print(f"H-alpha flux: {h_alpha_flux:.1f}")
465
print(f"[NII]/H-alpha ratio: {nii_ratio:.3f}")
466
```
467
468
### Parameter Constraints and Bounds
469
470
```python
471
# Model with parameter constraints
472
psf_model = models.Moffat2D(amplitude=1000, x_0=50, y_0=50, gamma=2, alpha=2)
473
474
# Set parameter bounds
475
psf_model.amplitude.bounds = (100, 10000)
476
psf_model.x_0.bounds = (48, 52)
477
psf_model.y_0.bounds = (48, 52)
478
psf_model.gamma.bounds = (0.5, 5)
479
psf_model.alpha.bounds = (1, 10)
480
481
# Fix certain parameters
482
psf_model.x_0.fixed = True
483
psf_model.y_0.fixed = True
484
485
# Tie parameters together (force circular PSF)
486
def tie_gamma_to_alpha(model):
487
return model.alpha * 0.5
488
489
psf_model.gamma.tied = tie_gamma_to_alpha
490
491
# Fit with constraints
492
fitter = fitting.LevMarLSQFitter()
493
# fitted_psf = fitter(psf_model, x, y, psf_data)
494
```
495
496
### Custom Model Creation
497
498
```python
499
from astropy.modeling import Fittable1DModel, Parameter
500
import numpy as np
501
502
class ExponentialCutoff(Fittable1DModel):
503
\"\"\"Custom exponentially cutoff power law model.\"\"\"
504
505
amplitude = Parameter(default=1, bounds=(0, None))
506
x_0 = Parameter(default=1, bounds=(0, None))
507
alpha = Parameter(default=1)
508
cutoff = Parameter(default=1, bounds=(0, None))
509
510
@staticmethod
511
def evaluate(x, amplitude, x_0, alpha, cutoff):
512
\"\"\"Model evaluation function.\"\"\"
513
return amplitude * (x / x_0)**(-alpha) * np.exp(-x / cutoff)
514
515
# Use custom model
516
x_data = np.logspace(0, 3, 50)
517
custom_model = ExponentialCutoff(amplitude=100, x_0=10, alpha=2, cutoff=200)
518
y_model = custom_model(x_data)
519
520
print(f"Custom model at x=50: {custom_model(50):.2f}")
521
```
522
523
### Model Arithmetic and Combinations
524
525
```python
526
# Create individual models
527
gaussian = models.Gaussian1D(amplitude=10, mean=0, stddev=1)
528
linear = models.Linear1D(slope=0.5, intercept=2)
529
constant = models.Const1D(amplitude=5)
530
531
# Combine models using arithmetic
532
combined_additive = gaussian + linear + constant
533
combined_multiplicative = gaussian * linear
534
combined_complex = (gaussian + constant) * linear
535
536
# Evaluate combined models
537
x_test = np.linspace(-5, 5, 100)
538
y_combined = combined_additive(x_test)
539
540
print(f"Number of parameters in combined model: {len(combined_additive.parameters)}")
541
print(f"Parameter names: {combined_additive.param_names}")
542
```
543
544
### Advanced Fitting with Uncertainties
545
546
```python
547
# Fitting with parameter uncertainties
548
fitter_with_errs = fitting.LevMarLSQFitter(calc_uncertainties=True)
549
550
# Generate data with known uncertainties
551
x_data = np.linspace(0, 10, 50)
552
true_model = models.Polynomial1D(degree=2, c0=1, c1=2, c2=-0.1)
553
y_true = true_model(x_data)
554
y_errors = 0.1 + 0.05 * np.abs(y_true) # Heteroscedastic errors
555
y_data = y_true + np.random.normal(0, y_errors)
556
557
# Fit with error weights
558
weights = 1.0 / y_errors**2
559
poly_model = models.Polynomial1D(degree=2)
560
fitted_poly = fitter_with_errs(poly_model, x_data, y_data, weights=weights)
561
562
# Extract parameter uncertainties (if supported by fitter)
563
if hasattr(fitter_with_errs, 'fit_info') and fitter_with_errs.fit_info.get('param_cov') is not None:
564
param_errors = np.sqrt(np.diag(fitter_with_errs.fit_info['param_cov']))
565
print("Parameter uncertainties:")
566
for name, value, error in zip(fitted_poly.param_names, fitted_poly.parameters, param_errors):
567
print(f" {name}: {value:.4f} ± {error:.4f}")
568
```