Metadata-Version: 2.4
Name: mcup
Version: 1.0.0
Summary: Monte Carlo Uncertainty Propagation for regression with measurement errors
Author-email: Daniel Herman <daniel.herman@protonmail.com>
License: MIT
Project-URL: Repository, https://github.com/detrin/MCUP
Keywords: physics,statistics,error,uncertainty,propagation,regression
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.21.0
Requires-Dist: scipy>=1.8.0
Requires-Dist: numdifftools>=0.9.39
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: pytest-cov>=4.0; extra == "dev"
Requires-Dist: ruff>=0.4.0; extra == "dev"
Requires-Dist: mypy>=1.8; extra == "dev"
Requires-Dist: twine>=4.0; extra == "dev"
Requires-Dist: build>=0.10; extra == "dev"
Provides-Extra: docs
Requires-Dist: mkdocs>=1.4; extra == "docs"
Requires-Dist: mkdocs-material>=9.0; extra == "docs"
Requires-Dist: mkdocstrings[python]>=0.20; extra == "docs"
Dynamic: license-file

# MCUP

MCUP (Monte Carlo Uncertainty Propagation) is a Python library for regression with measurement errors. It provides three sklearn-like estimators that correctly propagate x and y measurement uncertainties into parameter confidence intervals.

[![PyPI pyversions](https://img.shields.io/pypi/pyversions/mcup.svg)](https://pypi.org/project/mcup/) [![PyPI version](https://img.shields.io/pypi/v/mcup.svg)](https://pypi.org/project/mcup/) [![CI](https://github.com/detrin/MCUP/actions/workflows/package-main.yml/badge.svg)](https://github.com/detrin/MCUP/actions/workflows/package-main.yml) [![codecov](https://codecov.io/gh/detrin/MCUP/branch/master/graph/badge.svg?token=Dx6elQkztR)](https://codecov.io/gh/detrin/MCUP)

## Install

```bash
uv add mcup
```

Or with pip:

```bash
pip install mcup
```

## Quick start

The core idea: you have data where measurement noise is not uniform, or x itself is measured. MCUP gives you honest parameter uncertainties in both cases.

**Case 1 — only y has errors (heteroscedastic noise)**

A photodetector where noise grows with signal: points at high intensity are less reliable. OLS doesn't know that and gives overconfident slope uncertainty. `WeightedRegressor` down-weights noisy points and produces calibrated intervals.

```python
import numpy as np
from mcup import WeightedRegressor

rng = np.random.default_rng(42)
x = np.linspace(1, 10, 30)
y_err = 0.1 * x                  # noise grows with x
y = 2.0 * x + 1.0 + rng.normal(0, y_err)

def line(x, p):
    return p[0] + p[1] * x

# Uniform weights (wrong — ignores that high-x points are noisier)
ols = WeightedRegressor(line, method="analytical")
ols.fit(x, y, y_err=np.ones_like(x), p0=[0.0, 1.0])

# Correct weights from measurement errors
wls = WeightedRegressor(line, method="analytical")
wls.fit(x, y, y_err=y_err, p0=[0.0, 1.0])

print(f"OLS:      slope = {ols.params_[1]:.3f} ± {ols.params_std_[1]:.4f}  ← overconfident")
print(f"Weighted: slope = {wls.params_[1]:.3f} ± {wls.params_std_[1]:.4f}  ← calibrated")
# true slope = 2.0
```

**Case 2 — both x and y have errors**

A spring balance where both extension (x) and force (y) are measured with error. Ignoring x-errors causes attenuation bias (slope pulled toward zero) and intervals that are far too narrow. `XYWeightedRegressor` propagates both error sources.

```python
from mcup import XYWeightedRegressor

rng = np.random.default_rng(0)
x_true = np.linspace(0.1, 2.0, 25)
x_err, y_err = 0.05 * np.ones(25), 0.15 * np.ones(25)
x_obs = x_true + rng.normal(0, x_err)
y = 8.0 * x_true + rng.normal(0, y_err)   # true spring constant k=8

# Ignoring x-errors (wrong)
bad = WeightedRegressor(line, method="analytical")
bad.fit(x_obs, y, y_err=y_err, p0=[0.0, 1.0])

# Propagating both errors (correct)
est = XYWeightedRegressor(line, method="analytical")
est.fit(x_obs, y, x_err=x_err, y_err=y_err, p0=[0.0, 1.0])

print(f"Ignoring x-err: k = {bad.params_[1]:.3f} ± {bad.params_std_[1]:.3f}  ← biased low, too narrow")
print(f"XYWeighted:     k = {est.params_[1]:.3f} ± {est.params_std_[1]:.3f}  ← unbiased, calibrated")
# true k = 8.0
```

## Why MCUP

Standard least squares (OLS) assumes all observations are equally reliable. Real experiments break this in two common ways:

- **Heteroscedastic y-errors** — measurement noise varies across the range. OLS overweights noisy points, biasing the fit and producing overconfident uncertainties.
- **Errors in x** — when the independent variable is itself measured (time, concentration, displacement), ignoring those errors causes attenuation bias: slopes are pulled toward zero, and uncertainty intervals shrink below their true size.

**Why not just use the covariance matrix from the optimizer?**

When measurement errors are large, the standard approach of reading off `sqrt(diag(cov))` from the fit residuals underestimates the true parameter uncertainty. The covariance matrix tells you how well the optimizer converged — it does not propagate the uncertainty that came *in* with your data. MCUP propagates measurement noise directly through the model so that `params_std_` reflects both fit quality and input uncertainty. For a worked example comparing both approaches, see this [Kaggle notebook on measurement error in regression](https://www.kaggle.com/code/jetakow/measurement-error-in-regression).

MCUP fixes both problems. The figure below illustrates the effect for a linear calibration with heteroscedastic y-errors:

![Comparison of OLS vs WeightedRegressor](docs/assets/comparison_linear.png)

*Left: OLS (red) fits the same data differently from weighted regression (blue) because it treats all points equally regardless of σ_y. Right: over 500 simulated experiments, OLS coverage deviates from the nominal 68.3% — WeightedRegressor stays calibrated.*

## Estimators

| Estimator | Use when | Error model |
|-----------|----------|-------------|
| `WeightedRegressor` | Only y has measurement errors | `Σ (y − f(x))² / σ_y²` |
| `XYWeightedRegressor` | Both x and y have errors (nonlinear) | Combined variance via error propagation (IRLS) |
| `DemingRegressor` | Both x and y have errors (linear only) | Joint optimisation over parameters + latent true x |

Each estimator supports:
- `method="analytical"` — weighted LS + `(J^T W J)^{-1}` covariance (fast, exact for well-posed problems)
- `method="mc"` — Monte Carlo with Welford covariance (robust cross-check for nonlinear models)

## Benchmark summary

Validated across 13 physical scenarios (200 independent parameter configurations each). The analytical solver achieves well-calibrated 1σ uncertainty intervals on all scenarios:

| Scenario | Estimator | Bias | RMSE | Coverage |
|----------|-----------|------|------|----------|
| Linear calibration (homo) | WeightedRegressor | +0.3% | 12.8% | ✓ 68% |
| Linear calibration (hetero) | WeightedRegressor | +0.5% | 7.2% | ✓ 71% |
| Radioactive decay | WeightedRegressor | −0.0% | 2.6% | ✓ 64% |
| Power law (diffusion) | WeightedRegressor | +0.0% | 4.6% | ✓ 68% |
| Gaussian spectral peak | WeightedRegressor | −0.1% | 1.7% | ✓ 66% |
| Damped oscillator | WeightedRegressor | −0.4% | 7.2% | ✓ 67% |
| Exp decay + timing errors | **XYWeightedRegressor** | −1.2% | 5.0% | ✓ 64% |
| Hooke's law (x+y errors) | **XYWeightedRegressor** | −1.0% | 54% | ✓ 75% |
| Beer-Lambert (x+y errors) | **XYWeightedRegressor** | +46% | 220% | ✓ 68% |
| Method comparison | **DemingRegressor** | +14% | 111% | ✓ 64% |
| Isotope ratio MS | **DemingRegressor** | +3.2% | 420% | ✓ 72% |
| Small sample (n=8) | WeightedRegressor | −2.7% | 29% | ✓ 69% |
| Low SNR | WeightedRegressor | −1.9% | 136% | ✓ 67% |

Bias and RMSE are relative to the true parameter values. Large RMSE on near-zero intercepts (Beer-Lambert baseline, isotope intercept) reflects small absolute values — the coverage column is the reliable calibration metric.

### OLS baseline: when ignoring measurement errors breaks uncertainty estimation

Plain OLS (no error weighting) estimates parameter uncertainty from fit residuals alone — `σ² = SSR/(n−p)`. This works when noise is truly uniform. When noise varies across the range, OLS produces miscalibrated intervals even though the parameter estimates themselves may look reasonable.

| Scenario | OLS coverage | WeightedRegressor coverage | What goes wrong |
|----------|:------------:|:--------------------------:|-----------------|
| S1 Linear (homo σ_y=0.5) | ✓ 68%/70% | ✓ 68%/70% | — OLS works; noise is uniform |
| S2 Linear (hetero σ_y=0.1+0.1·x) | ~ 86%/72% | ✓ 71%/72% | Intervals too wide; pooled σ² inflated by noisy high-x points |
| S3 Radioactive decay (Poisson √A) | ✗ 32%/42% | ✓ 64%/68% | Badly overconfident; large early-time counts dominate residuals |
| S4 Power law (8% relative noise) | ✓ 66%/66% | ✓ 68%/69% | — OLS approximately ok here |
| S5 Gaussian peak (Poisson counts) | ✗ 39%/54% | ✓ 66%/70% | Overconfident; amplitude and center poorly constrained |
| S6 Damped oscillator (uniform σ_y) | ✓ 64%/71% | ✓ 67%/72% | — OLS works; noise is uniform |

**The pattern:** OLS coverage is correct only when σ_y is constant across the range (S1, S6). As soon as noise scales with signal — Poisson counting (S3, S5) or percentage-of-reading errors (S2, S4) — the pooled residual variance is a poor proxy for per-point noise, and uncertainty intervals become unreliable. The parameter estimates themselves are often similar; it is the *uncertainty* that OLS gets wrong.

### Using the wrong estimator when x has errors

| Scenario | Wrong estimator | Coverage | Correct estimator | Coverage |
|----------|-----------------|:--------:|-------------------|:--------:|
| Exp decay + timing errors | WeightedRegressor | ✗ 30% | XYWeightedRegressor | ✓ 64% |
| Beer-Lambert | WeightedRegressor | ✗ 7% | XYWeightedRegressor | ✓ 68% |
| Method comparison | WeightedRegressor (OLS) | ✗ 32% | DemingRegressor | ✓ 66% |

See [DEVELOPING.md](DEVELOPING.md) for contributing, running tests, and building docs.
