Fitting

fit(problem::FitProblem; cost=:auto, maxiters=500, tol=nothing,
    scale_covariance=:auto, initial_guesses=nothing, multistart=1,
    solver=nothing) -> FitResult

Fit a validated Gaussian FitProblem.

Keyword contracts:

• cost: :auto, :chi2, or :gaussian_likelihood. :auto selects :chi2 for a static covariance, and the normalized Gaussian -2 log(L) cost when the effective covariance depends on the parameters (x uncertainties, or error components with mode=:model_relative).
• maxiters: positive solver budget used for every candidate.
• tol: positive stopping tolerance; nothing selects default_fit_tolerance for the solver and derivative mode.
• scale_covariance: :auto, :never, or :always. :auto multiplies the parameter covariance by chi2/ndf, and only when neither y nor x observation uncertainty was supplied; :never reports the unscaled covariance and :always forces the factor. Profile intervals and the stationarity check use the same scale, so they stay consistent with param_stderr. :always is rejected for cost=:gaussian_likelihood, where the covariance comes from the cost Hessian.
• initial_guesses: additional complete parameter vectors in p0 order. Explicit guesses are always tried; the candidate budget grows to cover them.
• multistart: total candidate budget, including problem.p0 and every explicit guess. Values beyond that add generated candidates. Generated candidates are deterministic: interval points when finite bounds exist, otherwise fixed rescalings of p0. No random number generation is involved.
• solver: nothing keeps the automatic choice: LsqFit exactly when static chi-square least squares represents the complete problem, the scalar solver otherwise. Pass OptimizationSolver(algorithm), NativeMinuitSolver(), or one of the shorthands :lbfgs, :ipnewton, :nelder_mead for an explicit choice. The same solver and settings are used for multistart and profile refits; result.backend records what actually solved the fit.

The candidate with the lowest finite cost is returned; convergence status only breaks exact cost ties. A non-converged winner is returned with converged == false even when another candidate converged at a higher cost — run diagnostic_dashboard(result) before interpreting such a result. If no candidate produces a finite result, the last model, validation, or solver error is rethrown.

fit_model(model, x, y; p0, kwargs...) -> FitResult

Fit scalar observations y measured at x with a model satisfying model(x, p) -> yhat. x, y, and p0 are copied to finite Float64 storage; length(x) must equal length(y) and the model must return one finite prediction per observation.

Observation uncertainty: supply at most one of sigma_y/cov_y and at most one of sigma_x/cov_x; y- and x-uncertainties and named error_components may be combined. A whitening operator supplies the complete covariance and excludes every other observation-uncertainty keyword. With no uncertainty, unweighted least squares is used.

For fits with x uncertainty, x_derivative=(x, p) -> dy_dx supplies a vectorized model derivative with respect to x. This avoids the default point-by-point AD path and is the preferred route for large datasets.

For large static correlated datasets, whitening=WhiteningOperator(...) supplies the complete covariance through a matrix-free operation. The operator must accept generic AbstractVector inputs and AD element types when the general optimizer is used.

An optional analytic jacobian(x, p) receives the full parameter vector (fixed values filled in) and returns the length(x) × length(p0) matrix J[i, j] = ∂ŷᵢ/∂pⱼ, including columns for fixed parameters; the free columns are selected internally. Prediction bands evaluate it on their own x grid, so it must accept any coordinate vector.

Set inplace=true for model!(out, x, p). The unbounded least-squares backend uses LsqFit's native in-place model interface; generic optimizer paths preserve the same contract with a type-correct output buffer. With inplace=true, an analytic Jacobian instead uses jacobian!(J, x, p) with the same shape.

Parameter control uses bounds (a (lower, upper) tuple of vectors with one entry per parameter, -Inf/Inf for unbounded sides), constraints, parameter_priors, parameter_constraints, and fixed_parameters; the accepted forms for the latter four are documented under FitProblem. Solver keywords are forwarded to fit(::FitProblem) with defaults cost=:auto, maxiters=500, tol=nothing, scale_covariance=:auto, initial_guesses=nothing, multistart=1, and solver=nothing. Omitting tol selects the solver-specific default_fit_tolerance.

Use derivatives=:finite for models implemented outside Julia or restricted to ordinary floating-point inputs. The policy also controls covariance, profiles, and predictions; see FitProblem for its numerical assumptions. With LsqFit/Optimization this changes the default tolerance to 1e-6 to account for differenced-gradient noise. NativeMinuit retains its native EDM (estimated distance to minimum) stopping criterion. An explicitly supplied tolerance is never relaxed.

Returns a FitResult. Invalid inputs raise an error before a result is constructed: DimensionMismatch for length mismatches (x vs. y, sigma vectors, bounds), DomainError for out-of-range values (non-positive standard deviations, reversed bounds), and ArgumentError for all remaining invalid inputs (non-finite values, contradictory uncertainty inputs, invalid covariance matrices, malformed parameter controls).

Example

model(x, p) = @. p[1] * x + p[2]
result = fit_model(model, [0.0, 1.0, 2.0], [0.1, 1.2, 1.9];
    p0=[1.0, 0.0], sigma_y=fill(0.2, 3))

FitProblem(model, x, y; p0, sigma_y, sigma_x, cov_y, cov_x, whitening,
           error_components, bounds, constraints, parameter_priors,
           parameter_constraints, fixed_parameters, jacobian,
           x_derivative, inplace=false, derivatives=:auto)

Build a fit problem for 1D x and scalar y observations.

model(x, p) receives the full x vector and the complete parameter vector p (ordered as in p0) and returns one prediction per observation, for example model(x, p) = @. p[1] * x + p[2]. p0 is the vector of initial parameter values; its ordering defines the parameter indices used by every other keyword.

sigma_y and sigma_x are per-observation 1-sigma standard deviations (length n, entries > 0). cov_y and cov_x are full n x n symmetric positive-definite observation covariance matrices. sigma_y and cov_y are mutually exclusive, as are sigma_x and cov_x.

bounds = (lower, upper) takes two vectors of length length(p0); use -Inf/Inf entries for one-sided or absent bounds.

jacobian(x, p) optionally supplies the analytic length(x) x length(p0) matrix of derivatives dmodel_i/dp_j; when omitted, the Jacobian follows the derivatives mode. See inplace below for the mutating signatures.

constraints accepts either ConstraintSpec or a NamedTuple with:

• ineq = p -> vector interpreted as ineq(p) <= 0
• eq = p -> vector interpreted as eq(p) == 0

Both callbacks receive the complete parameter vector in p0 order, including entries removed from optimization by fixed_parameters.

parameter_priors adds Gaussian penalties in parameter space and accepts a single NamedTuple or vector of NamedTuples:

• (index=i, mean=mu, sigma=sigma)
• (index=i, mean=mu, sigma_minus=sminus, sigma_plus=splus)

parameter_constraints adds correlated Gaussian constraints:

• (indices=[i, j], mean=[mu_i, mu_j], covariance=cov)

error_components adds named y/x uncertainty sources:

• (name=:stat, target=:y, mode=:absolute, values=sigma)
• (name=:scale, target=:y, mode=:relative, values=0.02)
• (name=:model_scale, target=:y, mode=:model_relative, values=0.02)
• (name=:corr, target=:y, mode=:covariance, values=cov)

fixed_parameters removes parameters from the optimizer and accepts:

• i => value
• (index=i, value=value)
• (index=i, value=value, sigma=sigma)
• (index=i, value=value, sigma_minus=sminus, sigma_plus=splus)

x_derivative(x, p) optionally supplies the vector derivative dy/dx used for effective x-uncertainty propagation. If omitted, ScientificFitting differentiates the model with respect to each x value using the selected differentiation mode.

derivatives=:auto preserves the native solver defaults and ForwardDiff for post-fit derivatives. Use :finite for Float64-only models (including Python callbacks). This selects numerical differentiation throughout fitting, covariance estimation, profiles, and prediction bands. Analytic jacobian and x_derivative callbacks take precedence where applicable. Numerical differences require a smooth model in a neighborhood of the evaluation point, including at parameter bounds; they do not make discontinuous objectives differentiable.

whitening=WhiteningOperator(...) supplies the complete static observation covariance through a matrix-free whitening operation. It is mutually exclusive with y/x uncertainties and active error_components; combining covariance models without an explicit derivation would change the statistical model.

Set inplace=true when the model has the signature model!(out, x, p). On the unbounded least-squares path, ScientificFitting forwards this contract to LsqFit's native in-place solver interface. Other solver paths use the same model through a type-preserving output buffer. If jacobian is also supplied, its in-place signature must be jacobian!(J, x, p). Mutating functions used with bounds, constraints, or parameter-dependent covariance must accept buffers whose element type is chosen by automatic differentiation; avoid Float64-specific method signatures.

FitOptions(; cost=:auto, maxiters=500, tol=1e-10, scale_covariance=:auto,
            multistart=1, parameter_covariance=:auto, solver=nothing)

Normalized solver and covariance options stored in a FitResult. User-facing fit functions expose these as keyword arguments; constructing FitOptions directly is mainly useful for lower-level workflows and tests. See fit for the meaning of each option; allowed values are cost in (:auto, :chi2, :gaussian_likelihood), scale_covariance in (:auto, :always, :never), parameter_covariance in (:auto, :hessian, :none), and multistart >= 1 candidate starting points. Invalid iteration, tolerance, covariance-scaling, and multistart settings fail during construction rather than inside a solver. solver is either nothing (the automatic choice) or a concrete AbstractFitSolver; the resolved choice is preserved by profile refits, and the backend that actually solved a fit is recorded on the result. Unlike the public fit keywords, tol here is a resolved positive number, not nothing; see default_fit_tolerance.

fit(problem::LikelihoodFitProblem; maxiters=1000, tol=nothing,
    initial_guesses=nothing, multistart=1,
    parameter_covariance=:auto, solver=nothing) -> LikelihoodFitResult

Minimize a validated likelihood-scale or custom objective problem. Bounds, fixed parameters, Gaussian parameter terms, nonlinear constraints, observation count, and cost name are already stored in problem.

initial_guesses may contain additional complete parameter vectors in p0 order. Explicit guesses are always tried; the candidate budget grows to cover them. multistart is the total candidate budget including p0 and every explicit guess. Remaining slots are filled with deterministic generated candidates: with declared bounds, the midpoints of the finite intervals plus interior points at 25% and 75% of each two-sided finite interval (parameters without such an interval keep their p0 value); without declared bounds, the rescalings 0.5*p0, 2*p0, and -p0 of the free parameters. A candidate is usable only when it yields a finite cost (for likelihood objectives: lies inside the distribution's support). ScientificFitting returns the candidate with the lowest finite cost; convergence status only breaks exact cost ties, and a non-converged winner is returned with converged == false. If no candidate produces a finite result, the last objective, validation, or solver error is raised.

The solver choice defaults to LBFGS, or IPNewton with nonlinear constraints. Pass solver=OptimizationSolver(algorithm) or solver=NativeMinuitSolver() to choose explicitly; solver objects and their native settings are preserved in profile refits, and every solver declares its capabilities, so incompatible combinations (for example nonlinear constraints with a derivative-free method) raise an error instead of being dropped. The shorthands solver=:lbfgs, solver=:ipnewton, and solver=:nelder_mead name the corresponding OptimizationSolver algorithms. All methods are local optimizers of continuous parameters, not guarantees of a global minimum.

parameter_covariance=:auto chooses :none for derivative-free solvers (such as :nelder_mead) and :hessian otherwise. Explicit :hessian computes 2 * inv(H) for the complete cost; use it only for locally smooth likelihoods on the documented -2 log(L) scale. :none skips curvature and stores NaN for free-parameter errors; fixed parameters keep their declared external uncertainty (zero when none was given), independent of this choice. Profile refits preserve both options. Supply explicit profile ranges when local errors are unavailable; profile thresholds still require statistical justification for non-regular models.

maxiters and explicit tol must be positive. tol=nothing selects default_fit_tolerance for the solver and derivative mode. With Nelder-Mead, maxiters limits objective evaluations, tol sets NLopt's absolute/relative parameter tolerances, and iterations is missing, not an invented count. The absolute tolerance acts in parameter units and one tol covers every parameter, so bring parameters to comparable magnitudes or set tol for the smallest relevant scale; tol is not a statistical error. Function-value stopping is disabled because equal costs need not mean a contracted simplex. A budget-limited solve returns converged == false.

fit_custom(objective; p0, nobs, gof=nothing, cost_name=:custom,
           kwargs...) -> LikelihoodFitResult

Fit a user-defined scalar objective. objective(p) is minimized directly. If gof(p) is supplied it is used for reduced goodness-of-fit statistics and p-values; otherwise these fields are NaN. nobs controls degrees of freedom and the BIC sample-size term, and must represent the number of statistically independent observations used by the objective; Gaussian parameter terms are counted as additional observations in both.

When requested, ScientificFitting computes local covariance as 2 * inv(H), where H is the Hessian of the complete minimized cost (the objective plus any parameter priors and parameter constraints), as for fit(::LikelihoodFitProblem). AIC, BIC, and that covariance therefore have their standard interpretation only when objective is a normalized -2log(L) cost (Gaussian chi-square is on the same scale). For an arbitrary loss, the optimizer result remains usable but those inferential fields are only arithmetic summaries.

Common keywords are bounds, constraints, parameter_priors, parameter_constraints, fixed_parameters, parameter_names, maxiters, tol, initial_guesses, multistart, and derivatives=:auto (or :finite). See ConstraintSpec, ParameterPrior, ParameterConstraint, and FixedParameter for the accepted forms of the corresponding keywords. solver and parameter_covariance independently control minimization and local errors; see fit(::LikelihoodFitProblem) for choices and limits. Parameter callbacks receive the complete vector in p0 order. nobs must be positive. Invalid parameter controls raise DimensionMismatch, DomainError, or ArgumentError depending on the violation; a non-finite objective fails with an error rather than producing a reportable result.

Example

cost(p) = ((p[1] - 2.0) / 0.3)^2
result = fit_custom(cost; p0=[0.0], nobs=1, cost_name=:calibration_chi2)

fit_likelihood_model(model, x, y; logprob=nothing, error=nothing, p0, kwargs...)
    -> LikelihoodFitResult

Fit independent observations with a user-defined continuous or discrete observation distribution. model(x, p) returns one prediction per observation; logprob(y, prediction, p) returns a vector of normalized log densities or log probability masses. The two callbacks are evaluated once per objective evaluation, allowing vectorized Julia or Python models. Capture per-observation scales, trial counts, or other known inputs in the logprob closure.

Alternatively, error=distribution models the additive residual y - model(x,p) directly with a fixed Distributions.jl distribution. A univariate object applies independently to every residual; a vector of objects gives point-specific errors. A multivariate object describes the complete residual vector jointly, including dependence between measurements. Its dimension must equal the number of values. Use exactly one of error and logprob; no Gaussian approximation is introduced. For prediction-dependent distributions such as Poisson counts, keep logprob; these describe observations, not additive continuous errors.

The core minimizes -2 * sum(logprob(...)). Include all normalization terms, especially those depending on fitted parameters. A log density may be positive; -Inf denotes zero probability and yields infinite cost, without clipping. NaN or positive-infinity log probabilities, non-finite predictions, and non-finite observations raise ArgumentError; wrong input or output dimensions raise DimensionMismatch. Begin at finite cost inside the distribution's support. The package cannot verify that an arbitrary callback is normalized.

The default gradient/Hessian backend requires a smooth objective near evaluated parameters, even for discrete observations. Choose solver=:nelder_mead for objectives that are not smooth or whose support boundary moves with the parameters; it does not compute Hessian errors unless parameter_covariance=:hessian is explicitly requested. ScientificFitting fits continuous parameters only. Use fit_custom for dependent observations with a joint likelihood rather than multiplying their marginal probabilities.

Parameter controls, derivatives, and solver options follow fit_custom. Requested local covariance and profiles use the same complete likelihood. There is no universal chi-square goodness statistic: chi2, chi2_ndf, and pvalue are NaN unless a justified gof(p) is supplied. Profile coverage is asymptotic, not automatically guaranteed for every distribution or parameter boundary.

Example

using Distributions
model(x, p) = p[1] .* x .+ p[2]
sigma = [0.2, 0.3, 0.2, 0.4]
# Student-t errors allow heavier tails; sigma is a scale, not a standard deviation.
logprob(y, mu, p) = logpdf.(TDist(4), (y .- mu) ./ sigma) .- log.(sigma)
result = fit_likelihood_model(model, [0., 1., 2., 3.], [0.1, 1.2, 1.9, 3.4];
    logprob=logprob, p0=[1., 0.])

fit_distribution(make_distribution, data; p0, obsdim=nothing, kwargs...)
    -> LikelihoodFitResult

Fit independent events using make_distribution(p), which returns a normalized Distributions.jl distribution. The factory runs once per objective evaluation, not once per event: expensive parameter-dependent normalization belongs there. ScientificFitting uses the upstream loglikelihood/logpdf implementation and minimizes -2 log(L). Mixtures whose components have different types (heterogeneous) batch the upstream component log densities in bounded event blocks, so temporary AD arrays do not grow with the sample size; mixtures of a single component type (homogeneous) keep their specialized native loop. Every event contributes; there is no subsampling. When a free mixture weight reaches exactly zero, first and second derivatives with respect to that weight are retained. A distribution exposing only pdf uses log(pdf); extreme-tail underflow then follows that upstream density implementation. Univariate PDF-only components also work inside native mixtures and products: an internal adapter supplies logpdf without changing the upstream model. No parameters are inferred from a distribution instance; declare which are fitted explicitly in the factory and p0.

For univariate continuous or discrete observations, pass a vector. For joint multivariate events pass a matrix and specify obsdim=1 (rows are events) or obsdim=2 (columns are events). Events are independent, but components of each event may be correlated. The observation count used for degrees of freedom and the BIC sample-size term counts events, not scalar coordinates of a multivariate event. Discrete observations do not enable discrete fitted parameters. Data are copied once.

Distribution support, normalization and derivatives come from the upstream package. Zero probability yields infinite cost, without clipping. Start inside the support; use derivatives=:finite if a constructor does not support dual numbers, or a derivative-free solver for nonsmooth parameter dependence. Bounds, parameter terms, solvers, multistart and covariance options follow fit_custom. Unbinned maximum likelihood carries no universal goodness-of-fit statistic, so the result's goodness-of-fit fields (chi2, chi2_ndf, pvalue) are NaN.

After using DistributionsHEP, factories returning ExtendedMixtureModel are also supported. These retain the extended Poisson yield term rather than being treated as just a normalized mixture. Component distributions must describe the observed domain (use upstream truncation for a selected interval).

With using BuildConstructors, an AbstractConstructor can replace the factory. Its names, starts, bounds and fixed/shared state are read through the public metadata API; conflicting shared descriptors are rejected. An optional named p0=(mu=0.2,) overrides free starts or supplies missing ones. Names, bounds and fixed state have one source of truth: set them on the constructor, not through duplicate fit keywords. Numerical metadata uncertainties are not priors or reported fixed-parameter errors. The fit holds a private constructor snapshot; fitted_model(result) reconstructs the fitted model without changing the original constructor. BuildConstructors.parameter_values(result) returns a named tuple for an explicit update! if desired. No Minuit dependency is needed.

Example

using Distributions
result = fit_distribution(p -> Normal(p[1], exp(p[2])), [-0.5, 0.2, 0.8, 1.1];
    p0=[0.0, 0.0], parameter_names=["mean", "log_scale"])

fit_distribution(make_distribution, edges, counts;
    p0, total_count=nothing, integration=:auto, rtol=1e-8, kwargs...)
    -> LikelihoodFitResult

Fit independent Poisson bin counts using a univariate distribution factory. edges must be finite and strictly increasing. Bin i is (edges[i], edges[i+1]]; for integer data, half-integer edges avoid any endpoint convention ambiguity. Counts must be nonnegative integers. The factory runs once per objective call.

For a normalized distribution, total_count is required: it is the known expected event count on the distribution's full support, not an automatically fitted or conditioned observed total. For a fitted rate use an ExtendedMixtureModel from DistributionsHEP and omit total_count; its component yields supply the rates. There is no implicit renormalization over the histogram window. Truncate the upstream distribution explicitly when it describes a selected sample instead.

integration=:auto uses upstream CDF/log-CDF differences where available and adaptive QuadGK integration otherwise. Roundoff-sized CDF range/order violations are reintegrated from the PDF, not clipped; larger violations raise an error. :cdf requires valid CDFs and does not apply this recovery. :quadgk forces quadrature for continuous components. Its estimated error uses relative tolerance rtol in the maximum norm of the value and all carried AD coefficients. Moving truncation bounds retain their derivatives through a fixed-interval map. Mixture components are integrated in batches, with dispatch outside the bin loop and temporary storage linear in the number of bins. Their normalization and parameter derivatives are retained. If an upstream implementation rejects dual numbers, the fit raises an error; pass derivatives=:finite explicitly. Derivatives are never dropped silently.

The objective is normalized Poisson -2log(L), with tail bin probabilities kept in log space. Bins with zero predicted probability contribute zero cost when their observed count is zero and infinite cost otherwise. Goodness-of-fit uses the asymptotic Poisson deviance only when every expected bin count is strictly positive; otherwise chi2, chi2_ndf, and pvalue are NaN. Low counts and fitted boundaries can also invalidate that approximation. Other fit controls and fitted_model work as for unbinned distribution fits. BuildConstructors metadata can be used with either form.

fit_poisson_model(model, x, counts; p0, kwargs...) -> LikelihoodFitResult

Fit count data with a Poisson likelihood. model(x, p) must return the nonnegative expected counts for each observation. counts must contain finite non-negative integer-valued observations; x must be finite and have the same length.

The minimized objective is normalized Poisson -2 log(L). stats.chi2 stores the Poisson deviance and its p-value uses the asymptotic chi-square reference. At an exactly zero expectation this reference is not regular, so the goodness-of-fit fields are NaN. A zero expectation with a positive observation has infinite cost; a zero observation contributes 2mu, without a probability floor. Common parameter-control and solver keywords are listed under fit_custom. Length mismatches raise DimensionMismatch; negative counts or expectations raise DomainError; other invalid input, such as non-finite values or non-integer counts, raises ArgumentError. Returns LikelihoodFitResult.

fit_histogram_model(expected_counts, edges, counts; p0, kwargs...)
    -> LikelihoodFitResult

Fit binned counts. expected_counts(edges, p) must return one nonnegative expected count per bin. Histogram edges must be finite and strictly increasing; counts must contain finite non-negative integer-valued observations.

The objective is normalized Poisson -2 log(L) and stats.chi2 is Poisson deviance. The caller is responsible for integrating any continuous density over the bins; use fit_histogram_density when ScientificFitting should perform that integration. Common parameter-control and solver keywords are listed under fit_custom. Length mismatches raise DimensionMismatch; edges that are not strictly increasing, negative counts, or negative expectations raise DomainError; other invalid input raises ArgumentError. Returns LikelihoodFitResult. Exactly zero expectations follow the same support and goodness-of-fit rules as fit_poisson_model.

fit_histogram_density(pdf, edges, counts; p0, total_count=sum(counts),
                      rtol=1e-8, vectorized=false, kwargs...) -> LikelihoodFitResult

Fit binned counts from a probability density. Expected bin counts are computed with adaptive Gauss-Kronrod quadrature. pdf(x, p) must be normalized on its intended physical domain; ScientificFitting does not renormalize it over the supplied bins. total_count and rtol must be finite and positive.

With vectorized=true, pdf(xs, p) receives a vector of quadrature nodes and must return one value per node. QuadGK batches these evaluations while retaining adaptive error control for each bin; the scalar callback remains the default. During automatic differentiation, the quadrature error norm includes the value, gradient and nested Hessian coefficients, not just the density value itself.

Each expectation is total_count * integral(pdf, edges[i], edges[i+1]). If pdf assigns probability mass outside [first(edges), last(edges)], these expectations sum to less than total_count; with the default total_count=sum(counts) the model then predicts fewer events than observed, and the fit distorts the shape parameters to compensate. Either pass the true full-support event count as total_count, or supply a density normalized on the binned range. Common parameter-control and solver keywords are listed under fit_custom. Quadrature/model failures are propagated. Length mismatches raise DimensionMismatch; non-increasing edges, negative counts, and non-positive total_count or rtol raise DomainError; other invalid input raises ArgumentError. Returns LikelihoodFitResult.

fit_unbinned_model(pdf, data; p0, vectorized=false, kwargs...) -> LikelihoodFitResult

Fit independent unbinned observations with a normalized positive density pdf(x, p). ScientificFitting checks positivity at the observations but cannot infer or verify the normalization domain. Observations must be finite. With vectorized=true, pdf(data, p) evaluates all observations in one call and must return one density per observation. The scalar callback is the default.

The objective is -2 * sum(log(pdf(x_i, p))). No universal chi-square goodness-of-fit statistic exists, so chi2, chi2_ndf, and pvalue are NaN. Common parameter-control and solver keywords are listed under fit_custom. Non-positive densities and non-finite observations raise ArgumentError. Returns LikelihoodFitResult.

fit_extended_unbinned_model(rate, data, domain; p0, rtol=1e-8,
                            vectorized=false, kwargs...) -> LikelihoodFitResult

Fit an inhomogeneous Poisson point process. rate(x, p) is the event intensity, not a normalized density. domain=(a, b) defines the integration range for the expected total event count. The domain endpoints must be finite, and all observations must lie inside the domain. With vectorized=true, rate(xs, p) returns one intensity per input point: all observations are evaluated together, and QuadGK batches the quadrature nodes.

The objective is 2 * integral(rate, domain) - 2 * sum(log(rate(x_i, p))). The parameter-independent term 2*log(n!) of the extended likelihood, with n the number of observations, is omitted: fits, covariance, and AIC differences on the same data are unaffected, but absolute -2 log(L) values differ from fully normalized implementations by that constant. No generic chi-square p-value is reported. Positive rtol controls the estimated Gauss-Kronrod error in the maximum norm, including all carried AD coefficients. Common parameter-control and solver keywords are listed under fit_custom. Invalid domains, non-positive rates, and integration/model failures raise an error. Returns LikelihoodFitResult.

fit_indexed_model(model, indices, y; p0, sigma_y=nothing, cov_y=nothing,
                  kwargs...) -> LikelihoodFitResult

Fit observations addressed by arbitrary indices. model(indices, p) must return one model value per index. This is useful when the independent variable is not a numeric 1D x-axis. Optional sigma_y entries must be finite and positive; optional cov_y must be a finite symmetric positive-definite covariance matrix.

The minimized objective is chi-square without additive Gaussian normalization constants. AIC/BIC may compare models fit to the same observations and the same uncertainty model, but not different uncertainty scales or datasets.

Use either sigma_y or cov_y, never both. With neither, the objective is the raw sum of squared residuals: chi2, chi2_ndf, and pvalue then treat every observation as having variance one in the units of y, and param_stderr is likewise computed for unit variance. Report these statistics only when the data truly have unit variance; otherwise supply sigma_y or cov_y. Common parameter-control and solver keywords are listed under fit_custom. Mismatched lengths of indices, y, sigma_y, cov_y, or the model output raise DimensionMismatch; non-positive sigma_y entries raise DomainError; other invalid uncertainty or covariance input raises ArgumentError. Returns LikelihoodFitResult.

fit_multi_model(models, xs, ys; p0, sigma_y=nothing,
                parameter_map=nothing, kwargs...) -> LikelihoodFitResult

Fit multiple datasets simultaneously with one shared global parameter vector. By default each models[i](xs[i], p) receives the full parameter vector p. With parameter_map, model i receives p[parameter_map[i]] instead. Per-dataset sigma_y entries are treated as physical standard deviations and must be finite and positive.

The minimized objective is the summed chi-square without additive Gaussian normalization constants. AIC/BIC may compare models fit to the same datasets and uncertainty model, but not different uncertainty scales or datasets.

sigma_y is either nothing or one uncertainty vector (or nothing) per dataset. A nothing entry leaves that dataset unweighted: its raw squared residuals enter the objective with unit weight relative to the weighted datasets, and chi2, chi2_ndf, and pvalue then treat those observations as having variance one. parameter_map[i] contains one-based indices into the global p0 and defines the local parameter order passed to model i. Common parameter-control and solver keywords are listed under fit_custom. Empty dataset collections, out-of-range parameter_map indices, and non-finite entries raise ArgumentError; mismatched collection, per-dataset x/y, sigma_y, or model-output lengths raise DimensionMismatch; non-positive uncertainties raise DomainError. Returns LikelihoodFitResult.

LikelihoodFitProblem(objective, gof, p0; nobs, cost_name, kwargs...)

Low-level public problem representation for likelihood and custom-objective fits. objective(p) is minimized directly, while optional gof(p) supplies a chi-square-like data goodness-of-fit statistic for reduced statistics and p-values. ScientificFitting adds the chi-square contributions from Gaussian parameter terms (the parameter_priors and parameter_constraints entries) to that statistic, matching their treatment as auxiliary observations in the degrees of freedom. The stored objective is the original callback: it does not add the Gaussian parameter terms or enforce bounds, fixed values or nonlinear constraints. For a normalized data -2log(L) callback, -objective(p)/2 is therefore the data log likelihood; external posterior models must supply priors and parameter support themselves. An arbitrary custom objective need not have this likelihood interpretation. Nonlinear constraint callbacks receive the same complete p vector, including fixed parameters; ScientificFitting handles the reduced optimizer coordinates internally. Starting values must be finite. Fixed values are validated against declared bounds while the problem is constructed. For covariance, profile thresholds, and information criteria to have their documented interpretation, objective must use the -2 log(L) scale. nobs counts the statistically independent observations entering objective; it must be positive and sets the degrees of freedom and the BIC sample size. Gaussian parameter terms are added to this count automatically, so do not include them in nobs. cost_name is a Symbol stored as the cost label in FitStatistics and echoed in fit reports; it does not affect the minimization. Most users should prefer fit_poisson_model, fit_histogram_model, fit_unbinned_model, fit_extended_unbinned_model, or fit_custom; construct this type directly when the objective must be stored, inspected, or refitted. derivatives=:finite selects numerical derivatives for foreign/Float64-only objectives and constraint callbacks, including the covariance Hessian and all profile refits. The default :auto uses ForwardDiff. Both modes require a smooth objective near the evaluation point. When the problem is later fitted, omitting tol there selects default_fit_tolerance for the solver and derivative mode; explicit tolerances are preserved.

fitted_model(result::LikelihoodFitResult)

Reconstruct the distribution or event-intensity model of a fit_distribution result using its fitted full parameter vector. The returned object has the upstream package's native interface. For BuildConstructors fits this uses the private constructor snapshot, even if the caller later updates the original. No refit is performed. Arbitrary fit_custom results have no stored model factory and raise ArgumentError rather than guessing how to reconstruct one.

AbstractFitSolver

Optional scalar-minimization adapter. Implement solver_capabilities and solve_fit; optionally specialize default_fit_tolerance. The statistical model stays in ScientificFitting.

OptimizationSolver(algorithm; kwargs...)

Use an Optimization.jl algorithm with its native solve keywords. Load the corresponding Optimization solver package first. Set maxiters and tol on the fit, not here. Other keywords (for example g_tol or store_trace) are preserved in profile refits. Unsupported bounds/constraints are rejected using SciML's algorithm traits. OptimizationSolver minimizes the cost as one scalar; it does not exploit least-squares structure (residual Jacobians), so plain chi-squared fits are usually faster with the automatic LsqFit path. For NLopt, pass an algorithm enum (for example NLopt.LN_NELDERMEAD), not a dimension-bound NLopt.Opt: profile refits change the number of free parameters.

NativeMinuitSolver(; steps=nothing, kwargs...)

Select MIGRAD after import NativeMinuit (Julia 1.11+); importing only the module avoids a name conflict with its exported profile. steps is a positive scalar or a vector in the original complete parameter order; it specifies initial numerical step sizes, not statistical errors or priors. Native Minuit constructor options such as strategy=2 and check_gradient=false are passed through. Bounds, fixed parameters, derivatives, and errordef=1 belong to the fit and cannot be overridden here, including through native fix_*, limit_* or error_* aliases. Set numerical initial steps with steps instead.

maxiters is the native function-call budget, not an iteration count; tol is MIGRAD's tolerance on the EDM (estimated distance to minimum), defaulting to the native 0.1 (target EDM 0.002 * tol, accepted up to ten times that goal, on our errordef=1 cost scale). An explicit tol is never rescaled or relaxed. A solve that fails, or that ScientificFitting's convergence validation rejects, may restart once from its last point, using only the remaining function-call budget. No solver, strategy or tolerance is silently changed. The native object from the retained attempt is available in result.solver_result.raw; its coordinates follow result.solver_result.parameter_indices. Symmetric errors in the ScientificFitting result still follow its covariance policy.

FitSolverResult(params; backend, converged, message, iterations=missing,
                parameter_indices=eachindex(params), raw=nothing)

Solver output in free coordinates, separate from statistical inference. parameter_indices maps these coordinates to the full scientific parameter vector. backend names the adapter family (for example :optimization); converged reports the solver's own termination status, which ScientificFitting validates independently; iterations stays missing when the solver reports none; message is the human-readable termination reason. raw preserves native status/covariance/trace information; mutating it does not update the already constructed ScientificFitting result.

solver_capabilities(solver) -> NamedTuple

Declare bounds, constraints, gradient, and hessian Boolean flags. The last two indicate required derivatives. Missing adapter implementations raise MethodError, rather than silently assuming capabilities.

default_fit_tolerance(solver, derivatives::Symbol) -> Real

Stopping tolerance used when a fit omits tol or passes tol=nothing. derivatives is :auto or :finite; solver=nothing keeps the automatic choice (LsqFit for static chi-square least squares, an Optimization.jl algorithm otherwise; see fit), which uses the generic fallback below. LsqFit/Optimization use 1e-10, or 1e-6 for noisier finite differences. NativeMinuit uses its native EDM (estimated distance to minimum) tolerance 0.1 in either derivative mode.

Third-party solvers may specialize this method to match their stopping rule. Return a finite positive value. The fit stores the resolved tolerance in result.options.tol and reuses it in multistart and profile refits. Explicit tolerances bypass this method. These are numerical criteria, not statistical error guarantees; equal values need not mean equal accuracy across solvers.

solve_fit(solver, problem::OptimizationProblem; maxiters, tol,
          parameter_indices, parameter_count, parameter_names) -> FitSolverResult

Public solver-extension boundary. problem.f(q, problem.p) is the complete validated cost; q contains only free parameters. Bounds and nonlinear constraints are already mapped to that order. Do not add priors, rescale the objective, modify inputs, or compute ScientificFitting's statistical summaries. parameter_count is the original full dimension (needed for parameter-specific solver settings). Retain unavailable iteration counts as missing and failed termination as converged=false. The core checks capabilities before dispatch. parameter_names are unique labels in free-coordinate order, suitable for named native solver operations. Unnamed problems use p1, p2, etc. tol is already a finite positive number, resolved by default_fit_tolerance unless the user supplied it explicitly.

ConstraintSpec(; ineq=nothing, eq=nothing)

Container for nonlinear parameter constraints passed to the general Optimization.jl backend. ineq and eq receive the complete parameter vector in the order defined by p0, including fixed parameters. Each callback returns a scalar or a vector; every component is constrained elementwise (ineq(p) .<= 0, eq(p) .== 0). Use this only when simple bounds, fixed parameters, or Gaussian parameter constraints are not expressive enough.

ParameterPrior(index, mean, sigma)
ParameterPrior(index, mean, sigma_minus, sigma_plus)

Gaussian prior term for one fitted parameter. The prior adds ((p[index] - mean) / sigma)^2 to the objective and counts as one observation in ndf; asymmetric uncertainties use sigma_minus below the mean and sigma_plus above the mean. With cost=:gaussian_likelihood, asymmetric scales define a continuous split-normal density with one shared normalization across both sides of mean.

FixedParameter(index, value[, sigma])
FixedParameter(index, value, sigma_minus, sigma_plus)

Fix one parameter to value during the fit. Optional uncertainties describe the externally known value. ScientificFitting stores them in reports and in that fixed parameter's diagonal covariance entry with zero fitted cross-covariances. With asymmetric uncertainties the diagonal entry is max(sigma_minus, sigma_plus)^2; a single covariance entry cannot represent asymmetry, so the larger side is used. They do not change the objective, make the parameter free, or propagate uncertainty into the fitted parameters; use ParameterPrior or ParameterConstraint for that coupled statistical treatment.

ParameterConstraint(indices, mean, covariance)

Correlated Gaussian constraint on several parameters. With delta = p[indices] .- mean, the constraint adds the penalty delta' * inv(covariance) * delta (evaluated via a Cholesky factorization, never an explicit inverse) to the chi-square objective: the selected parameters are treated as one correlated auxiliary measurement with the given mean and positive-definite covariance. With cost=:gaussian_likelihood, the parameter-independent normalization length(indices) * log(2pi) + logdet(covariance) is added as well. Each constrained component counts as one auxiliary observation in ndf; the External Parameter Information section of the statistics reference derives the scalar analogue and the ndf counting.

ErrorComponent(name, target, mode, values; active=true)

Named uncertainty contribution used by component-based covariance models. target is :x or :y. For mode=:absolute, values are scalar or pointwise standard deviations. :relative interprets them as fractions of the absolute measured values, i.e. sigma_i = values_i * |y_i| (or |x_i| for target=:x); the y-only :model_relative uses |model(x_i, p)| instead. For :covariance, a vector contains pointwise standard deviations and a matrix is the complete covariance contribution. Active components add in covariance space; active=false keeps a source documented but excludes it from the current fit.

WhiteningOperator(whiten!; logdet_covariance, marginal_sigma=nothing)

Represent a complete static data covariance without materializing it. whiten!(out, residual) must apply a linear operator W satisfying W'W = inv(C), where C is the observation covariance. The required logdet_covariance is log(det(C)); ScientificFitting uses it for the normalized Gaussian -2 log(L) cost, AIC, and BIC.

marginal_sigma is optional scalar or pointwise marginal standard deviation. It does not change the fit; plotting uses it for data error bars and pointwise prediction bands. Without it, confidence bands remain available, but requesting a prediction band (band=:prediction) raises an ArgumentError.

The mutating function must write every element of out and support the element types used by automatic differentiation when the general optimizer is needed, unless the fit uses derivatives=:finite for Float64-only callbacks. It must accept AbstractVector views because ScientificFitting applies the same operator columnwise to analytic Jacobians.