Development documentation · 0.2.3 30c2f8ba · Intro examples tested with 0.2.3 · Version details · Known limitations

Fitting#

Regression, generalized-least-squares, and MCMC fitting APIs.

gwexpy.fitting.fit_series(series, model, x_range=None, sigma=None, cov=None, cost_function=None, p0=None, limits=None, fixed=None, **kwargs)#

Fit a Series object using iminuit.

Supports real and complex valued Series (simultaneous Re/Im fit).

Parameters:

series (Series) – Data series to fit.
model (callable or str) – Model function or name (e.g., “gaussian”, “power_law”).
x_range (tuple, optional) – (x_min, x_max) to crop data before fitting.
sigma (array-like or scalar, optional) – Per-point error estimates. Ignored if cov or cost_function is provided.
cov (BifrequencyMap or 2D ndarray, optional) – Covariance matrix for Generalized Least Squares (GLS) fitting. If provided, overrides sigma and uses GLS χ² minimization. Ignored if cost_function is provided.
cost_function (callable, optional) – User-defined cost function for Minuit. If provided, takes priority over sigma, cov, and automatic cost function selection. Must be callable with signature cost_function(*params) -> float. Should have errordef attribute (default: 1.0 for least squares). Parameter names are extracted via iminuit.util.describe().
p0 (dict or list, optional) – Initial parameter values.
limits (dict, optional) – Parameter limits, e.g., {“A”: (0, 100)}.
fixed (list, optional) – List of parameter names to fix during fit.
**kwargs – Additional arguments passed to Minuit.

Returns:

Object containing fit results, parameters, and plotting methods.

Return type:

FitResult

class gwexpy.fitting.FitResult(minuit_obj, model, x, y, dy=None, cost_func=None, x_label=None, y_label=None, x_kind=None, x_data=None, y_data=None, dy_data=None, x_fit_range=None, cov_inv=None, cov=None, unit=None, x_unit=None)#

Bases: object

Store fit outputs, metadata, and plotting helpers for a completed fit.

bode_plot(ax=None, num_points=1000, **kwargs)#: Create a Bode plot (Magnitude and Phase) for the fit result.

See also

Tutorial: Transfer Function

property chi2#: Chi-square value (valid only for LeastSquares-like costs).

property errors#: Parameter errors (dict).

property mcmc_chain#

Get the full MCMC chain (not flattened, not discarded).

Returns:: Shape (n_steps, n_walkers, n_params). Returns None if MCMC not run.
Return type:: ndarray

property model#

Fitted model function.

Can be called as model(x) to use best-fit parameters, or as model(x, **params) to use specific parameters. Returns a Quantity with units if the original data had units.

property ndof#: Number of degrees of freedom.

property parameter_intervals#

Get parameter confidence intervals from MCMC samples.

Returns 16th, 50th, and 84th percentiles for each parameter, corresponding to median and ±1σ bounds.

Returns:: Dictionary mapping parameter names to (lower, median, upper) tuples.
Return type:: dict
Raises:: RuntimeError – If run_mcmc() has not been called.

property params#: Best fit parameters (dict of ParameterValue).

plot(ax=None, num_points=1000, **kwargs)#

Plot data and best-fit curve.

For complex data, delegates to bode_plot().

plot_corner(show_titles=True, quantiles=None, **kwargs)#

Plot corner plot of MCMC samples.

Parameters:

show_titles (bool, optional) – Whether to show parameter value titles. Default is True.
quantiles (list, optional) – Quantiles for title display. Default is [0.16, 0.5, 0.84].
**kwargs – Additional arguments passed to corner.corner().

Returns:

figure – The corner plot figure.

Return type:

matplotlib.figure.Figure

plot_fit_band(ax=None, num_points=200, n_samples=100, alpha=0.3, **kwargs)#

Plot the fit with uncertainty band from MCMC samples.

Parameters:

ax (matplotlib.axes.Axes, optional) – Axes to plot on. If None, creates new figure.
num_points (int, optional) – Number of points for model curve. Default is 200.
n_samples (int, optional) – Number of MCMC samples to use for band. Default is 100.
alpha (float, optional) – Alpha transparency for uncertainty band. Default is 0.3.
**kwargs – Additional arguments passed to ax.fill_between().

Returns:

ax – The axes with the plot.

Return type:

matplotlib.axes.Axes

property reduced_chi2#: Reduced Chi-square value.

run_mcmc(n_walkers=32, n_steps=3000, burn_in=500, progress=True)#

Run MCMC using emcee starting from the best-fit parameters.

This method supports both standard least squares and GLS (Generalized Least Squares) error structures. If cov_inv is available, the log probability is computed using the full covariance structure.

For complex-valued data (e.g., Transfer Functions), the cost function is computed using the magnitude of residuals or the Hermitian form in the case of GLS, ensuring correct handling of real and imaginary parts.

Parameters:

n_walkers (int, optional) – Number of MCMC walkers. Default is 32. Must be at least 2 * ndim, where ndim is the number of free (non-fixed) parameters, as required by emcee’s ensemble sampler.
n_steps (int, optional) – Number of MCMC steps per walker. Default is 3000.
burn_in (int, optional) – Number of initial steps to discard. Default is 500.
progress (bool, optional) – Whether to show progress bar. Default is True.

Returns:

sampler – The emcee sampler object containing the full chain.

Return type:

emcee.EnsembleSampler

Raises:

ValueError – If there are no free parameters, or if n_walkers is below 2 * ndim.

Notes

Fixed Covariance Assumption

This implementation assumes that the covariance matrix Σ (via cov_inv) is fixed (parameter-independent). Under this assumption, the log determinant term log|Σ| is constant and can be omitted from the log likelihood:

\[\log p(y|\theta) = -\frac{1}{2} r^T \Sigma^{-1} r + \text{const}\]

If Σ depends on the model parameters θ, the full log likelihood including the log|Σ| term must be used:

\[\log p(y|\theta) = -\frac{1}{2} r^T \Sigma^{-1} r - \frac{1}{2} \log|\Sigma| - \frac{N}{2}\log 2\pi\]

This assumption has been validated through unit tests and independent technical review. All models agreed that the current implementation is correct for fixed-covariance use cases.

Complex Residuals

For complex-valued data, the Hermitian form r.conj() @ cov_inv @ r is used, which assumes a circular complex Gaussian distribution (i.e., real and imaginary parts have equal variance and are uncorrelated).

References

class gwexpy.fitting.Fitter(model: Any)#

Bases: object

Wrap fit_series in a small stateful helper class.

fit(series: Any, **kwargs: Any) → FitResult#: Fit the provided series to the model.

class gwexpy.fitting.GeneralizedLeastSquares(x: ndarray, y: ndarray, cov_inv: ndarray, model: Callable[[...], Any], cov: ndarray | None = None)#

Bases: object

Generalized Least Squares (GLS) cost function.

Minimizes χ² = r.T @ cov_inv @ r where r = y - model(x, **params).

This cost function accounts for correlations between data points through the inverse covariance matrix.

Parameters:

x (array-like) – Independent variable (e.g., frequency array).
y (array-like) – Observed data (real-valued).
cov_inv (ndarray) – Inverse covariance matrix, shape (n, n) where n = len(y). Can be obtained from BifrequencyMap.inverse().value.
model (callable) – Model function with signature model(x, *params) -> y. The first argument must be x, followed by fit parameters.
cov (ndarray, optional) – Original covariance matrix. If provided, Cholesky decomposition is used for better numerical stability.

Notes

errordef is set to Minuit.LEAST_SQUARES (= 1.0) for iminuit.

Examples

>>> def linear(x, a, b):
...     return a * x + b
>>> gls = GeneralizedLeastSquares(x, y, cov_inv, linear)
>>> m = Minuit(gls, a=1, b=0)
>>> m.migrad()

errordef = 1.0#

property ndata: int#: Number of data points.

class gwexpy.fitting.GLS(X: ndarray, y: ndarray, cov: ndarray | None = None, cov_inv: ndarray | None = None)#

Bases: object

Direct solver for Generalized Least Squares problems (Linear).

Parameters:

X (array-like) – Design matrix (n_samples, n_params).
y (array-like) – Observation vector (n_samples,).
cov (array-like, optional) – Covariance matrix (n_samples, n_samples).
cov_inv (array-like, optional) – Inverse covariance matrix (n_samples, n_samples).

solve() → ndarray#

Solve the linear GLS problem.

beta = (X.T @ W @ X)^-1 @ X.T @ W @ y where W = cov_inv.

gwexpy.fitting.fit_bootstrap_spectrum(data_or_spectrogram: TimeSeries | Spectrogram, model_fn: Callable, freq_range: tuple[float, float] | None = None, method: str = 'median', rebin_width: float | None = None, block_size: float | str | None = None, ci: float = 0.68, window: str = 'hann', fftlength=None, overlap=None, nfft: int | None = None, noverlap: int | None = None, n_boot: int = 1000, initial_params: dict[str, float] | None = None, bounds: dict[str, tuple[float, float]] | None = None, fixed: list | None = None, run_mcmc: bool = False, mcmc_walkers: int = 32, mcmc_steps: int = 5000, mcmc_burn_in: int = 500, plot: bool = True, progress: bool = True, **kwargs) → FitResult#

Integrated spectral analysis pipeline with bootstrap, GLS fitting, and MCMC.

This function provides a unified workflow for: 1. Converting TimeSeries to Spectrogram (if needed) 2. Bootstrap resampling to estimate PSD and covariance 3. GLS fitting with proper frequency correlation 4. Optional MCMC for Bayesian parameter inference 5. Visualization of results

Parameters:

data_or_spectrogram (TimeSeries or Spectrogram) – Input data. If TimeSeries, a spectrogram will be computed automatically.
model_fn (callable) – Model function with signature model(f, *params) -> y. First argument must be frequency array. Example: lambda f, A, alpha: A * f**alpha
freq_range (tuple of (fmin, fmax), optional) – Frequency range for fitting. If None, use all frequencies.
method (str, optional) – Bootstrap averaging method: ‘median’ (default) or ‘mean’.
rebin_width (float, optional) – Frequency rebinning width in Hz. If None, no rebinning.
block_size (float, Quantity, or 'auto', optional) – Duration of blocks for block bootstrap in seconds. Can be specified as float (seconds), Quantity with time units, or ‘auto’. If None, standard bootstrap.
ci (float, optional) – Confidence interval for bootstrap errors. Default is 0.68 (1-sigma).
window (str, optional) – Window function for spectrogram and correlation correction. Default is ‘hann’.
fftlength (float or Quantity, optional) – FFT segment length in seconds (e.g. 1.0 or 1.0 * u.s). Used to generate the spectrogram from a TimeSeries and for VIF overlap-correlation correction. If None and a TimeSeries is given, GWpy chooses a default covering the full duration.
overlap (float or Quantity, optional) – Overlap between FFT segments in seconds. If None, defaults to the recommended overlap for window (50 % for Hann). Cannot be used with noverlap.
nfft (int, optional) – FFT segment length in samples. Alternative to fftlength. Cannot be used with fftlength.
noverlap (int, optional) – Overlap length in samples. Must be used with nfft. Cannot be used with overlap.
n_boot (int, optional) – Number of bootstrap resamples. Default is 1000.
initial_params (dict, optional) – Initial parameter values for fitting, e.g., {“A”: 10, “alpha”: -1.5}.
bounds (dict, optional) – Parameter bounds, e.g., {“A”: (0, 100), “alpha”: (-5, 0)}.
fixed (list, optional) – List of parameter names to fix during fitting.
run_mcmc (bool, optional) – Whether to run MCMC after fitting. Default is False.
mcmc_walkers (int, optional) – Number of MCMC walkers. Default is 32.
mcmc_steps (int, optional) – Number of MCMC steps. Default is 5000.
mcmc_burn_in (int, optional) – MCMC burn-in steps to discard. Default is 500.
plot (bool, optional) – Whether to display plots. Default is True.
progress (bool, optional) – Whether to show progress bars for MCMC. Default is True.
**kwargs – Additional keyword arguments forwarded to lower-level fitting steps.

Returns:

Fit result object containing: - Best-fit parameters and errors - Chi-square and reduced chi-square - Covariance matrix (in cov_inv and accessible cov attribute) - MCMC samples and intervals (if run_mcmc=True) - Plotting methods

Return type:

FitResult

Examples

>>> from gwexpy.fitting.highlevel import fit_bootstrap_spectrum
>>>
>>> # Define model
>>> def power_law(f, A, alpha):
...     return A * f**alpha
>>>
>>> # Run pipeline
>>> result = fit_bootstrap_spectrum(
...     data,
...     model_fn=power_law,
...     freq_range=(5, 50),
...     fftlength=1.0,
...     overlap=0.5,
...     rebin_width=0.25,
...     block_size=2.0,  # 2 seconds
...     initial_params={"A": 10, "alpha": -1.5},
...     run_mcmc=True,
... )
>>>
>>> # Access results
>>> print(result.params)
>>> print(result.parameter_intervals)  # If MCMC was run

Notes

The covariance matrix from bootstrap is stored and used for GLS fitting, properly accounting for frequency correlations in the spectral estimate.