Python Interface

NumPy model functions, the Julia numerical core, editable Matplotlib figures. Python 3.10+ is required; Makie is neither installed nor loaded.

# [sparse]: SciPy distributions and sparse covariance matrices.
python -m pip install 'scientificfitting[plot,sparse]'

Omit the extras for fitting and reports only. JuliaCall provisions Julia and the compatible 0.3.x core on first use, with network access and compilation; the Julia runtime's size and startup cost remain. Importing scientificfitting alone does not start Julia.

PyPI. A Conda-forge recipe is under review; until it lands, the pip command also works inside a Conda environment.

Fit, Inspect, Plot

Code cells run in order in one Python session.

import numpy as np
import matplotlib.pyplot as plt
from scientificfitting import fit_model, plot_fit

def line(x, slope, offset):
    return slope * x + offset

x = np.array([0., 1., 2., 3.])
y = np.array([0.1, 1.2, 1.9, 3.2])
result = fit_model(line, x, y, p0={"slope": 1., "offset": 0.}, sigma_y=0.2)
print(result.report())
print(result.diagnose())

# A normal Matplotlib axis: adding a marker does not repeat the fit.
fig, ax = plot_fit(result, xlabel="x / mm", ylabel="U / V")
ax.axvline(1.5, color="black", linestyle="--")
fig.savefig("calibration.pdf")
plt.close(fig)

Callback contract: model(x, **parameters) returns one prediction per coordinate. x is read-only. p0 keys bind parameters by name and set result-array order. Use deterministic models; derivatives use finite differences, including profile refits. Optional jacobian and x_derivative callbacks differentiate the unweighted model.

# Read values directly; never parse the terminal report.
slope_value = result.values["slope"]
report = result.report(structured=True)
slope = report.parameters["slope"]
print(slope.value, slope.uncertainty, slope.fixed)

diagnosis = result.diagnose(structured=True, max_actions=3)
for finding in diagnosis.findings:
    print(finding.code, finding.recommendation)

# Local uncertainty of the fitted mean, excluding observation noise.
prediction, uncertainty = result.predict(x, uncertainty=True)

Arrays are independent, read-only snapshots. result.converged records solver convergence, not model validity; an unknown iteration count is None. diagnosis.status == "ok" means the implemented checks found no issue.

Choose The Error Model

The entry point is selected by the observation model; the table in Choose An Entry Point applies with the same function names, with two exceptions. fit_distribution is Julia-only: in Python, write the density explicitly and use fit_unbinned_model (unbinned observations) or fit_histogram_density (binned counts). The model-contract column shows the Julia convention model(x, p); Python callbacks receive parameters by keyword, model(x, **parameters), as above. Gaussian errors enter as sigma_x/sigma_y, dense or SciPy sparse covariance as cov_x/cov_y, and named error sources as ErrorComponent inputs (help(ErrorComponent)).

Student-t errors allow heavier tails; here the scales are known separately for each observation:

from scipy import stats
from scientificfitting import fit_likelihood_model

locations = np.array([-1., -0.4, 0., 0.5, 1.2, 1.8, 2.4])
readings = np.array([-1.1, -0.23, 0.31, 0.94, 2.04, 2.91, 6.9])
scales = np.array([0.12, 0.20, 0.14, 0.18, 0.11, 0.25, 0.20])

def measurement_logprob(y, prediction, slope, offset):
    # One normalized log density per independent observation, not a sum.
    return stats.t.logpdf(y, df=4, loc=prediction, scale=scales)

robust_result = fit_likelihood_model(
    line, locations, readings, logprob=measurement_logprob,
    p0={"slope": 1., "offset": 0.},
)
print(robust_result.report())

The Student-t standard deviation is $\sqrt{df/(df-2)}\times$ scale ($\sqrt{2}$ for df=4, defined for df > 2). Choose the distribution from the error process, not to conceal a wrong mean model. Discrete observations need a log mass (logpmf); no universal goodness-of-fit p-value is assigned to this custom distribution.

from scientificfitting import add_report, plot_style

with plt.rc_context(plot_style("sans")):
    fig, ax = plt.subplots(layout="constrained")
    grid = np.linspace(locations.min(), locations.max(), 300)
    ax.plot(grid, line(grid, **robust_result.values), color="#0072B2", label="fitted mean")
    # Error-distribution quantiles, NOT uncertainty of the fitted line.
    ax.errorbar(locations, readings, yerr=stats.t.ppf(0.84, df=4)*scales,
                fmt="o", color="black", markersize=3, elinewidth=0.8, capsize=2,
                label="data; 16-84% error range")
    ax.set(xlabel="reference setting", ylabel="response / V", title="Student-t errors")
    add_report(fig, robust_result, ax=ax, statistics=("cost_min",),
               statistic_labels={"cost_min": r"$-2\log L$"}, expand=True)
    fig.savefig("student_t_errors.pdf")
    plt.close(fig)

Non-Smooth Likelihoods

A Laplace location fit estimates the sample median; its cost has corners, so use a derivative-free solver:

readings_laplace = np.array([-1.2, -0.1, 0.2, 0.4, 0.8, 1.3, 5.0])
def constant(x, location):
    return np.full_like(x, location)

laplace_result = fit_likelihood_model(
    constant, np.arange(len(readings_laplace)), readings_laplace,
    logprob=lambda y, mu, location: stats.laplace.logpdf(y, loc=mu, scale=1),
    p0={"location": 0.1}, solver="nelder_mead", tol=1e-10,
)
print(laplace_result.report())
# Nelder-Mead reports no local errors, so the automatic profile grid
# is unavailable: pass explicit scan values.
scan = laplace_result.profile("location", values=[-0.2, 0., 0.4, 0.7, 1.])

Nelder-Mead supports bounds, fixed values and Gaussian parameter terms, but not nonlinear constraints. It defaults to parameter_covariance="none" (free-parameter errors are reported as NaN, not zero); maxiters limits function evaluations. Select "hessian" only for a locally smooth cost. Zero probability returns -np.inf; start at a finite likelihood. Moving support can invalidate standard profile thresholds: see the support-boundary calculation.

Profiles And Contours

Four exposures with one observed event give an asymmetric rate interval:

from scientificfitting import fit_poisson_model

counts = np.array([0, 0, 1, 0])
count_fit = fit_poisson_model(
    lambda x, rate: np.full_like(x, rate), np.arange(4), counts,
    p0={"rate": 0.3}, bounds={"rate": (0.001, 10)},
)
interval = count_fit.profile_interval("rate", npoints=61, nsigma=4)
print(interval.lower, interval.upper)  # expected counts per exposure

The default threshold=1 cut on the profile cost has an asymptotic 68.27% interpretation; at low counts it is not exact. Two-parameter 68.27%/95.45% regions use [2.30, 6.18] (Profiles and Contours). Missing crossings remain NaN; failed refits remain gaps.

from scientificfitting import (plot_contour, plot_diagnostics, plot_profile,
                              plot_profile_matrix)

# Explicit selection and grid sizes bound the number of nuisance refits.
matrix = result.profile_matrix(
    ["slope", "offset"], npoints_profile=9, npoints_contour=7, nsigma=3,
)
pair = matrix.contours["slope", "offset"]
review = {k: v for k, v in matrix.panel_status.items() if v != "ok"}

# Rendering completed scans does not minimize again.
fig, ax = plot_profile(interval.profile_result, delta_max=5)
plt.close(fig)
fig, ax = plot_contour(pair)
plt.close(fig)
fig, axes = plot_profile_matrix(matrix)
fig.savefig("calibration_profiles.pdf")
plt.close(fig)

# Reassess stored scans without changing their confidence thresholds.
profile_review = interval.profile_result.diagnose(structured=True)
pair_review = pair.diagnose(tolerance=0.5, structured=True)

tolerance above measures departure from the local parabola/ellipse, not confidence level. Arrays store delta_cost[i, j] at (x[i], y[j]); direct Matplotlib calls need pair.delta_cost.T. result.report(errors="profile") computes asymmetric errors with additional fits; its covariance remains local.

Customize With Matplotlib

Style ("sans" / "tex") and report visibility (panel=True / False) are independent. "tex" uses bundled STIX fonts and MathText, not external TeX. rc_context leaves global settings unchanged.

with plt.rc_context(plot_style("tex")):
    fig, ax = plot_fit(result, panel=False, xlabel="x / mm", ylabel="U / V")
    # Add labeled artists before constructing the report legend.
    ax.axvline(1.5, color="black", linestyle=":", label="reference position")
    panel = add_report(
        fig, result, ax=ax, expand=True,
        parameter_labels={"slope": r"$m$", "offset": r"$b$"},
        statistics=("chi2_ndf", "pvalue"),
    )
    fig.savefig("calibration_with_reference.pdf")
    plt.close(fig)

fig, axes = plot_diagnostics(result, kinds=("residual", "pull"), xlabel="x / mm")
fig.savefig("calibration_diagnostics.pdf")
plt.close(fig)

Use layout="constrained" for outside reports on your figures. expand=True allows canvas growth; an explicitly fixed tiny canvas stays too small. plot_fit and x-y residual helpers target Gaussian results; other fit families use ordinary Matplotlib plus add_report. Correlated pulls are whitened residuals (Residuals and Pulls); their unit bands are reference guides, not coverage intervals.

Large Datasets And In-Place Models

SciPy sparse cov_x/cov_y inputs are not densified. Static y-covariance reuses its factorization; parameter-dependent covariance does not. Sparse factors can still fill in: the factorization of a sparse covariance can be much denser than the matrix itself, so memory use can grow. A known whitening operator avoids storing the matrix (Structured Whitening):

from scientificfitting import WhiteningOperator

# Stationary AR(1): C[i,j] = sigma**2 * rho**abs(i-j), equal sample spacing.
sigma, rho, n = 0.2, 0.6, len(y)
def whiten(out, residual):
    out[0] = residual[0] / sigma
    # Remove the predictable neighbor contribution, then standardize.
    out[1:] = (residual[1:] - rho * residual[:-1]) / (sigma * np.sqrt(1 - rho**2))

noise = WhiteningOperator(
    whiten, 2*n*np.log(sigma) + (n-1)*np.log1p(-rho**2),  # log(det(C))
    marginal_sigma=sigma, inplace=True,  # marginal_sigma is only plot metadata
)
correlated = fit_model(line, x, y, p0={"slope": 1., "offset": 0.}, whitening=noise)

def line_inplace(out, x, slope, offset):
    np.multiply(x, slope, out=out)
    out += offset

inplace_result = fit_model(line_inplace, x, y, p0={"slope": 1., "offset": 0.},
                           sigma_y=0.2, inplace=True)

For AR(1), $\det(C) = \sigma^{2n}(1-\rho^2)^{n-1}$, which gives the log-determinant above. Verify a custom operator's log-determinant against np.linalg.slogdet of a small dense $C$ before using it at scale.

Whitening replaces other observation errors. In-place callbacks fill every entry and return None; never retain the borrowed arrays. An in-place Jacobian fills an (n, k) matrix, columns in p0 order. Without x_derivative, the model must also be defined near the measured coordinates for finite differencing.

Event Densities And Complete Workflows

For uncensored positive waiting times (Unbinned and Extended Likelihoods), evaluate the exponential density in one NumPy call:

from scientificfitting import fit_unbinned_model

waiting_times = np.array([0.12, 0.28, 0.51, 0.62, 0.75, 1.3, 1.8])
def waiting_pdf(t, tau):
    return np.exp(-t/tau)/tau

waiting_fit = fit_unbinned_model(
    waiting_pdf, waiting_times, p0={"tau": 0.5}, bounds={"tau": (0.01, 5.)},
    vectorized=True,  # one array call, not one call per event
)

tau estimates the mean waiting time. A detection threshold requires a correspondingly normalized truncated density. vectorized=True also supports extended likelihoods and adaptive histogram integration; rtol controls the integral accuracy.

The count bands drawn in likelihood_workflows.py are Poisson quantile bands: they describe count fluctuations conditional on the fitted mean, not parameter uncertainty. For backgrounds near zero, use profile intervals; local symmetric errors can cross the physical lower bound.

Development Setup

Only contributors need a local Julia checkout. From the repository root, in a separate Python environment:

python -m pip install -e './python[plot,test]'
python python/develop.py
python examples/python/numpy_matplotlib.py
python -m pytest python/tests

develop.py persistently selects this checkout in its JuliaPkg environment. Use a fresh environment for registry-installation checks. Installed-wheel checks live in python/tests/check_install.py; register the Julia core before publishing a Python release that depends on it. Python startup and restart-latency measurements live in python/benchmarks/startup.py (run it after installation and precompilation; an empty JuliaPkg environment additionally times provisioning), callback overhead in python/benchmarks/callbacks.py. The Julia core's reproducible benchmarks and fresh-process loading gate are described under Performance Checks.