clophfit.fitting.noise_calibration#

Noise-model calibration from fit residuals.

Estimators that turn a canonical residual table into a PlateNoiseModel: per-label floor, photon gain, and proportional error, plus the plate slope helpers used to propagate x-axis noise.

Functions#

`calibrate_noise_robust`(residuals, sigma_floor, *[, ...])	Calibrate a per-label noise model from outlier-screened residuals.
`fit_rel_error_from_residuals`(df, sigma_floor)	Estimate proportional error (alpha) per label via moment estimator.
`fit_gain_from_residuals`(df, sigma_floor)	Estimate Poisson gain per label via moment estimator.
`fit_noise_model_nnls`(df[, sigma_floor_fixed, ...])	Fit heteroscedastic noise model via non-negative least squares.
`fit_noise_model_merged`(df, *[, n_bins])	Fit one gain and one alpha per label, on every plate's residuals at once.
`dof_scale`(n_points, n_params)	Return `sqrt(n / (n - p))`, the factor undoing what a fit absorbed.
`fit_gain_nnls`(df, *[, fit_floor])	Fit the gain-only noise model to residuals pooled over wells and plates.
`fit_gain_floor_fixed`(df, sigma_floor)	Estimate one plate's gain per label, floor held and alpha at 0.
`fit_gain_and_floor`(df)	Estimate one plate's gain and floor per label, alpha at 0.
`compute_binding_slope`(ph, pka, s0, s1)	Compute \|dS/dpH\| for the Henderson-Hasselbalch equation.
`compute_plate_slopes`(results)	Compute per-well per-label `∂S/∂pH` from pass-1 fit results.
`fit_ph_slope_noise`(df, noise_model, plate_slopes)	Fit global `sigma_ph` from excess variance after per-label model.

Module Contents#

clophfit.fitting.noise_calibration.calibrate_noise_robust(residuals, sigma_floor, *, p_threshold=0.9, min_keep=3)#

Calibrate a per-label noise model from outlier-screened residuals.

Drops points whose posterior outlier probability exceeds p_threshold (from a PyMC mixture fit) and then estimates gain and alpha per label with the single-term moment estimators (fit_gain_from_residuals(), fit_rel_error_from_residuals()) on the retained points. Screening with the mixture’s p_outlier keeps outliers from inflating the estimate, while the two single-term estimators avoid the gain/alpha collinearity of the joint NNLS over narrow titration ranges.

Parameters:

residuals (pd.DataFrame) – Canonical residual table (e.g. MultiFitResult.residuals from a mixture fit). Must have label, raw_res, yhat columns; a p_outlier column enables screening (otherwise all points are used), and a dof_scale column corrects the residuals for the parameters the fit spent. Supply it: a mixture fit has no residual degrees of freedom this function could infer.
sigma_floor (dict[str, float]) – Known read-noise floor per label, e.g. tit.bg_noise. Used as the fixed floor and copied into the returned model.
p_threshold (float) – Posterior outlier probability above which a point is dropped.
min_keep (int) – Per label, if screening would retain fewer than this many points the full (unscreened) set is used instead.

Returns:

Per-label model with sigma_floor from sigma_floor and calibrated gain/alpha.

Return type:

PlateNoiseModel

clophfit.fitting.noise_calibration.fit_rel_error_from_residuals(df, sigma_floor)#

Estimate proportional error (alpha) per label via moment estimator.

Assumes the simplified noise model sigma^2 = floor^2 + alpha^2 * yhat^2 (no Poisson gain term). With floor known from buffer measurements and using model-predicted values yhat in the denominator to avoid noise-in-variables bias, the closed-form moment estimator is:

\[\hat{\alpha}^2 = \frac{\overline{r^2} - \sigma_{\text{floor}}^2}{\overline{\hat{y}^2}}\]

Parameters:

df (pd.DataFrame) – DataFrame with columns label (str), raw_res (float), and yhat (float – the model-predicted signal at each point), and optionally dof_scale (see _squared_residuals()). Typically from clophfit.fitting.model_validation.residuals_from_fit_results().
sigma_floor (dict[str, float]) – Known read-noise floor per label, e.g. from tit.bg_noise.

Returns:

Per-label proportional error estimate alpha (non-negative).

Return type:

dict[str, float]

Examples

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> y_pred = np.linspace(50, 500, 200)
>>> floor, true_alpha = 5.0, 0.02
>>> sigma = np.sqrt(floor**2 + (true_alpha * y_pred) ** 2)
>>> resid = sigma * rng.standard_normal(200)
>>> df = pd.DataFrame({"label": "1", "raw_res": resid, "yhat": y_pred})
>>> alpha = fit_rel_error_from_residuals(df, sigma_floor={"1": floor})
>>> round(alpha["1"], 2)  # should be close to true_alpha=0.02
0.02

clophfit.fitting.noise_calibration.fit_gain_from_residuals(df, sigma_floor)#

Estimate Poisson gain per label via moment estimator.

Symmetric counterpart to fit_rel_error_from_residuals(). Assumes the Poisson-only noise model sigma^2 = floor^2 + gain * yhat (no proportional term), which sidesteps the gain/alpha collinearity of the joint fit. With floor known from buffer measurements and using model-predicted values yhat, the closed-form moment estimator is:

\[\hat{\text{gain}} = \frac{\overline{r^2} - \sigma_{\text{floor}}^2}{\overline{\hat{y}}}\]

Parameters:

df (pd.DataFrame) – DataFrame with columns label (str), raw_res (float), and yhat (float – the model-predicted signal at each point), and optionally dof_scale (see _squared_residuals()). Typically from clophfit.fitting.model_validation.residuals_from_fit_results().
sigma_floor (dict[str, float]) – Known read-noise floor per label, e.g. from tit.bg_noise.

Returns:

Per-label Poisson gain estimate (non-negative).

Return type:

dict[str, float]

Examples

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> y_pred = np.linspace(50, 500, 400)
>>> floor, true_gain = 5.0, 0.8
>>> sigma = np.sqrt(floor**2 + true_gain * y_pred)
>>> resid = sigma * rng.standard_normal(400)
>>> df = pd.DataFrame({"label": "1", "raw_res": resid, "yhat": y_pred})
>>> gain = fit_gain_from_residuals(df, sigma_floor={"1": floor})
>>> round(gain["1"], 1)  # should be close to true_gain=0.8
0.8

clophfit.fitting.noise_calibration.fit_noise_model_nnls(df, sigma_floor_fixed=None, rel_error_fixed=None, gain_fixed=None)#

Fit heteroscedastic noise model via non-negative least squares.

Model: \(\sigma^2 = \sigma_\text{floor}^2 + \text{gain} \cdot y + \alpha^2 \cdot y^2\)

Uses scipy.optimize.nnls() to enforce non-negativity on all parameters, which stabilises estimates when \(y\) and \(y^2\) are highly collinear (typical for narrow-range titrations).

Parameters:

df (pd.DataFrame) – Residual DataFrame with columns label, raw_res, yhat, and optionally dof_scale, sqrt(n / (n - p)) per row, which undoes the variance the fit absorbed (see _squared_residuals()). Every caller holding a fit should pass it; without it the terms come out low by that factor squared.
sigma_floor_fixed (dict[str, float] | None) – If given, fix floor per label and only fit gain and alpha.
rel_error_fixed (dict[str, float] | None) – If given, fix alpha per label and only fit the remaining terms.
gain_fixed (dict[str, float] | None) – If given, fix gain per label and only fit the remaining terms. Any subset may be fixed, including all three, which simply returns them.

Returns:

(sigma_floor, gain, alpha) per label — all non-negative.

Return type:

tuple[dict[str, float], dict[str, float], dict[str, float]]

clophfit.fitting.noise_calibration.fit_noise_model_merged(df, *, n_bins=10)#

Fit one gain and one alpha per label, on every plate’s residuals at once.

The pooled estimator in its literal sense: the residual tables from all plates are concatenated and a single gain and alpha are fitted to the lot, rather than fitting each plate and averaging the answers. Fitting per plate and summarising cannot recover a shared value when the per-plate fits are themselves at a boundary, which they routinely are - y and y**2 are near-collinear over the range one plate covers.

The floor is supplied, never fitted, and differs per plate when it is scaled to that plate’s reader Gain. Each plate’s own floor variance is therefore removed before pooling, and what is fitted is the remainder:

\[\sigma^2 - \text{floor}^2 = \text{gain} \cdot a \cdot y + (\alpha \cdot y)^2\]

where \(a\) is gain_amp. Binning happens within a plate, so a bin never mixes plates of different brightness, and bins are weighted by count.

Parameters:

df (pd.DataFrame) – Residual table with columns plate, label, yhat, raw_res and sigma_floor. Two columns are optional: gain_amp, the amplification the gain term is quoted against (default 1.0), and dof_scale (default 1.0). Pass gain_amp only for a label whose noise actually tracks reader Gain - scaling a label whose floor is supplied unscaled quotes the two terms at different references.
n_bins (int) – Equal-count signal bins per plate and label.

Returns:

One row per label: label, n_bins, n_plates, gain and alpha (both at gain_amp == 1), log_misfit and negative_bins.

Return type:

pd.DataFrame

Raises:

ValueError – If a required column is missing. The floor in particular is an input, because a fit that frees all three terms drives one to its boundary.

clophfit.fitting.noise_calibration.dof_scale(n_points, n_params)#

Return sqrt(n / (n - p)), the factor undoing what a fit absorbed.

A least-squares fit pulls its own residuals towards the data, so their mean square is (n - p) / n of the noise variance. For a two-label pH well (16 points, K plus two plateaus per label) that is 0.69, and a gain read from the uncorrected residuals comes out a third low; on simulated wells this factor restores it to within 2%, as closely as the full per-point leverage.

Parameters:

n_points (int) – Observations the fit used.
n_params (int) – Parameters it varied.

Returns:

The factor, or nan when the fit has no residual degrees of freedom.

Return type:

float

clophfit.fitting.noise_calibration.fit_gain_nnls(df, *, fit_floor=False)#

Fit the gain-only noise model to residuals pooled over wells and plates.

\[\sigma^2 = \text{floor}_p^2 + \text{gain} \cdot a_p \cdot \hat{y}\]

with alpha held at 0, which removes the gain/alpha ridge that y and y**2 form over a titration’s narrow range. One gain per label is shared by every plate in df; the floor is each plate’s own, either as supplied (fit_floor=False) or fitted alongside the gain (fit_floor=True). With a single plate this is the per-plate calibration, with many it is the pooled one, and nothing else differs.

Squared residuals are first multiplied by dof_scale**2, since a fit’s residuals under-state the noise by the degrees of freedom it used (see dof_scale()). The fit is non-negative least squares on the squared residuals, reweighted by 1 / sigma^2 - the inverse standard deviation of a squared Gaussian residual - until the terms stop moving: the Gaussian maximum-likelihood fit of the variance (see _variance_ml()).

Parameters:

df (pd.DataFrame) – Residual table with label, yhat and raw_res. Optional: plate (default one plate, ""), dof_scale (default 1), gain_amp (default 1; the amplification the gain is quoted against, for a label whose shot noise tracks reader Gain) and sigma_floor, required unless fit_floor.
fit_floor (bool) – Fit each plate’s floor with the gain instead of holding it. Over a range where the gain term dominates, the floor is poorly determined and often lands at 0.

Returns:

One row per label and plate: label, plate, gain (at gain_amp == 1, shared by the label’s plates), sigma_floor, floor_fitted and n.

Return type:

pd.DataFrame

Raises:

ValueError – If a required column is missing, or the floor is to be held but a row has none.

Examples

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> y = np.linspace(50, 500, 400)
>>> r = np.sqrt(3.0**2 + 0.8 * y) * rng.standard_normal(400)
>>> df = pd.DataFrame({"label": "1", "yhat": y, "raw_res": r, "sigma_floor": 3.0})
>>> round(float(fit_gain_nnls(df)["gain"].iloc[0]), 1)
0.8

clophfit.fitting.noise_calibration.fit_gain_floor_fixed(df, sigma_floor)#

Estimate one plate’s gain per label, floor held and alpha at 0.

Parameters:

df (pd.DataFrame) – One plate’s residual table (label, yhat, raw_res, and optionally dof_scale).
sigma_floor (Mapping[str, float]) – The measured read-noise floor per label.

Returns:

Gain per label.

Return type:

dict[str, float]