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).
- パラメータ:
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().
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.
- 戻り値:
Object containing fit results, parameters, and plotting methods.
- 戻り値の型:
- 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)#
ベースクラス:
objectStore 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.
- 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).
- 戻り値:
Shape (n_steps, n_walkers, n_params). Returns None if MCMC not run.
- 戻り値の型:
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.
- 戻り値:
Dictionary mapping parameter names to (lower, median, upper) tuples.
- 戻り値の型:
- 例外:
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.
- パラメータ:
- 戻り値:
figure -- The corner plot 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.
- パラメータ:
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().
- 戻り値:
ax -- The axes with the plot.
- 戻り値の型:
- 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.
- パラメータ:
n_walkers (int, optional) -- Number of MCMC walkers. Default is 32. Must be at least
2 * ndim, wherendimis 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.
- 戻り値:
sampler -- The emcee sampler object containing the full chain.
- 戻り値の型:
emcee.EnsembleSampler
- 例外:
ValueError -- If there are no free parameters, or if
n_walkersis below2 * ndim.
メモ
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 @ ris used, which assumes a circular complex Gaussian distribution (i.e., real and imaginary parts have equal variance and are uncorrelated).参照
- class gwexpy.fitting.Fitter(model: Any)#
ベースクラス:
objectWrap fit_series in a small stateful helper class.
- class gwexpy.fitting.GeneralizedLeastSquares(x: ndarray, y: ndarray, cov_inv: ndarray, model: Callable[[...], Any], cov: ndarray | None = None)#
ベースクラス:
objectGeneralized 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.
- パラメータ:
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.
メモ
errordef is set to Minuit.LEAST_SQUARES (= 1.0) for iminuit.
サンプル
>>> 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#
- class gwexpy.fitting.GLS(X: ndarray, y: ndarray, cov: ndarray | None = None, cov_inv: ndarray | None = None)#
ベースクラス:
objectDirect solver for Generalized Least Squares problems (Linear).
- パラメータ:
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).
- 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
- パラメータ:
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.0or1.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.
- 戻り値:
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
- 戻り値の型:
サンプル
>>> 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
メモ
The covariance matrix from bootstrap is stored and used for GLS fitting, properly accounting for frequency correlations in the spectral estimate.
参考
gwexpy.spectral.bootstrap_spectrogramBootstrap resampling function
gwexpy.fitting.fit_seriesLower-level fitting function
gwexpy.fitting.GeneralizedLeastSquaresGLS cost function