Fitting engine

Asymmetry fits μSR data with an iminuit-based engine that exposes the same machinery through the GUI fit panels and through the Python API. The minimiser is MINUIT [1]: the engine minimises the least-squares cost \(\chi^2 = \sum_i [(A_i - f(t_i))/\sigma_i]^2\) with migrad, and takes its symmetric parameter uncertainties from the HESSE step that migrad runs on convergence. This chapter documents the engine itself — the workflow, the parameter and results objects, the statistics, the global-fit API, and the conventions that apply across every model. The physical-model catalogue is split out into dedicated chapters: see Composite models for the expression grammar that drives the GUI builder, Fit wizard and Global fit wizard for the guided model-selection workflows, and the per-function chapters for the depolarisation forms themselves (Relaxation, Oscillation, Muonium, Kubo–Toyabe, Nuclear dipolar, Frequency domain, Field-dependent transport models (DiffusionLF / BallisticLF), Superconducting penetration depth models).

When fitting through the GUI, the plot panel’s current bunch factor is applied before the dataset is handed to the fitter. If the plot is showing a rebinned dataset, the next GUI fit will use that rebinned dataset. To fit the full-resolution data and only simplify the view afterwards, run the fit with Bunch = 1 and then increase the bunch factor; the fit curve remains overlaid on the plot.

Basic workflow

The minimal end-to-end fit is a model from the MODELS registry, a ParameterSet of initial values and bounds, and a FitEngine call:

from asymmetry.core.fitting.engine import FitEngine
from asymmetry.core.fitting.models import MODELS
from asymmetry.core.fitting.parameters import Parameter, ParameterSet
from asymmetry.core.io import load

dataset = load("data.nxs")

model = MODELS["ExponentialRelaxation"]
params = ParameterSet([
    Parameter("A0", value=25.0, min=0.0),
    Parameter("Lambda", value=0.5, min=0.0),
    Parameter("baseline", value=0.0),
])

result = FitEngine().fit(dataset, model.function, params)

print(f"χ²ᵣ = {result.reduced_chi_squared:.3f}")
for p in result.parameters:
    err = result.uncertainties.get(p.name, 0.0)
    print(f"  {p.name} = {p.value:.6f} ± {err:.6f}")

For anything richer than a single-component model — multiple amplitudes, shared baselines across detector groups, multiplicative envelopes, fraction-constrained additive groups — use the composite-model path in Composite models. The MODELS registry is intentionally short and contains only the standalone variants with explicit baselines; composite expressions are the canonical way to assemble realistic muSR models.

Parameter control

Bounds, fixing, and initial values are set on individual Parameter objects:

from asymmetry.core.fitting.parameters import Parameter

p = Parameter(name="Lambda", value=0.5, min=0.0, max=10.0)

# Hold a parameter constant during the fit
p_fixed = Parameter(name="A0", value=25.0, fixed=True)

# Or toggle later
p_fixed.fixed = True

In the GUI parameter table, the Fix checkbox does the same job; bounds are entered directly in the min and max columns. Parameter metadata — symbol, unit, default lower bound, physical description — comes from PARAM_INFO_REGISTRY; symbols and units displayed in the GUI, in fit reports, and in exported GLE figures all derive from this registry, so they are guaranteed to be consistent with what the model function actually expects.

Fit statistics

The FitResult exposes the standard goodness-of-fit numbers:

print(f"χ²  = {result.chi_squared:.4f}")
print(f"χ²ᵣ = {result.reduced_chi_squared:.4f}")
print(f"converged: {result.success}")
print(f"message:   {result.message}")

A reduced \(\chi^2\) near 1 is the usual target. Values much greater than 1 indicate that either the model is missing physics, the parameter errors on the input asymmetry are underestimated, or both. Values much less than 1 generally mean the input errors are overestimated. Neither case is automatically fatal — in the first, the next move is usually to add or substitute a component motivated by the physics; in the second, the parameter uncertainties on the fitted values will be inflated correspondingly and should be re-derived against a calibrated error model.

Advisory warnings (result.warnings)

Before it minimises, the engine runs two advisory guards that flag the setup traps most likely to produce a converged-but-wrong fit, and carries their messages on result.warnings (they are also emitted through the Python warnings system). The guards never raise and never change the outcome:

  • AsymmetryScaleWarning — the data and the seeded model sit on different asymmetry scales (the percent-vs-fraction mix-up; a loaded MuonDataset.asymmetry is on the percent scale).

  • FixedFrequencyFieldMismatchWarning — a precession frequency held fixed more than ~2 % from \(\gamma_\mu B\) implied by the run’s field; the misfit then leaks into the damping term and inflates the fitted Gaussian sigma (the vortex-state transverse-field trap).

The GUI fit panel surfaces these beneath the converged summary (deduplicated for a batch/global fit). The science and the fix for each are in the cookbook — see the transverse-field frequency entry and the asymmetry-scale entry.

A stuck χ² ≈ 200? Free the background amplitude (TF data)

A frequent cause of a huge, stubborn reduced \(\chi^2\) (often \(\chi^2_r \approx 200\)) on muSR transverse-field data is a background amplitude pinned at a non-physical lower bound. The constant background term A_bg (and the standalone baseline) defaults to a zero lower bound, but for TF data the fitted background is routinely negative — commonly around \(-20\%\) of the full asymmetry — because of detector-pair imbalance and the asymmetry baseline convention. Clamped at zero, the minimiser cannot reach the true minimum and the fit parks at a large \(\chi^2\).

The fix is to initialise the background unbounded (or with an explicitly negative lower bound) so it can go negative:

from asymmetry.core.fitting.parameters import Parameter

# WRONG for TF data: a zero floor traps A_bg → χ²ᵣ stuck near 200
bad = Parameter("A_bg", value=0.0, min=0.0)

# RIGHT: leave it unbounded so the background can settle negative
bg = Parameter("A_bg", value=0.0)              # min=-inf, max=+inf by default
# ...or pin a generous negative floor if you want some guard rails:
bg = Parameter("A_bg", value=-0.05, min=-0.5, max=0.5)

In the GUI parameter table, clear the min cell for the background row (or set it negative) before fitting TF data. If a TF fit converges with a glued-to-zero background and a \(\chi^2_r\) in the hundreds, this bound is the first thing to check.

Parameter values and Hessian uncertainties come from result.parameters and result.uncertainties:

for p in result.parameters:
    err = result.uncertainties.get(p.name, 0.0)
    print(f"{p.name}: {p.value:.6f} ± {err:.6f}")

Hessian (HESSE) errors assume the local quadratic approximation around the minimum is a good description of the likelihood, i.e. that the parameter uncertainties are symmetric. For well-conditioned fits this holds. For ill-conditioned problems — strongly correlated parameters, fits near a bound, multi-modal likelihoods — the interval is skewed and the symmetric error misrepresents it; use MINUIT’s asymmetric MINOS intervals [1] (Asymmetric (MINOS) errors) or the re-simulation estimate below (Monte Carlo error estimation).

Residuals

Inspecting the residual time series is the most reliable check that the model is structurally appropriate:

import matplotlib.pyplot as plt

fit_values = model.function(
    dataset.time, **{p.name: p.value for p in result.parameters}
)
residuals = (dataset.asymmetry - fit_values) / dataset.error

plt.plot(dataset.time, residuals, "o", alpha=0.6)
plt.axhline(0, color="k", linestyle="--")
plt.xlabel("Time (μs)")
plt.ylabel("Normalised residual (σ)")

A good fit produces residuals that scatter randomly around zero with most points within \(\pm 2\sigma\) and no systematic trend over time. Persistent structure — a slow oscillation, a turn-up at long times, a visible knee — almost always points at a missing component rather than at noise.

Restricting the fit range

The t_min and t_max keyword arguments to FitEngine.fit restrict which time bins enter the cost function:

result = engine.fit(
    dataset, model.function, params, t_min=0.1, t_max=5.0,
)

The two common reasons to clip the range are (a) early-time bins contaminated by deadtime or finite-pulse-width effects, and (b) late-time bins where the signal has decayed into the background and adds only noise to the fit. The same range is honoured by the fit wizard and the global-fit workflow.

Global fitting

The global-fit interface accepts a list of datasets, the model function, and explicit lists of parameter names to treat as global (shared across all runs) or local (independent per run):

from asymmetry.core.fitting.engine import FitEngine
from asymmetry.core.fitting.models import MODELS
from asymmetry.core.fitting.parameters import Parameter, ParameterSet
from asymmetry.core.io import load

datasets = [load(f"run_{i}.nxs") for i in range(1, 4)]
model = MODELS["ExponentialRelaxation"]
engine = FitEngine()

initial_params = {}
for ds in datasets:
    ps = ParameterSet([
        Parameter("A0", value=25.0, min=0.0),
        Parameter("Lambda", value=0.5, min=0.0),
        Parameter("baseline", value=0.0),
    ])
    initial_params[ds.run_number] = ps

results_dict, global_result = engine.global_fit(
    datasets=datasets,
    model_fn=model.function,
    global_params=["A0"],
    local_params=["Lambda", "baseline"],
    initial_params=initial_params,
    t_min=0.1,
    t_max=10.0,
)

for p in global_result:
    print(f"[global] {p.name} = {p.value:.6f}")
for run_number, run_result in results_dict.items():
    print(f"run {run_number}: χ²ᵣ = {run_result.reduced_chi_squared:.3f}")

The GUI Batch tab automates the same workflow: select multiple datasets, mark parameters as Global (shared across runs) or Local (per-run) in the parameter table, and click Run Batch Fit. A fit where at least one parameter is Global is a global fit; otherwise each run is fitted independently but the results are collected into one trendable series. Results land in the Fit Parameters panel where they can be browsed, exported to TSV, or passed into the parameter-trending fit framework documented in Parameter trending.

The same engine is also used in the Frequency workspace for displayed Fourier spectra. In that mode the Fit dock switches labels from A(t) to S(ν), uses MHz fit ranges, and offers peak/background components documented in Frequency-domain fitting.

The Global fit wizard automates model selection on ordered field or temperature series; use that wizard before constructing a hand-built global fit if you do not yet know which composite model the data prefer.

Monte Carlo error estimation

For parameters whose Hessian uncertainties are not trustworthy — heavy correlations, fits sitting near a bound, multi-modal likelihoods — re-simulating the data and refitting gives an empirical distribution. The snippet below draws each replica by adding Gaussian noise at the per-bin error to the fitted data (a parametric Monte Carlo), so the spread it reports is only as good as that error model: it inherits the Gaussian-error assumption and will not capture non-Gaussian tails, and with a small replica count the estimated width itself carries appreciable scatter (raise n_bootstrap and seed the generator for a reproducible number). It is an approximation to the error, not an exact interval — for the definitive error-bar check re-simulate at matched statistics and inspect the pull distribution (Validating error bars: the pull distribution).

import numpy as np
from asymmetry.core.data import MuonDataset

rng = np.random.default_rng(0)
n_bootstrap = 500
bootstrap_results = []
for _ in range(n_bootstrap):
    resampled = dataset.asymmetry + rng.standard_normal(len(dataset.asymmetry)) * dataset.error
    boot_dataset = MuonDataset(
        time=dataset.time, asymmetry=resampled, error=dataset.error,
    )
    bootstrap_results.append(engine.fit(boot_dataset, model.function, params))

lambda_dist = np.array([
    next(p.value for p in r.parameters if p.name == "Lambda")
    for r in bootstrap_results
])
print(f"Lambda = {lambda_dist.mean():.4f} ± {lambda_dist.std():.4f}")

For a single parameter whose likelihood is skewed, the engine’s own asymmetric MINOS interval (FitEngine.fit(..., minos=True), reported in result.minos_errors) is usually the cheaper and more direct check; see Asymmetric (MINOS) errors. The definitive test of whether the reported uncertainties are calibrated is the pull diagnostic (Validating error bars: the pull distribution).

Custom models

Any callable with signature f(t, **params) -> array is a valid model function:

import numpy as np

def my_model(t, amplitude, tau, frequency):
    return amplitude * np.exp(-t / tau) * np.sin(2 * np.pi * frequency * t)

params = ParameterSet([
    Parameter("amplitude", value=20.0),
    Parameter("tau", value=1.0, min=0.0),
    Parameter("frequency", value=5.0, min=0.0),
])

result = engine.fit(dataset, my_model, params)

Custom callables do not get the symbol/unit metadata that registered components carry, so the GUI parameter table will show raw names and no units. For anything used more than once, a proper component or ModelDefinition registration is preferable.

References

[1] F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).

See also