diff --git a/CHANGELOG.md b/CHANGELOG.md index 0574516..a9a30d7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,11 +15,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 > **[0.0.6]** corresponds to the legacy `weakl` package available on PyPI. --- -## [Unreleased] +## [1.4.1] - 2026-10-08 + +Patch release for probabilistic forecasting: new methods and options only, plus one fix. Existing calls give the same forecasts; the only change of behaviour is that a distributional or mixture `StaticTAM` with the default log target now raises on a training target with values `<= 0` instead of returning a meaningless model (see "Fixed"). - The `DCO` sign-off is required on every commit of a pull request (`DCO`, `CONTRIBUTING.md`, `docs/contributing/dco_howto.md`); a `REUSE` workflow runs `reuse lint` on each push and pull request. - The `pid()` page of `architecture/core/06_the_spectrum_api.md` says when `d_pen` acts and that the derivative is computed along the frame given to `predict()`; three tests pin `pid` behaviour. +### Added +- **Limiting the CPU threads is documented and tested** (`architecture/core/04_hardware_memory.md`, README): set `OMP_NUM_THREADS` and `MKL_NUM_THREADS` before starting Python. The page also describes a hang of `torch.set_num_threads(n)` seen on one PyTorch / oneMKL build, which reproduces without `tam`; the results of `tam` are unchanged. +- **`loss="negative_binomial"` and `loss="tweedie"`** on `StaticTAM` (log link): overdispersed counts, and positive targets with exact zeros (sparse retail, rain, night-time solar). The negative binomial dispersion is estimated by maximum likelihood after the IRLS fit (or fixed with `loss_kwargs={"dispersion": alpha}`); the Tweedie power is `loss_kwargs={"power": p}` with `1 < p < 2` (1.5 by default) and its dispersion is the Pearson estimate. The estimate is `model.dispersion_`. A negative target raises a `ValueError`. +- **`StaticTAM.predict_quantiles(data, taus)` for `poisson`, `negative_binomial` and `tweedie`**: discrete quantiles for the counts and the quantiles of the compound Poisson-gamma law for Tweedie, which are `0` where the probability of zero exceeds the level; non-negative and non-decreasing in the level by construction. The Poisson loss is otherwise unchanged. +- **`KalmanTAM.predict_quantiles(data, taus)`**: Gaussian quantiles of the online forecast from the predictive variance of the filter, `x' P x + R` in the target scale (for `horizon_steps > 1`, the covariance of the delayed state grown by the process noise of the delay). Returns a DataFrame with `q` columns, like the distributional `StaticTAM`; like `predict_online`, `data` holds the target. The constructor, `predict_online`, `predict` and `fit` are unchanged. The intervals are as calibrated as `observation_noise_var` and `process_noise_var` are. +- **`AdaptiveTAM.predict_quantiles(data, taus)`**: online quantile forecasts when the `base_model` is a distributional `StaticTAM` (a dict formula). The adaptive correction moves the location; the scale of the base is rescaled at each update by the root mean square of the standardized residuals of the last `training_window_periods` rows that are at least `horizon_steps` old (past residuals only), and the quantiles go back through the target transform of the base. Returns a DataFrame with `q` columns, like the distributional `StaticTAM`; like `predict_online`, `data` holds the target. A plain base raises a `ValueError`. `predict_online`, `predict` and `fit` give the same results as before for a plain base, and now also accept a distributional one. +- **`dist_kwargs={"target_transform": "asinh" | "logit"}`** on a distributional `StaticTAM`: the location-scale law sits on `arsinh((y - c) / s)` (centre and robust scale estimated on the training target; a skewed, heavy-tailed target with negative values, such as electricity prices) or on the logit of `(y - a) / (b - a)` (a target bounded in `bounds=(a, b)`, estimated from the training range when not given). `predict_median`, `predict_quantile(s)`, `cdf`, `crps`, `anomaly_score` and the conformal intervals work through the transform, and quantiles stay inside the bounds for the logit. `target_transform="log"` is the default and the same as `log_target=True`; `"none"` the same as `log_target=False`. A mixture accepts `"none"` and `"log"` only. +- **`f(x, ..., period=...)`**: the range of one cycle of a Fourier input, on which the feature is normalised in `StaticTAM` and in each rolling window of `AdaptiveTAM`. `period=None` (default): the minimum and maximum of the training rows, as before. `period='auto'`: first value to last value plus one step, the step read from the training rows (hours 0 to 23 give (0, 24), months 1 to 12 give (1, 13)). `period=(low, high)`: given, for a history covering only part of the cycle (`(0, 1)` for a time of year trained on half a year). Without it the training range is still taken as one period, which glues the last value of an integer-coded position onto the first (`f(hour, cyclic=True)` with hours 0 to 23 forces 23:00 = 0:00; Sunday = Monday for 0 to 6; December = January for 1 to 12) and reads a history covering half a year (or a 60-day adaptive window) as a whole year: use `'auto'` for the first case and a given range for the second. Results without `period` are unchanged. +- **`OperaTAM(..., loss_type='pinball', tau=0.9)`**: online aggregation of quantile forecasts with the pinball loss of level `tau` (MLpol and EWA). `tam.model.opera.aggregate_quantiles(df, target_col, expert_cols_by_level)` aggregates every level and sorts them on each row. The square and absolute losses are unchanged. +- **`SafetyTAM.predict_quantiles(y_pred, levels, scale=None, signed=False)`**: conformal quantiles at several levels in one call (`q` columns, sorted on each row). By default the symmetric bands of `predict_intervals`; `signed=True` calibrates each tail on the signed residuals kept by `calibrate()` (one-sided split conformal), for skewed errors. `calibrate(..., verbose=False)` silences the sample-count line, printed by default as before. + +### Fixed +- **A distributional or mixture `StaticTAM` with the default log target rejects a training target with values `<= 0`** (`ValueError` naming the count and `dist_kwargs={"log_target": False}`). They were clipped to a tiny constant before the log, i.e. silently turned into huge negative outliers, and the fitted location, scale and quantiles were meaningless. Positive targets are unchanged. + ## [1.4.0] - 2026-10-06 First release on PyPI and Zenodo since 1.3.1; it also carries everything of 1.3.2 to 1.3.4, which were tags only. **`AdaptiveTAM`, `KalmanTAM` and `NeuralTAM` [EXP] forecasts change** (see 1.3.3, 1.3.4 and the sections below); `StaticTAM` changes only for `te()` with GCV (-2%), a mixture mean (-0.6%) and a prediction without the group column. To get the 1.3.1 behaviour: `pip install tam-ml==1.3.1`. @@ -308,7 +324,8 @@ This version introduced the Formula API and the first object-oriented refactorin ### Added * Initial project setup based on the original `weakl` v0.0.6 package. -[Unreleased]: https://github.com/EDF-Lab/tam/compare/v1.4.0...HEAD +[Unreleased]: https://github.com/EDF-Lab/tam/compare/v1.4.1...HEAD +[1.4.1]: https://github.com/EDF-Lab/tam/releases/tag/v1.4.1 [1.4.0]: https://github.com/EDF-Lab/tam/releases/tag/v1.4.0 [1.3.4]: https://github.com/EDF-Lab/tam/releases/tag/v1.3.4 [1.3.3]: https://github.com/EDF-Lab/tam/releases/tag/v1.3.3 diff --git a/CITATION.cff b/CITATION.cff index a49e7d1..fd16a4a 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -21,8 +21,8 @@ identifiers: repository-code: "https://github.com/EDF-Lab/tam" -version: 1.4.0 -date-released: 2026-10-06 +version: 1.4.1 +date-released: 2026-10-08 keywords: - "Python" diff --git a/README.md b/README.md index 908e8be..d93008c 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ [🧪 `EXAMPLES`](EXAMPLES.md) | [📚 `THEORY`](THEORY.md) | [📄 `PAPER`](paper.md) -**Version:** 1.4.0 +**Version:** 1.4.1 **License:** LGPL-3.0-or-later. **Copyright (c):** 2023-2026 EDF (Electricité De France). **Copyright (c):** 2023-2025 Sorbonne Université. @@ -79,6 +79,8 @@ Raw data flows through dynamic normalizers before being processed by the `Static pip install tam-ml ``` +> **Threads:** to limit the CPU threads, set `OMP_NUM_THREADS` and `MKL_NUM_THREADS` in the environment before starting Python (see `architecture/core/04_hardware_memory.md`). + > ⚠️ **Hardware Note:** TAM relies heavily on PyTorch for VRAM management and tensor acceleration. Ensure you have installed the correct PyTorch distribution (CUDA for NVIDIA, or MPS for Apple Silicon) for your specific hardware to fully leverage the framework's GPU capabilities. --- @@ -167,10 +169,10 @@ If you use these packages in your research, please cite them using their permane **TAM** ```bibtex @misc{tam2026package, - title={TAM: Time series Additive Model (v1.4.0)}, + title={TAM: Time series Additive Model (v1.4.1)}, author={Allioux, Yann and Goude, Yannig}, year={2026}, - doi={10.5281/zenodo.23171729}, + doi={10.5281/zenodo.23248249}, note={With foundational WeaKL research by Nathan Doumèche, and tensorisation and adaptive model contributions by Éloi Bedek} } ``` diff --git a/architecture/core/04_hardware_memory.md b/architecture/core/04_hardware_memory.md index e738d68..5f29efa 100644 --- a/architecture/core/04_hardware_memory.md +++ b/architecture/core/04_hardware_memory.md @@ -24,6 +24,18 @@ Crucially, it manages disaster recovery via the `handle_oom` method. When an ope --- +## Limiting the CPU threads + +To share a machine between several fits (a benchmark, an AutoTAM search run in parallel), limit the threads **through the environment, before Python starts**: + +```bash +OMP_NUM_THREADS=2 MKL_NUM_THREADS=2 python my_script.py +``` + +On Windows PowerShell: `$env:OMP_NUM_THREADS=2; $env:MKL_NUM_THREADS=2; python my_script.py`. `torch.get_num_threads()` then returns 2 (`tests/test_thread_limits.py` runs a grouped fit this way in a fresh process). + +**Known issue on one build.** With the PyTorch 2.14.0 CPU wheel on Windows (Intel oneMKL 2026.1), calling `torch.set_num_threads(n)` with `n >= 2` made every later batched LU factorisation of matrices of a few hundred rows (`torch.linalg.solve`, `lu_factor`, `inv`, which the grouped solves use) print `Intel oneMKL ERROR: Parameter 6 was incorrect on entry to DLASWP` and hang. It reproduces with PyTorch alone, without `tam`, so it is not specific to this library (Cholesky factorisations, unbatched solves and `n = 1` were not affected). Other platforms and builds were not tested. If a fit prints that message on your setup, remove the `torch.set_num_threads` call and use the environment variables above. + ## The Memory Oracle and Safe Chunking To proactively avoid invoking the OOM handler, the framework utilizes `_memory.py` as an advanced predictive oracle. diff --git a/architecture/core/06_the_spectrum_api.md b/architecture/core/06_the_spectrum_api.md index 2aea02f..aee4a93 100644 --- a/architecture/core/06_the_spectrum_api.md +++ b/architecture/core/06_the_spectrum_api.md @@ -65,6 +65,19 @@ The linear effect applies direct spatial scaling, and its penalty is a standard :end-before: "#: " ``` +**Declared period.** `f(x, ..., period=...)` makes the feature normalised on the range of one cycle instead of the minimum and +maximum of the training rows: `None` (default) keeps the training range, `'auto'` uses (first value, last value + one step), the +step being the smallest gap between the distinct training values (resolved once on all training rows by `_resolve_fixed_ranges`, +so every group and adaptive window shares it), and `(low, high)` is used as given (a history covering part of the cycle). This applies in `StaticTAM` and in every rolling window +of `AdaptiveTAM` (`KalmanTAM` fits through a `StaticTAM`). Without it the training range is the period, as before: an integer-coded position +has its last value glued onto its first (hours 0 to 23 give 23:00 = 0:00; use `period='auto'` or `(0, 24)`), and a history covering half a cycle is read as a whole cycle. `fixed_ranges` collects the declared periods (tensor products included) and +`_fit_normalization_params` / `_transform_data_adaptive` apply them; two different periods for one feature raise. The default +`extrapolate='continue'` lets an input beyond the period wrap around the cycle (`'constant'` would clamp it). +```{literalinclude} ../../../../src/tam/model/spectrum/_factory.py +:language: python +:pyobject: fixed_ranges +``` + **2. Feature Map** Evaluates batched trigonometric functions, dynamically scaling angular frequencies based on the `cyclic` boundary flag to prevent endpoint distortions. ```{literalinclude} ../../../../src/tam/model/spectrum/_fourier.py diff --git a/architecture/meta/01_adaptive_code.md b/architecture/meta/01_adaptive_code.md index d60687e..b2e0daa 100644 --- a/architecture/meta/01_adaptive_code.md +++ b/architecture/meta/01_adaptive_code.md @@ -96,6 +96,28 @@ The `grid_search_fit()` method implements the "Multi-Start" hyperparameter solve :caption: src/tam/model/adaptative.py (Coordinate Descent Algorithm for Hyperparameters) ``` +## Online Quantiles from a Distributional Base + +`AdaptiveTAM(base_model=...)` also accepts a fitted distributional `StaticTAM` (a dict formula). The correction is then fitted on the model scale of its location sub-model (the target `__mu__`, the target after its transform): the constructor redirects the base to that sub-model, maps `Residual` in `adaptive_formula` to `Residual__mu__`, and `_with_model_target` adds the model-scale column to any frame that holds the target. The forecasts of a plain base are untouched. + +`predict_quantiles(data, taus)` is the probabilistic counterpart of `predict_online`, so `data` holds the target. It runs the online simulation, then rescales the scale of the base with `_online_scale`, a rolling root mean square of the standardized residuals computed per group from cumulative sums (vectorised), and goes back to the response scale with the transform of the base: + +```{literalinclude} ../../../../src/tam/model/adaptative.py +:language: python +:start-after: "#: " +:end-before: "#: " +:caption: src/tam/model/adaptative.py (rolling scale of the standardized residuals) +``` + +```{literalinclude} ../../../../src/tam/model/adaptative.py +:language: python +:start-after: "#: " +:end-before: "#: " +:caption: src/tam/model/adaptative.py (online quantiles) +``` + +The first window forecasts the row `training_window_periods + horizon_steps - 1` of a group, and the scale is renewed every `update_interval_periods` rows; each renewal reads the residuals of the `training_window_periods` rows that are at least `horizon_steps` old, so no residual of a later row enters an earlier quantile. A renewal with fewer than five usable residuals keeps the scale of the base. A plain base raises a `ValueError`. + ## Separation of Concerns: Simulation vs. Inference To guarantee speed and safety in operational production pipelines, `AdaptiveTAM` strictly separates the continuous historical simulation from out-of-sample inference. diff --git a/architecture/meta/02_kalman_torchscript.md b/architecture/meta/02_kalman_torchscript.md index ac2c8c6..d9023c8 100644 --- a/architecture/meta/02_kalman_torchscript.md +++ b/architecture/meta/02_kalman_torchscript.md @@ -18,3 +18,16 @@ **Process noise per term (verified by `tests/test_kalman_per_term_noise.py`):** the TorchScript kernel takes the diagonal of $Q$ as a tensor (`process_noise_diag`, one variance per design column) instead of a scalar. `KalmanTAM._term_spans` maps each formula term (the names of `decomposition_names`, plus `"offset"` for the two constant columns: the filter's systematic offset and the formula's intercept) to its design columns, and `_noise_diagonal` fills the diagonal from `process_noise_var` (a float, or a dict `{term: q, "offset": q, "default": q}`). The boost on the offset's noise (`offset_q_boost`) is `offset_boost` unless the offset has a variance of its own; the initial covariance of the offset keeps `offset_boost`. A float call builds the same diagonal as before, so its forecasts are unchanged. **Base-model components (verified by `tests/test_base_effect_columns.py`):** with `add_base_effects=True`, one `l(effect_)` term per base-model component is appended to the Kalman formula, with the column names `decompose_prediction` emits (`decomposition_names()`). + +## Predictive quantiles + +The compiled loop can also return the predictive variance of every forecast row (`with_variance`, off by default: the filter and its outputs are then identical to before). Per block it adds one batched product `X P Xᵀ` on the a priori covariance, and for `horizon_steps > 1` the covariance of the delayed state grown by `(h - 1) Q`, written at the row it forecasts. `_run_filter` keeps it in `_predictive_variance_`; `_online` is the former body of `predict_online`, which now calls it, and `predict_quantiles` converts the standard deviation to the target scale with the reference scale and adds `z_tau` standard deviations to the online forecast: + +```{literalinclude} ../../../../src/tam/model/kalman.py +:language: python +:start-after: "#: " +:end-before: "#: " +:caption: src/tam/model/kalman.py (Gaussian quantiles from the filter variance) +``` + +The constructor and the outputs of `predict_online`, `predict` and `fit` are unchanged. One limit of the filter carries over: a row without a target leaves the state unchanged but still reduces the covariance, so the variance after a long gap of missing targets is understated. diff --git a/architecture/meta/05_opera_gpu.md b/architecture/meta/05_opera_gpu.md index 205e36a..e56a63e 100644 --- a/architecture/meta/05_opera_gpu.md +++ b/architecture/meta/05_opera_gpu.md @@ -60,6 +60,16 @@ A `torch.where` mask safely intercepts any groups with perfectly zero targets to :caption: src/tam/model/opera.py (TorchScript Compiled EWA Loop) ``` +`loss_type='pinball'` with `tau` aggregates quantile forecasts: the MLpol loop uses the subgradient `1{y_hat > y} - tau` in place +of the square or absolute one, the EWA loop the pinball loss of each expert. `aggregate_quantiles` runs one `OperaTAM` per level (each +with its own expert columns) and sorts the levels on each row: + +```{literalinclude} ../../../../src/tam/model/opera.py +:language: python +:start-after: "#: " +:end-before: "#: " +``` + --- ## Causal Boundary Enforcement (Horizon Shifting) diff --git a/architecture/statistics/01_reweighted_estimation_code.md b/architecture/statistics/01_reweighted_estimation_code.md index 301bcb4..d0b9411 100644 --- a/architecture/statistics/01_reweighted_estimation_code.md +++ b/architecture/statistics/01_reweighted_estimation_code.md @@ -16,6 +16,28 @@ The `estimation/` package mirrors `spectrum/`: a base contract, one file per los :end-before: "#: " ``` +## Negative binomial and Tweedie + +`negative_binomial` and `tweedie` are `GLMFamily` subclasses in `_count_families.py`, registered in `build_strategy` (`StaticTAM(..., loss="negative_binomial", loss_kwargs={"dispersion": 0.5})` fixes the dispersion, `loss_kwargs={"power": 1.3}` sets the Tweedie power). The negative binomial variance reads the current `dispersion`, so the same IRLS loop serves any value of it: + +```{literalinclude} ../../../../src/tam/model/statistics/estimation/_count_families.py +:language: python +:start-after: "#: " +:end-before: "#: " +``` + +After the first fit, `BaseTAM.fit` asks a family that has `estimate_dispersion` for the dispersion at the fitted means and, if it moved by more than 0.1%, refits the means with the new variance (at most ten rounds). The Tweedie dispersion does not enter the weights, so its estimate returns without a refit. `validate_target` rejects a negative target before the fit. The estimate is exposed as `StaticTAM.dispersion_` (None for the other losses). + +The Tweedie family computes the compound Poisson-gamma distribution function with the Poisson jump probabilities up to a tail of 1e-10, and its quantiles by bisection on that function (cost proportional to rows times jumps times 45 steps; a quantile below the mass at zero is exactly 0): + +```{literalinclude} ../../../../src/tam/model/statistics/estimation/_count_families.py +:language: python +:start-after: "#: " +:end-before: "#: " +``` + +`StaticTAM.predict_quantiles(data, taus)` on a plain model whose loss has a `quantiles` method (`poisson`, `negative_binomial`, `tweedie`) returns the `q` columns of `predict_count_quantiles`; on any other plain loss it raises as before. + ## Robust M-estimators and expectiles ```{literalinclude} ../../../../src/tam/model/statistics/estimation/_robust.py diff --git a/architecture/statistics/02_distributional_code.md b/architecture/statistics/02_distributional_code.md index 8cc7fb5..1f02fed 100644 --- a/architecture/statistics/02_distributional_code.md +++ b/architecture/statistics/02_distributional_code.md @@ -36,6 +36,37 @@ never numpy's silent overflow warning: :pyobject: _to_response_scale ``` +The log target is the default (`dist_kwargs={"log_target": True}`). `fit` first checks the training target: a value `<= 0` has no +logarithm, and clipping it to a tiny constant would turn it into a huge negative outlier, so the fit raises a `ValueError` naming the +count and the option `log_target=False` instead. The mixture mode (`_mixture.py`) runs the same check. Only the training target is +checked: the CDF of a new observation `<= 0` under a log-scale law is ~0, a valid answer. + +```{literalinclude} ../../../../src/tam/model/statistics/estimation/_distributional.py +:language: python +:pyobject: require_positive_target +``` + +**Target transform.** `dist_kwargs={"target_transform": ...}` takes `"none"`, `"log"` (the default, same as `log_target=True`), `"asinh"` or `"logit"`; giving it together with a contradicting `log_target` raises a `ValueError`, and a mixture accepts `"none"` and `"log"` only. `resolve_target_transform` reads the options once, in the constructor: + +```{literalinclude} ../../../../src/tam/model/statistics/estimation/_distributional.py +:language: python +:pyobject: resolve_target_transform +``` + +`fit` estimates the parameters of the transform on the training target and stores them in `_transform_params_`: for `"asinh"` the median and a robust scale (normalised MAD, else the standard deviation); for `"logit"` the bounds, from `dist_kwargs={"bounds": (low, high)}` or the observed range widened by 5% on each side, and a training value outside given bounds raises a `ValueError`. + +```{literalinclude} ../../../../src/tam/model/statistics/estimation/_distributional.py +:language: python +:pyobject: fit_target_transform +``` + +The forward map is applied to the target everywhere it is read (fit, CDF, CRPS, anomaly score), and `from_model_scale` is the single back-transform used by `predict_median`, `predict_quantile` and `anomaly_score`, and by the studentized conformal intervals through the model method `_from_model_scale`: + +```{literalinclude} ../../../../src/tam/model/statistics/estimation/_distributional.py +:language: python +:pyobject: from_model_scale +``` + ## The `additive.py` bridge (one-line delegates) ```{literalinclude} ../../../../src/tam/model/additive.py diff --git a/architecture/statistics/04_risk_conformal_evt_code.md b/architecture/statistics/04_risk_conformal_evt_code.md index 927ad23..bd6a4ed 100644 --- a/architecture/statistics/04_risk_conformal_evt_code.md +++ b/architecture/statistics/04_risk_conformal_evt_code.md @@ -15,6 +15,16 @@ One source of truth for conformal: the static engine `safety.py::SafetyTAM` hold :pyobject: SafetyTAM.conformal_quantile ``` +`calibrate()` keeps the signed residuals (`signed_scores_`, studentized when a scale is given) next to the absolute scores, and prints the +sample count unless `verbose=False`. `predict_quantiles()` returns several levels at once, `q` columns: symmetric bands from +`conformal_quantile(2 * tau)` by default (the bands of `predict_intervals`), or one-sided bounds per tail from the signed scores with +`signed=True`; the levels are sorted on each row. + +```{literalinclude} ../../../../src/tam/model/safety.py +:language: python +:pyobject: SafetyTAM.predict_quantiles +``` + ## Adaptive Conformal Inference (streaming, model-agnostic) ```{literalinclude} ../../../../src/tam/model/statistics/risk/aci.py diff --git a/docs/source/conf.py b/docs/source/conf.py index a02b878..7429745 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -20,7 +20,7 @@ ) author = 'Yann Allioux' -release = '1.4.0' +release = '1.4.1' # -- General Configuration --------------------------------------------------- extensions = [ diff --git a/math/meta/01_adaptive_online.md b/math/meta/01_adaptive_online.md index fe160f8..0bc1503 100644 --- a/math/meta/01_adaptive_online.md +++ b/math/meta/01_adaptive_online.md @@ -53,6 +53,18 @@ $$\hat{\theta}_{adapt, t} = \left( \Phi_{adapt, W}^\top \Lambda_{W}^\top \Lambda By flattening the group and window dimensions into a single massive 3D tensor batch, the framework computes these dynamic corrector weights simultaneously across all groups and time steps. This isolates the regularization scale for each distinct topological window and neutralizes concept drift, bypassing the computational overhead typically associated with rolling dense matrix inversions. +### Online Scale for Quantile Forecasts + +When the long-term expert is a location-scale law with location $\mu_0$ and scale $\sigma_0$, the corrected location $\hat\mu_t=\mu_0+\hat\epsilon_t$ moves with the data but the scale does not, so the intervals do not follow a change of variance. The standardized residual $r_s=(y_s-\hat\mu_s)/\sigma_{0,s}$ has unit scale when the expert is calibrated; its root mean square over the most recent training window, + +$$ s_t = \sqrt{\frac{1}{|\mathcal{W}_t|}\sum_{u\in\mathcal{W}_t} r_u^2}, $$ + +measures the miscalibration, and the quantile of level $\tau$ is + +$$ Q_\tau(t) = \hat\mu_t + s_t\,\sigma_{0,t}\,F^{-1}(\tau), $$ + +with $F$ the standardized law of the expert. The window $\mathcal{W}_t$ holds only residuals observed before the update that serves $t$, at least one forecast horizon back, so the construction is causal. After a break from variance $1$ to variance $9$, $s_t$ rises to $3$ within one window and the interval widens by the same factor, which restores the nominal coverage. + ### Equivalence with a Rolling StaticTAM Each window is normalised with the minimum and maximum of its own training rows only, and no bound fitted on the whole period is applied. The forecasts of a window are therefore those of the estimator fitted on that window's training rows and applied to its forecast rows: the adaptive model is the operational loop "refit on the last window, forecast the next rows", computed in one batch. A feature whose range moves between windows is handled as a refit would handle it. \ No newline at end of file diff --git a/math/meta/02_kalman_filter.md b/math/meta/02_kalman_filter.md index 00e72e1..44d33d1 100644 --- a/math/meta/02_kalman_filter.md +++ b/math/meta/02_kalman_filter.md @@ -31,3 +31,15 @@ The tracked coefficients follow a random walk $\theta_{t+1} = \theta_t + w_t$, $ $$Q = \mathrm{diag}\big(q_{\text{offset}},\; q_{1}\,\mathbf{1}_{d_1},\; q_{2}\,\mathbf{1}_{d_2},\; \dots\big)$$ All the $d_k$ design columns of term $k$ share its variance $q_k$. Variational Bayesian variance tracking {cite:p}`vilmarest2024viking` estimates such variances dynamically, online; the variances of `KalmanTAM` stay static: set by the user or tuned offline on a calibration period. A variance of $0$ lets the covariance of the term shrink with every observation, so the term follows a changing coefficient less and less. The variances are best chosen on a calibration period that ends before the period that is reported. + +## Predictive variance and quantiles + +The covariance recursion of the filter does not depend on the observations: it is a function of the design rows $x_t$, the observation noise $R$ and the process noise $Q$. Before the update at step $t$, the state has mean $\hat\theta_t$ and covariance $P_t$, so the one-step forecast of the standardised target is Gaussian, + +$$ y_t \mid y_{y\}-\tau$. +At $\tau=\tfrac12$ it is half the subgradient of the absolute loss, and the adaptive learning rates below make MLpol insensitive to a +constant factor on the regrets. Each level is aggregated on its own {cite:p}`gaillard2016opera`; the aggregated levels are then +sorted on each time step so that the quantile forecasts are monotone in $\tau$ (they never cross). + According to the polynomial potential bounds established by Cesa-Bianchi and Lugosi {cite:p}`cesa2006prediction`, the optimal weight update is proportional to the gradient of the squared positive regret. The TAM implementation strictly translates this as applying a Rectified Linear Unit (ReLU) to the cumulative regret, multiplied by an adaptive learning rate $\eta_{k,t}$: $$w_{k,t+1} = \frac{\eta_{k,t} \max(0, R_{k,t})}{\sum_{j=1}^K \eta_{j,t} \max(0, R_{j,t})}$$ diff --git a/math/spectrum/FOURIER.md b/math/spectrum/FOURIER.md index 1bdcb90..c051105 100644 --- a/math/spectrum/FOURIER.md +++ b/math/spectrum/FOURIER.md @@ -14,7 +14,13 @@ To achieve Gigadata scalability while strictly preserving the theoretical guaran ## Formula Definition: The Real Spectral Basis ($\Phi$) -As established in the data preparation phase, the input variable $x$ is normalized to the domain $[-1, 1]$. To properly evaluate the trigonometric basis, this domain is linearly rescaled: +As established in the data preparation phase, the input variable $x$ is normalized to the domain $[-1, 1]$. By default the normalization +maps the training range $[\min x, \max x]$ onto $[-1, 1]$; for a periodic phenomenon this takes the training rows to span exactly one +period. For a position sampled on a regular grid $x_00$: $\log y$ is undefined at and + below zero, so a target taking such values is modelled on its own scale, $t=y$. 2. **Scale.** Read the squared residuals $r^2=(t-\hat\mu)^2$ and fit $\sigma(x)$ as a **Gamma-GLM** atom on $r^2$ (or an L2 atom on $\log r^2$ with a $\log\chi^2_1$ bias correction). ## The tail law and quantiles @@ -18,4 +19,15 @@ The standardized residual law $F$ (Normal or Student-t) is selected from the res $$ Q_\tau(x) = \exp\!\Big(\hat\mu(x) + \hat\sigma(x)\,F^{-1}(\tau)\Big). $$ +## Monotone target transforms + +The location-scale law can sit on a transformed target $t=g(y)$ instead of the logarithm, for any strictly increasing $g$. Quantiles of $y$ follow exactly, since a monotone map preserves order: + +$$ Q_\tau^{(y)}(x) = g^{-1}\Big(\hat\mu(x) + \hat\sigma(x)\,F^{-1}(\tau)\Big). $$ + +Two maps cover the targets that the logarithm cannot take. + +- **Real-valued, skewed, heavy-tailed targets with negative values** (electricity prices, spikes and returns): the area hyperbolic sine $g(y)=\operatorname{arsinh}\big((y-c)/s\big)$, with the centre $c$ the median and the scale $s$ a robust spread of the training target. It behaves like the identity near $c$ and like a logarithm in both tails, and it is the standard variance-stabilizing choice for electricity prices {cite:p}`uniejewski2018variance`. With a Normal law on $t$, $y$ follows the Johnson $S_U$ family {cite:p}`johnson1949systems`: skewed, with heavier tails than a Normal. +- **Targets bounded in an interval $(a,b)$** (shares, humidity, capacity-capped power): the logit $g(y)=\log\big(u/(1-u)\big)$ with $u=(y-a)/(b-a)$. Every quantile lies strictly inside $(a,b)$, which a Normal law on the response scale cannot guarantee. A value exactly on a bound has an infinite logit and is moved inside by a tiny fraction of the interval. + An optional convex `scale_shrinkage` $\lambda\in[0,1]$ blends the conditional variance toward the global one. The probability integral transform $F\big((t-\hat\mu)/\hat\sigma\big)$ gives the CDF, a two-sided anomaly score, and the CRPS (closed form for the Normal, else a Gauss-Legendre quantile integral). diff --git a/math/statistics/04_risk_conformal_evt.md b/math/statistics/04_risk_conformal_evt.md index 14ada92..14a813d 100644 --- a/math/statistics/04_risk_conformal_evt.md +++ b/math/statistics/04_risk_conformal_evt.md @@ -23,6 +23,13 @@ The framework implements the **Split Conformal** method {cite:p}`lei2018distribu $$ \hat q = \mathrm{Quantile}\!\left(\{s_i\},\ (1-\alpha)\big(1+\tfrac1n\big)\right) $$ 4. **Prediction Interval:** For any new observation at time $t$, the prediction interval is symmetrically defined around the baseline by $\hat{q}$. +**Several quantile levels and skewed errors.** A quantile at level $\tau<\tfrac12$ is the lower end of the symmetric interval with $\alpha=2\tau$, +and $1-\tau$ its upper end. With skewed errors one radius cannot fit both tails: the band is too wide on the short side and too narrow on the +long one. Keeping the signed scores $r_i=(Y_i-\hat Y_i)/\hat s_i$ and calibrating each tail on its own gives one-sided split-conformal +bounds: with $r_{(1)}\le\dots\le r_{(n)}$, the level-$\tau$ quantile is $\hat Y+\hat s\,r_{(k)}$ with $k=\lceil (n+1)\tau\rceil$ for $\tau\ge\tfrac12$ +and $k=\lfloor (n+1)\tau\rfloor$ for $\tau<\tfrac12$, so that $P(Y\le \hat Q_\tau)\ge\tau$ above the median and $P(Y<\hat Q_\tau)\le\tau$ below it +under exchangeability {cite:p}`lei2018distribution`. + --- ## Adaptive Conformal Inference (ACI) diff --git a/references.bib b/references.bib index a4bc0c2..675bf3e 100644 --- a/references.bib +++ b/references.bib @@ -1178,4 +1178,44 @@ @article{green1984iteratively pages={149--170}, year={1984}, publisher={Wiley Online Library} -} \ No newline at end of file +} + +@article{uniejewski2018variance, + title = {Variance stabilizing transformations for electricity spot price forecasting}, + author = {Uniejewski, Bartosz and Weron, Rafa{\l} and Ziel, Florian}, + journal = {IEEE Transactions on Power Systems}, + volume = {33}, + number = {2}, + pages = {2219--2229}, + year = {2018} +} + +@article{johnson1949systems, + title = {Systems of frequency curves generated by methods of translation}, + author = {Johnson, Norman L.}, + journal = {Biometrika}, + volume = {36}, + number = {1/2}, + pages = {149--176}, + year = {1949} +} + +@article{jorgensen1987exponential, + title = {Exponential dispersion models}, + author = {J{\o}rgensen, Bent}, + journal = {Journal of the Royal Statistical Society: Series B}, + volume = {49}, + number = {2}, + pages = {127--162}, + year = {1987} +} + +@article{lawless1987negative, + title = {Negative binomial and mixed {P}oisson regression}, + author = {Lawless, Jerald F.}, + journal = {The Canadian Journal of Statistics}, + volume = {15}, + number = {3}, + pages = {209--225}, + year = {1987} +} diff --git a/src/tam/__init__.py b/src/tam/__init__.py index 3a18934..05ea6b6 100644 --- a/src/tam/__init__.py +++ b/src/tam/__init__.py @@ -27,7 +27,7 @@ from .evaluation.tracker import BenchmarkTracker from .common.plotting import plot_component -__version__ = "1.4.0" +__version__ = "1.4.1" __all__ = [ "StaticTAM", diff --git a/src/tam/model/_base.py b/src/tam/model/_base.py index cd41368..ee9f68e 100644 --- a/src/tam/model/_base.py +++ b/src/tam/model/_base.py @@ -24,6 +24,9 @@ from ._data import _reassemble_predictions, _check_known_groups, _coefficients_of_groups from ._dispatcher import smart_solve +# Outer refits of the means when a family estimates a dispersion (negative binomial): it converges in a few. +_DISPERSION_REFITS = 10 + #: class BaseTAM(ABC): r""" @@ -247,12 +250,23 @@ def fit(self, data_train: pd.DataFrame) -> 'BaseTAM': if reweighting_strategy is None: self.coefficients_ = self._solve_pwls_step(x_train, y_train, weights=None) else: - from .statistics.estimation._reweighting import reweighted_penalized_fit - self.coefficients_ = reweighted_penalized_fit( - self, x_train, y_train, reweighting_strategy, - max_iter=getattr(self, "_irls_max_iter_", 25), - tol=getattr(self, "_irls_tol_", 1e-6), - ) + from .statistics.estimation._reweighting import reweighted_penalized_fit, _linear_predictor + validate = getattr(reweighting_strategy, "validate_target", None) + if validate is not None: + validate(y_train) + irls = dict(max_iter=getattr(self, "_irls_max_iter_", 25), tol=getattr(self, "_irls_tol_", 1e-6)) + self.coefficients_ = reweighted_penalized_fit(self, x_train, y_train, reweighting_strategy, **irls) + # A family with a dispersion (negative binomial, Tweedie) estimates it from the fitted means and refits the means + # with the new variance until it stops moving. + estimate = getattr(reweighting_strategy, "estimate_dispersion", None) + if estimate is not None: + y_fit = y_train.to(torch.get_default_dtype()) + for _ in range(_DISPERSION_REFITS): + mean = reweighting_strategy.inverse_link(_linear_predictor(self, x_train, self.coefficients_)) + if not estimate(y_fit, mean): + break + self.coefficients_ = reweighted_penalized_fit(self, x_train, y_train, reweighting_strategy, **irls) + self.dispersion_ = getattr(reweighting_strategy, "dispersion", None) return self #: diff --git a/src/tam/model/_data.py b/src/tam/model/_data.py index b44bc14..c646905 100644 --- a/src/tam/model/_data.py +++ b/src/tam/model/_data.py @@ -28,7 +28,8 @@ def _fit_normalization_params( data: pd.DataFrame, features: List[str], group_col: str, - categorical_levels: Optional[Dict[str, int]] = None + categorical_levels: Optional[Dict[str, int]] = None, + fixed_ranges: Optional[Dict[str, Tuple[float, float]]] = None ) -> Tuple[Dict, List]: r""" Calculates the min/max normalization parameters for features, computed per group. @@ -39,6 +40,8 @@ def _fit_normalization_params( group_col: The column name used to group the data. categorical_levels: ``{feature: n_cat}`` of the categorical features: their range is the full level range (``categorical_range``), whatever levels the training rows hold. + fixed_ranges: ``{feature: (low, high)}`` of the features read on a fixed domain (a cyclic Fourier term reads [-1, 1]): + that range is used whatever values the training rows hold. Returns: A tuple (norm_params, unique_groups): @@ -63,9 +66,41 @@ def _fit_normalization_params( column = grouped.get_group(group_name)[feature] params['min'][feature], params['max'][feature] = categorical_range( n_cat, params['min'][feature], params['max'][feature], bool((column.dropna() == np.round(column.dropna())).all())) - + for feature, (low, high) in (fixed_ranges or {}).items(): + if feature in params['min'].index: + params['min'][feature], params['max'][feature] = float(low), float(high) + return norm_params, unique_groups +def _resolve_fixed_ranges(data: pd.DataFrame, ranges: Optional[Dict]) -> Dict[str, Tuple[float, float]]: + r""" + Turns the ``"auto"`` periods of ``fixed_ranges`` into ``(low, high)`` read from the training rows; given ranges pass through. + + ``"auto"`` is (first value, last value + one step), the step being the smallest gap between the distinct values (rounded + to 10 decimals): hours 0..23 give (0, 24), months 1..12 give (1, 13), a time of year sampled once a day gives (0, 1). + It is computed once on all the rows given (all groups, the whole series), so every group and every window shares it. + A feature with a single distinct value keeps the min/max rule. + + Args: + data: The training rows (or the whole series of an adaptive simulation). + ranges: ``{feature: (low, high) or "auto"}``. + + Returns: + ``{feature: (low, high)}``. + """ + resolved: Dict[str, Tuple[float, float]] = {} + for feature, value in (ranges or {}).items(): + if value != "auto": + resolved[feature] = (float(value[0]), float(value[1])) + continue + if feature not in data.columns: + continue + distinct = np.unique(np.round(data[feature].dropna().to_numpy(dtype=float), 10)) + if distinct.size < 2: + continue + resolved[feature] = (float(distinct[0]), float(distinct[-1] + np.min(np.diff(distinct)))) + return resolved + #: def normalize(df_to_normalize: pd.DataFrame, params: dict) -> pd.DataFrame: @@ -374,7 +409,8 @@ def _transform_data_adaptive( steps_per_period: int, horizon_steps: int = 1, date_col: Optional[str] = None, - categorical_levels: Optional[Dict[str, int]] = None + categorical_levels: Optional[Dict[str, int]] = None, + fixed_ranges: Optional[Dict[str, Tuple[float, float]]] = None ) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor, _AdaptiveWindows]: r""" Cuts the data into the windows of a rolling refit, and normalises each window on its own training rows. @@ -395,6 +431,7 @@ def _transform_data_adaptive( steps_per_period: Rows of one group in a period. horizon_steps: Horizon of forecasting. categorical_levels: ``{feature: n_cat}`` of the categorical features: normalised on their full level range. + fixed_ranges: ``{feature: (low, high)}`` of the features read on a fixed domain: every window uses that range. Returns: (x_stacked, y_stacked, x_to_predict, windows): (n_groups, n_windows, L, F), (n_groups, n_windows, L, 1), @@ -407,6 +444,7 @@ def _transform_data_adaptive( dtype = torch.get_default_dtype() n_features = len(features) categorical = categorical_levels or {} + fixed = fixed_ranges or {} x_all = data[features].to_numpy(dtype=np.float64) y_all = data[target_col].to_numpy(dtype=np.float64) @@ -444,6 +482,8 @@ def _transform_data_adaptive( start = torch.where(high[..., j] <= n_cat - 1, torch.zeros_like(low[..., j]), low[..., j]) low[..., j] = torch.where(integer, start, low[..., j]) high[..., j] = torch.where(integer, start + (n_cat - 1), high[..., j]) + elif name in fixed: + low[..., j], high[..., j] = fixed[name] amplitude = high - low amplitude = torch.where(amplitude == 0, torch.ones_like(amplitude), amplitude) center = (high + low) / 2 diff --git a/src/tam/model/adaptative.py b/src/tam/model/adaptative.py index bd12ecd..9e7fde5 100644 --- a/src/tam/model/adaptative.py +++ b/src/tam/model/adaptative.py @@ -19,7 +19,8 @@ by mapping the base model's decomposed effects to its current error. """ -from typing import Dict, List, Any, Union, Tuple, Optional +import re +from typing import Dict, List, Any, Sequence, Union, Tuple, Optional import torch import pandas as pd import numpy as np @@ -35,6 +36,7 @@ from ._memory import get_safe_window_batch_size from ._data import ( + _resolve_fixed_ranges, _transform_data_adaptive, _transform_data_stacked, _window_starts, @@ -52,6 +54,7 @@ BaseEffect, create_effects_from_parsed_terms, categorical_ranges, + fixed_ranges, categorical_features, extrapolating_features, initialize_effects, @@ -116,6 +119,56 @@ def rolling_windows( #: +_MODEL_TARGET = "__mu__" +# A scale needs a few residuals: below this count the base scale stands (s = 1). +_MIN_RESIDUALS_FOR_SCALE = 5 + + +#: +def _online_scale(standardized: np.ndarray, group: np.ndarray, order: np.ndarray, window: int, interval: int, horizon: int) -> np.ndarray: + r""" + Rolling root mean square of the standardized residuals, one value per row, using past residuals only. + + Rows of a group are taken in the order of ``order``. The first window forecasts the row ``window + horizon - 1`` (zero-based); + from there the scale is renewed every ``interval`` rows, at an origin ``o``, from the finite residuals of the rows + ``o - window ... o - horizon``. A row before the first origin, or whose origin has fewer than ``_MIN_RESIDUALS_FOR_SCALE`` + finite residuals, gets 1. + + Args: + standardized: ``(y - location) / sigma_base`` per row (NaN where unknown). + group: Group label per row. + order: Sort key within a group (the dates). + window: ``training_window_periods``. + interval: ``update_interval_periods``. + horizon: ``horizon_steps``. + + Returns: + The scale per row, in the row order of the input. + """ + scale = np.ones(len(standardized)) + first = window + horizon - 1 + frame = pd.DataFrame({"g": group, "o": order, "pos": np.arange(len(standardized))}).sort_values(["g", "o", "pos"], kind="stable") + for _, rows in frame.groupby("g", sort=False): + positions = rows["pos"].to_numpy() + residual = standardized[positions] + finite = np.isfinite(residual) + squares = np.where(finite, residual ** 2, 0.0) + cum_squares = np.concatenate([[0.0], np.cumsum(squares)]) + cum_count = np.concatenate([[0], np.cumsum(finite)]) + rank = np.arange(len(positions)) + eligible = rank >= first + origin = first + ((rank - first) // interval) * interval + low = np.maximum(origin - window, 0) + high = origin - horizon + 1 # exclusive end of the rows known at the origin + usable = eligible & (high > low) + count = np.where(usable, cum_count[np.clip(high, 0, None)] - cum_count[low], 0) + total = np.where(usable, cum_squares[np.clip(high, 0, None)] - cum_squares[low], 0.0) + enough = usable & (count >= _MIN_RESIDUALS_FOR_SCALE) + scale[positions] = np.where(enough, np.sqrt(total / np.maximum(count, 1)), 1.0) + return scale +#: + + class AdaptiveTAM: r""" Initializes the AdaptiveTAM model. @@ -139,7 +192,8 @@ def __init__( Initializes the AdaptiveTAM model. Args: - base_model: A fitted StaticTAM model instance. + base_model: A fitted StaticTAM model instance. A fitted distributional StaticTAM (dict formula) is accepted: the correction is + fitted on its location on the model scale, and ``predict_quantiles`` rescales its scale online. adaptive_formula: Formula for the adaptive correction model. Features must be columns produced by base_model.decompose_prediction() (e.g., 'effect_temp'). @@ -161,6 +215,16 @@ def __init__( Raises: ValueError: If the base_model has not been fitted. """ + self.distributional_base_ = None + if base_model is not None and getattr(base_model, '_mode_', 'plain') == 'distributional': + # The correction is fitted on the model scale of the location sub-model (target ``__mu__``); the scale of the + # distributional base is kept for ``predict_quantiles``. + if getattr(base_model, '_location_submodel_', None) is None: + raise ValueError("The base_model must be fitted before initializing AdaptiveTAM.") + self.distributional_base_ = base_model + adaptive_formula = re.sub(rf"\bResidual{re.escape(base_model.target_col_)}\b", f"Residual{_MODEL_TARGET}", adaptive_formula) + base_model = base_model._location_submodel_ + if base_model is not None and getattr(base_model, 'coefficients_', None) is None: raise ValueError("The base_model must be fitted before initializing AdaptiveTAM.") @@ -217,6 +281,15 @@ def __init__( #: #: + def _with_model_target(self, data: pd.DataFrame) -> pd.DataFrame: + r"""Adds the model-scale target column of a distributional base (``__mu__``) when ``data`` holds the target.""" + base = self.distributional_base_ + if base is None or base.target_col_ not in data.columns: + return data + data = data.copy() + data[_MODEL_TARGET] = base._to_model_scale(data[base.target_col_].to_numpy()) + return data + def prepare_simulation(self, data: pd.DataFrame) -> 'AdaptiveTAM': r""" Prepares tensors for the adaptive sliding-window simulation. @@ -233,6 +306,7 @@ def prepare_simulation(self, data: pd.DataFrame) -> 'AdaptiveTAM': self: The instance with populated ``simulation_data_``. """ + data = self._with_model_target(data) if self.base_model_ is not None: data_bm = self.base_model_.decompose_prediction(data) data_pred = self.base_model_.predict(data) @@ -293,7 +367,8 @@ def prepare_simulation(self, data: pd.DataFrame) -> 'AdaptiveTAM': steps_per_period=self.steps_per_period_, horizon_steps=self.horizon_steps_, date_col=self.date_col_, - categorical_levels=categorical_ranges(self.adaptive_model_.effects_list_) + categorical_levels=categorical_ranges(self.adaptive_model_.effects_list_), + fixed_ranges=_resolve_fixed_ranges(real_data, fixed_ranges(self.adaptive_model_.effects_list_)) ) # Data-dependent state (spline knots) comes from the training windows, never from a memory probe. initialize_effects( @@ -549,6 +624,58 @@ def predict_online(self, data: pd.DataFrame) -> pd.DataFrame: self._free_simulation() return self.predictions_ +#: + def predict_quantiles(self, data: pd.DataFrame, taus: Sequence[float] = (0.05, 0.5, 0.95)) -> pd.DataFrame: + r""" + Online quantile forecasts of a distributional base, the probabilistic counterpart of ``predict_online``. + + The adaptive correction moves the location; the scale of the base is kept and rescaled online. At each update of a + group, ``s_t`` is the root mean square of the standardized residuals ``(y - adapted location) / sigma_base`` of the + last ``training_window_periods`` rows that are at least ``horizon_steps`` old, and the quantile of level ``tau`` is + ``g^-1(adapted location + s_t * sigma_base * z_tau)`` with ``z_tau`` from the tail family of the base and ``g`` its + target transform. A row before the first window, or without a usable residual, keeps the base scale (``s_t = 1``). + Only past residuals enter ``s_t``: changing a later ``y`` leaves earlier quantiles unchanged. + + Like ``predict_online``, ``data`` must hold the target (a backtest over history); the state of the model is + refitted on it. + + Args: + data: The rows to forecast, with the target. + taus: Quantile levels in (0, 1). + + Returns: + A DataFrame with one ``q`` column per level, on the rows of ``predict_online(data)``. + + Raises: + ValueError: If the base model is not a distributional StaticTAM, or a level is outside (0, 1). + """ + base = self.distributional_base_ + if base is None: + raise ValueError("predict_quantiles needs a distributional StaticTAM as base_model (a dict formula, such as " + "{'mu': 'y ~ s(x)', 'sigma': '~ s(x)'}); a plain base has no scale to rescale.") + levels = [float(tau) for tau in taus] + if any(not 0.0 < tau < 1.0 for tau in levels): + raise ValueError(f"taus must lie in (0, 1); got {list(taus)!r}") + + from .statistics.estimation import _distributional + + frame = self.predict_online(data) + location = frame[f"AdaptedEstimated{_MODEL_TARGET}"].to_numpy(dtype=float) + observed = frame[_MODEL_TARGET].to_numpy(dtype=float) + mu_base, sigma_base = base._mu_sigma(frame) + standardized = (observed - location) / sigma_base + group = frame[self.group_col_].to_numpy() if self.group_col_ in frame.columns else np.zeros(len(frame)) + order_key = frame[self.date_col_] if self.date_col_ in frame.columns else pd.Series(np.arange(len(frame)), index=frame.index) + scale = _online_scale(standardized, group, order_key.to_numpy(), self.training_window_periods_, + self.update_interval_periods_, self.horizon_steps_) + columns = { + f"q{tau}": base._from_model_scale(location + scale * sigma_base * _distributional._standardized_ppf(base, tau), + f"the quantile at level {tau} (predict_quantiles)") + for tau in levels + } + return pd.DataFrame(columns, index=frame.index) +#: + def fit(self, data: pd.DataFrame) -> 'AdaptiveTAM': r""" Fits the adaptive model by solving the regularized linear system @@ -588,6 +715,7 @@ def predict(self, df: pd.DataFrame) -> pd.DataFrame: raise RuntimeError("Call fit() first to train the final adaptive state.") # --- 1. Base Model Extraction --- + df = self._with_model_target(df) if self.base_model_ is not None: data_bm = self.base_model_.decompose_prediction(df) data_pred = self.base_model_.predict(df) diff --git a/src/tam/model/additive.py b/src/tam/model/additive.py index 3bbc585..a015937 100644 --- a/src/tam/model/additive.py +++ b/src/tam/model/additive.py @@ -28,6 +28,7 @@ from ._base import BaseTAM from .safety import SafetyTAM from ._data import ( + _resolve_fixed_ranges, _fit_normalization_params, _transform_data_stacked, _groups_in_data, @@ -48,12 +49,14 @@ create_effects_from_parsed_terms, initialize_effects, categorical_ranges, + fixed_ranges, categorical_features, extrapolating_features, build_phi_from_effects, build_penalty_from_effects ) from .statistics.estimation import build_strategy, _distributional, _mixture +from .statistics.estimation._count_families import predict_count_quantiles _TINY: float = 1e-12 @@ -104,6 +107,10 @@ def __init__( you must pass date_col explicitly to guarantee correct results on data that isn't already sorted by time. default_alpha_p: Default log10(lambda_p) regularization strength. + dist_kwargs: Options of the distributional and mixture modes: ``tail_family``, ``kurtosis_threshold``, + ``scale_shrinkage``, ``location_alpha_p``, ``scale_alpha_p``, and the target scale: ``target_transform`` + (``"log"`` by default, ``"none"``, ``"asinh"`` for a skewed target with negative values, ``"logit"`` for a + bounded one, with ``bounds=(low, high)`` or estimated) or the older ``log_target`` (True = ``"log"``). _internal_effects_list: (Internal) Used for restoring state during grid search. _internal_features_config: (Internal) Used for restoring state during grid search. """ @@ -118,7 +125,9 @@ def __init__( dk = dist_kwargs or {} # Extracted universally since both distributional and mixture modes use the target scale transformation. - self._log_target_ = bool(dk.get("log_target", True)) + self._target_transform_, self._log_target_, self._transform_bounds_ = _distributional.resolve_target_transform( + dk, mixture=mixture_components is not None + ) # Conformal state (set by calibrate_conformal; used by predict_intervals). self._safety_ = None self._conformal_studentized_ = False @@ -155,6 +164,7 @@ def __init__( # scalar-loss reweighting strategy. None (l2/gaussian/normal, or mixture mode) keeps the default # single unweighted solve; any other loss builds an IRLS strategy from loss_kwargs. self.loss_ = loss + self.dispersion_ = None if self._mixture_components_ is not None: self._reweighting_strategy_ = None else: @@ -242,6 +252,9 @@ def predict(self, data: pd.DataFrame) -> pd.DataFrame: def _to_model_scale(self, y: np.ndarray) -> np.ndarray: return _distributional._to_model_scale(self, y) + def _from_model_scale(self, z: np.ndarray, quantity: str = "the back-transformed values") -> np.ndarray: + return _distributional.from_model_scale(self, z, quantity) + def _mu_sigma(self, data: pd.DataFrame) -> Tuple[np.ndarray, np.ndarray]: return _distributional.mu_sigma(self, data) @@ -252,6 +265,12 @@ def predict_quantile(self, data: pd.DataFrame, tau: Union[float, Sequence[float] return _distributional.predict_quantile(self, data, tau) def predict_quantiles(self, data: pd.DataFrame, taus: Sequence[float] = (0.05, 0.5, 0.95)) -> pd.DataFrame: + """Quantiles on the response scale: of the location-scale law (distributional mode), or of the fitted law of a + ``poisson``, ``negative_binomial`` or ``tweedie`` loss (discrete quantiles for the counts, the compound Poisson-gamma + law for Tweedie; non-negative by construction).""" + strategy = getattr(self, "_reweighting_strategy_", None) + if getattr(self, "_mode_", "plain") != "distributional" and hasattr(strategy, "quantiles"): + return predict_count_quantiles(self, data, taus) return _distributional.predict_quantiles(self, data, taus) def cdf(self, data: pd.DataFrame) -> np.ndarray: @@ -489,7 +508,8 @@ def _prepare_data( data=data, features=self.features_config_["features"], group_col=self.group_col_, - categorical_levels=categorical_ranges(self.effects_list_) + categorical_levels=categorical_ranges(self.effects_list_), + fixed_ranges=_resolve_fixed_ranges(data, fixed_ranges(self.effects_list_)) ) if target_col is not None: @@ -959,8 +979,8 @@ def show_syntax(): { "Token": "f(x)", "Effect": "Fourier", - "Syntax Example": "f(doy, m=6, s=1, cyclic=True)", - "Specific Params": "m (harmonics), s (smoothness), cyclic (bool)" + "Syntax Example": "f(hour, m=6, s=1, cyclic=True, period='auto')", + "Specific Params": "m (harmonics), s (smoothness), cyclic (bool), period (None: training range; 'auto': first to last value + one step; (low, high): given)" }, { "Token": "c(x)", diff --git a/src/tam/model/kalman.py b/src/tam/model/kalman.py index 7347fe2..f2c0d3f 100644 --- a/src/tam/model/kalman.py +++ b/src/tam/model/kalman.py @@ -38,7 +38,8 @@ import numpy as np import itertools import warnings -from typing import Tuple, Dict, Any, Optional, Union +from typing import Tuple, Dict, Any, Optional, Sequence, Union +from scipy import stats from tam.common.utils import ( TORCH_DEVICE, _balance_groups, @@ -60,12 +61,19 @@ def _kalman_block_loop_optimized( process_noise_diag: torch.Tensor, eps: float, offset_boost: float, - offset_q_boost: float -) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]: + offset_q_boost: float, + with_variance: bool, + horizon_steps: int +) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]: r""" Compiled TorchScript loop for the Block Kalman Filter. Executes the Predict and Update phases across all groups simultaneously while keeping the temporal recursion entirely within the compiled C++ backend. + + With ``with_variance`` the fourth output holds, for every forecast row, the predictive variance ``x' P x + R`` of the + standardised target: ``P`` is the covariance of the state that forecasts the row (a priori at its block; for a delayed + state of ``horizon_steps > 1``, the covariance of the state ``horizon_steps - 1`` steps before, grown by that many steps + of process noise). Without it the fourth output is empty and the filter is unchanged. """ G, N, d = phi_matrix.shape dtype = phi_matrix.dtype @@ -85,6 +93,12 @@ def _kalman_block_loop_optimized( Q_matrix[:, 0, 0] = Q_matrix[:, 0, 0] * offset_q_boost predictions = torch.zeros((G, N, 1), device=device, dtype=dtype) + variance_history = torch.zeros((G, N if with_variance else 0, 1), device=device, dtype=dtype) + if with_variance and horizon_steps > 1: + # The first rows are forecast by the initial state (zero), whose covariance is the initial one. + for j in range(min(horizon_steps - 1, N)): + x_j = phi_matrix[:, j : j + 1, :] + variance_history[:, j : j + 1, :] = torch.bmm(torch.bmm(x_j, P_t), x_j.transpose(-2, -1)) + observation_noise_var num_blocks = (N + B - 1) // B state_history_gpu = torch.zeros((num_blocks, G, d, 1), device=device, dtype=dtype) @@ -103,6 +117,18 @@ def _kalman_block_loop_optimized( predictions[:, t : t + curr_B, :] = y_hat_B state_history_gpu[i] = theta_t + if with_variance: + if horizon_steps > 1: + forecast_row = t + horizon_steps - 1 + if forecast_row < N: + x_h = phi_matrix[:, forecast_row : forecast_row + 1, :] + P_h = P_t + Q_matrix * float(horizon_steps - 1) + variance_history[:, forecast_row : forecast_row + 1, :] = ( + torch.bmm(torch.bmm(x_h, P_h), x_h.transpose(-2, -1)) + observation_noise_var + ) + else: + variance_history[:, t : t + curr_B, :] = (torch.bmm(X_B, P_t) * X_B).sum(dim=-1, keepdim=True) + observation_noise_var + # --- UPDATE PHASE (A-Posteriori) --- innovations = Y_B - y_hat_B @@ -131,7 +157,7 @@ def _kalman_block_loop_optimized( # state_history_gpu holds the state each block was predicted with (before its update); theta_t is the state after # the last update: the one that forecasts the next, unseen steps. - return predictions, state_history_gpu, theta_t + return predictions, state_history_gpu, theta_t, variance_history #: class KalmanTAM: @@ -485,11 +511,15 @@ def _run_filter( prepared_data: Dict[str, Any], P_init_diag: float, observation_noise_var: float, - process_noise_var: Union[float, Dict[str, float]] + process_noise_var: Union[float, Dict[str, float]], + with_variance: bool = False ) -> torch.Tensor: """ Executes the core Kalman filter loop on normalized data and denormalizes the results before returning. + + With ``with_variance`` the predictive variance of every forecast row (standardised target space) is also kept in + ``_predictive_variance_``, shape ``(groups, rows, 1)``; the returned predictions are the same. """ phi_matrix = prepared_data["phi_matrix"] y_stacked = prepared_data["y_stacked"] @@ -498,7 +528,7 @@ def _run_filter( actual_block_size = 1 if self.horizon_steps_ > 1 else self.block_size_ noise_diag, offset_q_boost = self._noise_diagonal(process_noise_var, phi_matrix.shape[-1]) - predictions_norm, state_history_gpu, final_state_gpu = _kalman_block_loop_optimized( + predictions_norm, state_history_gpu, final_state_gpu, variance_history = _kalman_block_loop_optimized( phi_matrix=phi_matrix, y_stacked=y_stacked, base_pred_stacked=base_pred_stacked, @@ -508,8 +538,11 @@ def _run_filter( process_noise_diag=noise_diag, eps=self.eps, offset_boost=float(self.offset_boost), - offset_q_boost=offset_q_boost + offset_q_boost=offset_q_boost, + with_variance=bool(with_variance), + horizon_steps=int(self.horizon_steps_) ) + self._predictive_variance_ = variance_history if with_variance else None self.states_history_ = state_history_gpu.squeeze(-1).cpu() self.final_state_ = final_state_gpu.squeeze(-1).cpu() # (G, d): the state after the last observation @@ -632,12 +665,17 @@ def predict_online(self, data: pd.DataFrame, **kwargs) -> pd.DataFrame: Main interface to adapt a base model to current drift. Reassembles tensors into the user's original DataFrame shape. """ + return self._online(data) + + def _online(self, data: pd.DataFrame, with_variance: bool = False) -> pd.DataFrame: + """The online simulation of ``predict_online``; with ``with_variance`` the frame also holds the ``__sd__`` column.""" prepared_data = self.prepare_data(data) preds_denorm = self._run_filter( prepared_data, self.P_init_diag_, self.observation_noise_var_, - self.process_noise_var_ + self.process_noise_var_, + with_variance=with_variance ) df_final = _reassemble_predictions( @@ -651,7 +689,15 @@ def predict_online(self, data: pd.DataFrame, **kwargs) -> pd.DataFrame: res_col = f"Estimated{self.target_col_}" df_final = df_final.rename(columns={res_col: f"KalmanAdapted_{self.target_col_}"}) - + if with_variance: + # Standard deviation of the predictive distribution in the target scale (the variance is in standardised units). + sd_stacked = torch.sqrt(self._predictive_variance_.clamp_min(0.0)) * prepared_data["y_scale"] + sd_frame = _reassemble_predictions( + prepared_data["df_original"], sd_stacked.squeeze(-1).cpu(), self.group_col_, + prepared_data["unique_groups"], "PredictiveSd", date_col=self.date_col_ + ) + df_final["__sd__"] = sd_frame["EstimatedPredictiveSd"] + df_final_masked = df_final[prepared_data["mask"]] # Save frozen states and scales for out-of-sample inference @@ -663,6 +709,39 @@ def predict_online(self, data: pd.DataFrame, **kwargs) -> pd.DataFrame: return _cleanup_dummies(df_final_masked, self.group_col_, self.date_col_) +#: + def predict_quantiles(self, data: pd.DataFrame, taus: Sequence[float] = (0.05, 0.5, 0.95)) -> pd.DataFrame: + r""" + Gaussian quantiles of the online forecast, from the predictive variance of the filter. + + The forecast of a row is Gaussian with the mean of ``predict_online`` and the variance ``x' P x + R`` (``P`` the covariance + of the state that forecasts the row, ``R`` the observation noise; for ``horizon_steps > 1`` the covariance of the delayed + state grown by the process noise of the delay), converted back to the target scale. Like ``predict_online``, ``data`` must + hold the target and the state of the model is refitted on it. The variance is a function of the design and the noise + parameters only, so a later target never moves an earlier quantile. + + The intervals are as calibrated as ``observation_noise_var`` and ``process_noise_var`` are (``tune_hyperparameters`` + sets them). A row without a target does not move the state but still shrinks the covariance, as in the filter. + + Args: + data: The rows to forecast, with the target. + taus: Quantile levels in (0, 1). + + Returns: + A DataFrame with one ``q`` column per level, on the rows of ``predict_online(data)``. + + Raises: + ValueError: If a level is outside (0, 1). + """ + levels = [float(tau) for tau in taus] + if any(not 0.0 < tau < 1.0 for tau in levels): + raise ValueError(f"taus must lie in (0, 1); got {list(taus)!r}") + frame = self._online(data, with_variance=True) + mean = frame[f"KalmanAdapted_{self.target_col_}"].to_numpy(dtype=float) + sd = frame["__sd__"].to_numpy(dtype=float) + return pd.DataFrame({f"q{tau}": mean + sd * stats.norm.ppf(tau) for tau in levels}, index=frame.index) +#: + def compact(self) -> 'KalmanTAM': r""" Keeps the last row of ``states_history_`` (the state after the last block) and drops the rest, which grows with the number of rows. diff --git a/src/tam/model/opera.py b/src/tam/model/opera.py index c86e54e..3fa8010 100644 --- a/src/tam/model/opera.py +++ b/src/tam/model/opera.py @@ -35,7 +35,8 @@ def _mlpol_loop_optimized_3d( experts_tensor: torch.Tensor, y_true: torch.Tensor, - loss_type: str + loss_type: str, + tau: float = 0.5 ) -> Tuple[torch.Tensor, torch.Tensor]: r""" Compiled TorchScript loop for MLpol (Polynomial Minimax Strategy). @@ -45,7 +46,8 @@ def _mlpol_loop_optimized_3d( Args: experts_tensor: Predictions from experts, shape (B, T, K). y_true: Ground truth targets, shape (B, T, 1). - loss_type: The loss function to use ('square' or 'absolute'). + loss_type: The loss function to use ('square', 'absolute' or 'pinball'). + tau: Quantile level of the pinball loss (ignored by the other losses). Returns: Tuple containing the mixed predictions (B, T) and weights history (B, T, K). @@ -91,6 +93,10 @@ def _mlpol_loop_optimized_3d( if loss_type == 'square': r = 2.0 * (y_hat_scaled.unsqueeze(1) - yt_scaled.unsqueeze(1)) * (y_hat_scaled.unsqueeze(1) - xt_scaled) + elif loss_type == 'pinball': + # Subgradient of the pinball loss in the forecast: 1{y_hat > y} - tau. + above = (y_hat_scaled > yt_scaled).to(dtype) + r = (above - tau).unsqueeze(1) * (y_hat_scaled.unsqueeze(1) - xt_scaled) else: r = torch.sign(y_hat_scaled.unsqueeze(1) - yt_scaled.unsqueeze(1)) * (y_hat_scaled.unsqueeze(1) - xt_scaled) @@ -126,7 +132,8 @@ def _ewa_loop_optimized_3d( experts_tensor: torch.Tensor, y_true: torch.Tensor, loss_type: str, - eta: float + eta: float, + tau: float = 0.5 ) -> Tuple[torch.Tensor, torch.Tensor]: r""" Compiled TorchScript loop for EWA (Exponentially Weighted Aggregation). @@ -155,6 +162,9 @@ def _ewa_loop_optimized_3d( # Compute expert losses if loss_type == 'square': losses = (xt - yt.unsqueeze(1)) ** 2 + elif loss_type == 'pinball': + diff = yt.unsqueeze(1) - xt + losses = torch.maximum(tau * diff, (tau - 1.0) * diff) else: losses = torch.abs(xt - yt.unsqueeze(1)) @@ -189,7 +199,8 @@ class OperaTAM: formula (str): Aggregation topology (e.g., 'Target ~ l(Expert1) + l(Expert2)'). algorithm (str): The aggregation strategy to use ('EWA' or 'MLPOL'). eta (float): Learning rate parameter (primarily for EWA). - loss_type (str): The loss function to evaluate experts ('square' or 'absolute'). + loss_type (str): The loss function to evaluate experts ('square', 'absolute' or 'pinball'). + tau (float): Quantile level of the pinball loss. horizon_steps (int): Number of steps to shift the learned weights forward to enforce causal boundaries in multi-step forecasting. group_col (str, optional): Column name for independent group processing. date_col (str): Column name for time indexing to ensure chronological evaluation. @@ -208,7 +219,8 @@ def __init__( loss_type: str = 'square', horizon_steps: int = 1, group_col: Optional[str] = None, - date_col: Optional[str] = None + date_col: Optional[str] = None, + tau: float = 0.5 ): r""" Initializes the Native GPU Aggregator. @@ -219,10 +231,11 @@ def __init__( expert_cols (list, optional): List of expert column names (used if formula is not provided). algorithm (str): 'EWA' or 'MLPOL'. eta (float): Learning rate parameter (used by EWA). - loss_type (str): Loss formulation ('square' or 'absolute'). + loss_type (str): Loss formulation ('square', 'absolute', or 'pinball' to aggregate quantile forecasts of level tau). horizon_steps (int): Number of steps to shift the learned weights forward to enforce causal boundaries in multi-step forecasting. Default is 1. group_col (str, optional): If provided, runs aggregation independently per group. date_col (str, optional): Column name utilized for temporal sorting and padding. + tau (float): Quantile level in (0, 1) of the pinball loss (used only with loss_type='pinball'). """ if formula is None: @@ -240,12 +253,15 @@ def __init__( self.date_col = date_col or "__dummy_date__" self.loss_type = loss_type self.horizon_steps = horizon_steps + self.tau = float(tau) if self.algorithm not in ['EWA', 'MLPOL']: raise ValueError("Algorithm must be 'EWA' or 'MLPOL'.") - if self.loss_type not in ['square', 'absolute']: - raise ValueError("Unsupported loss_type. Use 'square' or 'absolute'.") + if self.loss_type not in ['square', 'absolute', 'pinball']: + raise ValueError("Unsupported loss_type. Use 'square', 'absolute' or 'pinball'.") + if self.loss_type == 'pinball' and not 0.0 < self.tau < 1.0: + raise ValueError(f"The pinball loss needs a quantile level tau in (0, 1); got {tau}.") if '~' not in self.formula: raise ValueError("Formula must contain '~' separating target and experts.") @@ -308,9 +324,9 @@ def predict_online(self, df: pd.DataFrame) -> pd.DataFrame: # Computation Phase: Execute the entire historical simulation in one GPU pass if self.algorithm == 'MLPOL': - preds_3d, weights_3d = _mlpol_loop_optimized_3d(X_tensor_3d, Y_tensor_3d, self.loss_type) + preds_3d, weights_3d = _mlpol_loop_optimized_3d(X_tensor_3d, Y_tensor_3d, self.loss_type, self.tau) else: - preds_3d, weights_3d = _ewa_loop_optimized_3d(X_tensor_3d, Y_tensor_3d, self.loss_type, self.eta) + preds_3d, weights_3d = _ewa_loop_optimized_3d(X_tensor_3d, Y_tensor_3d, self.loss_type, self.eta, self.tau) # Reassembly Phase preds_np = preds_3d.cpu().numpy() @@ -543,4 +559,47 @@ def plot_weights( plt.tight_layout(rect=[0, 0, 0.85, 0.95]) - plt.show() \ No newline at end of file + plt.show() + + +#: +def aggregate_quantiles( + df: pd.DataFrame, + target_col: str, + expert_cols_by_level: Dict[float, list], + algorithm: str = 'MLPOL', + eta: float = 1.0, + horizon_steps: int = 1, + group_col: Optional[str] = None, + date_col: Optional[str] = None, +) -> pd.DataFrame: + r""" + Online aggregation of quantile forecasts: one ``OperaTAM`` with the pinball loss per level, levels sorted on each row. + + Args: + df: Data holding the target and, for each level, the experts' forecasts of that quantile. + target_col: Target column. + expert_cols_by_level: ``{tau: [expert columns forecasting the tau-quantile]}``. + algorithm: 'MLPOL' or 'EWA'. + eta: Learning rate (EWA). + horizon_steps: Delay of the weights, as in ``OperaTAM``. + group_col: Optional group column (independent aggregation per group). + date_col: Optional date column (chronological order). + + Returns: + A DataFrame with one ``q`` column per level, in the row order of ``OperaTAM.predict_online`` and the level + order of ``expert_cols_by_level``; the aggregated levels are sorted on each row so they never cross. + """ + levels = list(expert_cols_by_level) + columns = [] + index = None + for tau in levels: + out = OperaTAM(target_col=target_col, expert_cols=list(expert_cols_by_level[tau]), algorithm=algorithm, eta=eta, + loss_type='pinball', tau=tau, horizon_steps=horizon_steps, group_col=group_col, + date_col=date_col).predict_online(df) + index = out.index if index is None else index + columns.append(out.loc[index, 'prediction_opera'].to_numpy(dtype=float)) + ordered = np.argsort(levels) + table = np.sort(np.column_stack(columns)[:, ordered], axis=1)[:, np.argsort(ordered)] + return pd.DataFrame(table, index=index, columns=[f"q{tau}" for tau in levels]) +#: diff --git a/src/tam/model/safety.py b/src/tam/model/safety.py index 7bed586..64761fd 100644 --- a/src/tam/model/safety.py +++ b/src/tam/model/safety.py @@ -33,6 +33,7 @@ def __init__(self, alpha: float = 0.1): """alpha is the target error rate (e.g. 0.1 for 90% coverage).""" self.alpha_target = alpha self.residuals_calib_ = None + self.signed_scores_ = None def calibrate_scores(self, scores: np.ndarray) -> "SafetyTAM": """Calibrate directly on a supplied array of nonconformity scores (the generic entry point).""" @@ -45,18 +46,35 @@ def calibrate( y_pred: np.ndarray, scale: Optional[np.ndarray] = None, scores: Optional[np.ndarray] = None, + verbose: bool = True, ) -> "SafetyTAM": """Calibrate on a hold-out set. The nonconformity score defaults to the absolute residual `|y - y_pred|`; passing scale normalizes it - to `|y - y_pred|` / scale (studentized conformal), and passing scores uses them directly. + to `|y - y_pred|` / scale (studentized conformal), and passing scores uses them directly. The signed + residuals `(y - y_pred)` / scale are kept too, for the per-tail quantiles of `predict_quantiles(..., signed=True)`. + + Args: + y_true: Observed values of the hold-out set. + y_pred: Point predictions for the same rows. + scale: Optional per-row scale (studentized scores). + scores: Optional nonconformity scores used as given (then no signed residuals are kept). + verbose: Print the number of calibration samples (the default, as before); False is silent. + + Returns: + self """ + self.signed_scores_ = None if scores is not None: self.residuals_calib_ = np.asarray(scores, dtype=float) else: - residual = np.abs(np.asarray(y_true, dtype=float) - np.asarray(y_pred, dtype=float)) - self.residuals_calib_ = residual / np.asarray(scale, dtype=float) if scale is not None else residual - print(f"Safety calibrated on {len(self.residuals_calib_)} samples.") + signed = np.asarray(y_true, dtype=float) - np.asarray(y_pred, dtype=float) + if scale is not None: + signed = signed / np.asarray(scale, dtype=float) + self.signed_scores_ = signed + self.residuals_calib_ = np.abs(signed) + if verbose: + print(f"Safety calibrated on {len(self.residuals_calib_)} samples.") return self def conformal_quantile(self, alpha: Optional[float] = None) -> float: @@ -117,3 +135,59 @@ def predict_intervals( result['Actual'] = np.asarray(y_true_online, dtype=float) result['Covered'] = (result['Actual'] >= result['Lower']) & (result['Actual'] <= result['Upper']) return result + + def predict_quantiles( + self, + y_pred: np.ndarray, + levels, + scale: Optional[np.ndarray] = None, + signed: bool = False, + ) -> pd.DataFrame: + """Conformal quantiles at several levels in one call. + + Symmetric mode (default): the level tau < 0.5 is `y_pred - q * scale` and 1 - tau is `y_pred + q * scale`, with + q = conformal_quantile(2 * tau) of the absolute scores (the bands of `predict_intervals`); 0.5 is `y_pred`. + Signed mode: each level is `y_pred + s_(k) * scale`, s_(k) the k-th smallest signed calibration score, with + k = ceil((n + 1) * tau) above the median and k = floor((n + 1) * tau) below it (one-sided split conformal), + so skewed errors get a short and a long tail. The levels are sorted on each row (no crossing). + + Args: + y_pred: Point predictions. + levels: Quantile levels in (0, 1). + scale: Optional per-row scale (the one used at calibration for studentized scores). + signed: Use the signed scores kept by `calibrate` (not available after `calibrate_scores`). + + Returns: + A DataFrame with one column per level, named `q`. + + Raises: + RuntimeError: If the engine is not calibrated. + ValueError: If `signed=True` but no signed scores were kept, or a level is outside (0, 1). + """ + if self.residuals_calib_ is None: + raise RuntimeError("You must call .calibrate() or .calibrate_scores() before predicting quantiles.") + levels = [float(t) for t in levels] + if any(not 0.0 < t < 1.0 for t in levels): + raise ValueError(f"Quantile levels must lie in (0, 1); got {levels}.") + if signed and self.signed_scores_ is None: + raise ValueError("signed=True needs the signed residuals kept by .calibrate(y_true, y_pred); " + ".calibrate_scores() keeps only absolute scores.") + y_pred = np.asarray(y_pred, dtype=float) + scale_vector = np.ones_like(y_pred) if scale is None else np.asarray(scale, dtype=float) + + columns = [] + if signed: + ordered = np.sort(self.signed_scores_) + n = len(ordered) + for tau in levels: + k = int(np.ceil((n + 1) * tau)) if tau >= 0.5 else int(np.floor((n + 1) * tau)) + columns.append(y_pred + ordered[min(max(k, 1), n) - 1] * scale_vector) + else: + for tau in levels: + if tau == 0.5: + columns.append(y_pred.copy()) + continue + radius = self.conformal_quantile(round(2.0 * min(tau, 1.0 - tau), 12)) + columns.append(y_pred + np.sign(tau - 0.5) * radius * scale_vector) + table = np.sort(np.column_stack(columns), axis=1) + return pd.DataFrame(table, columns=[f"q{t}" for t in levels]) diff --git a/src/tam/model/spectrum/__init__.py b/src/tam/model/spectrum/__init__.py index e78970a..15ddfe9 100644 --- a/src/tam/model/spectrum/__init__.py +++ b/src/tam/model/spectrum/__init__.py @@ -46,6 +46,7 @@ create_effects_from_parsed_terms, initialize_effects, categorical_ranges, + fixed_ranges, categorical_range, categorical_features, extrapolating_features, @@ -77,6 +78,7 @@ "create_effects_from_parsed_terms", "initialize_effects", "categorical_ranges", + "fixed_ranges", "categorical_range", "categorical_features", "extrapolating_features", diff --git a/src/tam/model/spectrum/_factory.py b/src/tam/model/spectrum/_factory.py index 1f8f17a..acf7c73 100644 --- a/src/tam/model/spectrum/_factory.py +++ b/src/tam/model/spectrum/_factory.py @@ -103,7 +103,8 @@ def create_effects_from_parsed_terms( cyclic_raw = params_resolved.get('cyclic', 'False') is_cyclic = str(cyclic_raw).strip().lower() in ['true', '1', 't', 'y', 'yes'] extrap_val = params_resolved.get('extrapolate', 'continue') - effects_list.append(FourierEffect(feature_name, m, s, lambda_p, is_cyclic, extrap_val)) + period = _parse_period(params_resolved.get('period'), feature_name) + effects_list.append(FourierEffect(feature_name, m, s, lambda_p, is_cyclic, extrap_val, period=period)) #: #: @@ -355,6 +356,64 @@ def visit(effect) -> None: return {name: int(next(iter(n))) for name, n in levels.items() if name not in other and len(n) == 1} +def _parse_period(raw, feature_name: str): + """The ``period`` argument of ``f()``: None when absent, ``"auto"``, or ``(low, high)`` (a tuple or the formula text ``"(0, 1)"``).""" + if raw is None: + return None + if isinstance(raw, str) and raw.strip().strip("'\"").lower() == "auto": + return "auto" + if isinstance(raw, (tuple, list)): + values = list(raw) + else: + values = [v for v in str(raw).strip().strip("()[]").split(",") if v.strip()] + try: + low, high = (float(v) for v in values) + except (TypeError, ValueError): + raise ValueError(f"Invalid 'period' for f({feature_name}): {raw!r}; expected 'auto' or two numbers, e.g. period=(0, 1).") + return low, high + + +def fixed_ranges(effects_list: List[BaseEffect]) -> Dict[str, tuple]: + """ + ``{feature: (low, high) or "auto"}`` for the features read by a Fourier effect with a ``period`` (``f(x, cyclic=True, period=(0, 1))``). + + Every feature is normalised to [-1, 1] with the minimum and maximum of its training rows, and a cyclic term takes [-1, 1] as + one period: a history covering half a cycle (January to June of a time of year in [0, 1]) would be read as a whole cycle and + the forecast of July would start the cycle again. A declared period is the range the feature is normalised on instead, so the + history is read as the part of the cycle it covers. ``"auto"`` is resolved on the training rows by ``_resolve_fixed_ranges`` + (first value, last value + one step). Without ``period`` nothing changes. The declared range also applies to + the other non-categorical effects reading the feature (an affine rescaling for them); a feature also read by a categorical + effect keeps the categorical rule. + + Args: + effects_list: The effects of a model, tensor products included. + + Returns: + The features to normalise on a declared range, with that range. + + Raises: + ValueError: If two Fourier effects declare different periods for the same feature. + """ + periods: Dict[str, tuple] = {} + categorical: set = set() + + def visit(effect) -> None: + if isinstance(effect, TensorProductEffect): + for sub in effect.effects: + visit(sub) + return + if isinstance(effect, CategoricalEffect): + categorical.add(effect.feature_name) + elif isinstance(effect, FourierEffect) and getattr(effect, 'period', None) is not None: + known = periods.setdefault(effect.feature_name, effect.period) + if known != effect.period: + raise ValueError(f"Two periods declared for '{effect.feature_name}': {known} and {effect.period}.") + + for effect in effects_list: + visit(effect) + return {name: rng for name, rng in sorted(periods.items()) if name not in categorical} + + def categorical_range(n_cat: int, low: float, high: float, integer_codes: bool = True) -> tuple: """ The ``(min, max)`` a categorical code is normalised with, from the smallest and largest code of the training rows. diff --git a/src/tam/model/spectrum/_fourier.py b/src/tam/model/spectrum/_fourier.py index b46433e..0c9e0e4 100644 --- a/src/tam/model/spectrum/_fourier.py +++ b/src/tam/model/spectrum/_fourier.py @@ -9,6 +9,8 @@ Projects data onto a truncated Fourier basis. """ +from typing import Optional, Tuple, Union + import torch import numpy as np from tam.common.utils import TORCH_DEVICE @@ -23,15 +25,29 @@ class FourierEffect(BaseEffect): If cyclic=True, it enforces a strict periodic boundary. """ - def __init__(self, feature_name: str, m: int, s: int, lambda_p: float, cyclic: bool, extrapolate: str): + def __init__(self, feature_name: str, m: int, s: int, lambda_p: float, cyclic: bool, extrapolate: str, + period: Optional[Union[str, Tuple[float, float]]] = None): + """ + Args: + period: The range of the input that makes one period. ``None`` (default): the training range (minimum, maximum), + as before. ``"auto"``: (first value, last value + one step), the step being the smallest gap between the + distinct training values, so hours 0..23 give (0, 24) and months 1..12 give (1, 13). ``(low, high)``: given, + for a history that covers only part of the cycle (``(0, 1)`` for a time of year). + """ if m <= 0 or s < 0: raise ValueError(f"Invalid params for FourierEffect: m={m}, s={s}") - + if isinstance(period, str): + if period != "auto": + raise ValueError(f"FourierEffect period must be None, 'auto' or (low, high); got {period!r}") + elif period is not None and not float(period[0]) < float(period[1]): + raise ValueError(f"FourierEffect period must be (low, high) with low < high; got {period}") + super().__init__(feature_name, "fourier", lambda_p, extrapolate) - + self.m = int(m) self.s = int(s) self.cyclic = cyclic + self.period = period if period is None or isinstance(period, str) else (float(period[0]), float(period[1])) def get_n_coeffs(self) -> int: return 2 * self.m diff --git a/src/tam/model/statistics/estimation/_count_families.py b/src/tam/model/statistics/estimation/_count_families.py new file mode 100644 index 0000000..57d5df5 --- /dev/null +++ b/src/tam/model/statistics/estimation/_count_families.py @@ -0,0 +1,216 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +""" +Negative binomial and Tweedie families (log link) and their quantiles. + +Both are GLM families like ``poisson_family``: the P-WLS atom is reweighted with a variance function, here with a dispersion that +is estimated after each IRLS fit (negative binomial: the shape of the variance ``mu + alpha mu^2``; Tweedie: the scale ``phi`` of +``phi mu^p``). A predicted mean gives a whole law, so ``quantiles`` returns discrete quantiles for the counts (Poisson, negative +binomial) and the quantiles of the compound Poisson-gamma law for the Tweedie family, which carries a point mass at zero. +""" +from __future__ import annotations + +from typing import Optional, Sequence + +import numpy as np +import pandas as pd +import torch +from scipy import optimize, special, stats + +from ._glm import GLMFamily, LogLink, poisson_family + +_TINY: float = 1e-12 +# The dispersion is re-estimated until it moves by less than this relative amount. +_DISPERSION_TOL: float = 1e-3 +# Poisson mass beyond which the number of jumps of a compound Poisson-gamma is cut. +_JUMP_TAIL: float = 1e-10 +_MAX_JUMPS: int = 2000 +_BISECTION_STEPS: int = 45 +_CHUNK: int = 4000 + + +def _check_levels(taus: Sequence[float]) -> list: + levels = [float(tau) for tau in taus] + if any(not 0.0 < tau < 1.0 for tau in levels): + raise ValueError(f"taus must lie in (0, 1); got {list(taus)!r}") + return levels + + +def _check_non_negative(y: torch.Tensor, family: str) -> None: + n_bad = int(torch.count_nonzero(y < 0)) + if n_bad: + raise ValueError(f"The {family} target must be >= 0; {n_bad} negative values found.") + + +#: +class NegativeBinomialFamily(GLMFamily): + """Negative binomial (NB2) with log link: variance ``mu + alpha mu^2``, ``alpha = 1 / theta`` estimated by maximum likelihood.""" + + def __init__(self, dispersion: Optional[float] = None): + if dispersion is not None and not float(dispersion) > 0.0: + raise ValueError(f"dispersion must be > 0; got {dispersion!r}") + self.fixed_dispersion = None if dispersion is None else float(dispersion) + self.dispersion = self.fixed_dispersion if dispersion is not None else 0.1 + super().__init__(LogLink(), self._variance, self._unit_deviance, "negative_binomial") + + def _variance(self, mu: torch.Tensor) -> torch.Tensor: + return mu + self.dispersion * mu * mu + + def _unit_deviance(self, y: torch.Tensor, mu: torch.Tensor) -> torch.Tensor: + alpha = self.dispersion + y_pos, mu_pos = y.clamp_min(_TINY), mu.clamp_min(_TINY) + return 2.0 * (y * torch.log(y_pos / mu_pos) - (y + 1.0 / alpha) * torch.log((1.0 + alpha * y) / (1.0 + alpha * mu_pos))) + + def validate_target(self, y: torch.Tensor) -> None: + _check_non_negative(y, "negative binomial") + + def estimate_dispersion(self, y: torch.Tensor, mu: torch.Tensor) -> bool: + """Maximum-likelihood ``alpha`` given the fitted means; returns whether it moved enough to refit the mean.""" + if self.fixed_dispersion is not None: + return False + observed = y.detach().cpu().numpy().ravel().astype(float) + mean = np.clip(mu.detach().cpu().numpy().ravel().astype(float), _TINY, None) + + def negative_log_likelihood(log_theta: float) -> float: + theta = np.exp(log_theta) + return -float(np.sum( + special.gammaln(observed + theta) - special.gammaln(theta) + + theta * np.log(theta / (theta + mean)) + observed * np.log(mean / (theta + mean)) + )) + + best = optimize.minimize_scalar(negative_log_likelihood, bounds=(-8.0, 12.0), method="bounded", options={"xatol": 1e-6}) + updated = float(np.exp(-best.x)) + moved = abs(updated - self.dispersion) / self.dispersion > _DISPERSION_TOL + self.dispersion = updated + return bool(moved) + + def quantiles(self, mu: np.ndarray, taus: Sequence[float]) -> np.ndarray: + """Discrete quantiles (smallest ``k`` with ``F(k) >= tau``), shape ``(len(mu), len(taus))``.""" + theta = 1.0 / self.dispersion + mean = np.clip(np.asarray(mu, dtype=float), _TINY, None) + return np.stack([stats.nbinom.ppf(tau, theta, theta / (theta + mean)) for tau in _check_levels(taus)], axis=1) +#: + + +#: +class TweedieFamily(GLMFamily): + """Tweedie with power ``1 < p < 2`` and log link: variance ``phi mu^p``, a point mass at zero and a positive continuous part. + + It is the compound Poisson-gamma law: ``N ~ Poisson(lambda)`` jumps, each ``Gamma(a, b)``, with + ``lambda = mu^(2-p) / (phi (2-p))``, ``a = (2-p) / (p-1)`` and ``b = phi (p-1) mu^(p-1)``. ``phi`` is estimated from the Pearson + residuals (it does not change the weights, so no refit). + """ + + def __init__(self, power: float = 1.5): + if not 1.0 < float(power) < 2.0: + raise ValueError(f"The Tweedie power must lie in (1, 2) (compound Poisson-gamma); got {power!r}") + self.power = float(power) + self.dispersion = 1.0 + super().__init__(LogLink(), lambda mu: mu.clamp_min(_TINY) ** self.power, self._unit_deviance, "tweedie") + + def _unit_deviance(self, y: torch.Tensor, mu: torch.Tensor) -> torch.Tensor: + p = self.power + y_pos, mu_pos = y.clamp_min(0.0), mu.clamp_min(_TINY) + return 2.0 * (y_pos ** (2.0 - p) / ((1.0 - p) * (2.0 - p)) - y * mu_pos ** (1.0 - p) / (1.0 - p) + mu_pos ** (2.0 - p) / (2.0 - p)) + + def validate_target(self, y: torch.Tensor) -> None: + _check_non_negative(y, "Tweedie") + + def estimate_dispersion(self, y: torch.Tensor, mu: torch.Tensor) -> bool: + """Pearson estimate of ``phi``; the mean is not refitted (a constant factor of the weights).""" + observed = y.detach().cpu().numpy().ravel().astype(float) + mean = np.clip(mu.detach().cpu().numpy().ravel().astype(float), _TINY, None) + self.dispersion = float(max(np.mean((observed - mean) ** 2 / mean ** self.power), _TINY)) + return False + + def _jump_parameters(self, mu: np.ndarray): + p, phi = self.power, self.dispersion + rate = mu ** (2.0 - p) / (phi * (2.0 - p)) + shape = (2.0 - p) / (p - 1.0) + scale = phi * (p - 1.0) * mu ** (p - 1.0) + return rate, shape, scale + + def probability_of_zero(self, mu: np.ndarray) -> np.ndarray: + """``P(y = 0) = exp(-lambda)``.""" + rate, _, _ = self._jump_parameters(np.clip(np.asarray(mu, dtype=float), _TINY, None)) + return np.exp(-rate) + + def cdf(self, y: np.ndarray, mu: np.ndarray) -> np.ndarray: + """``F(y) = P(0) + sum_n P(N = n) GammaCDF(y; n a, b)`` for ``y >= 0``, summed over the jumps that carry mass.""" + mu = np.clip(np.asarray(mu, dtype=float), _TINY, None) + y = np.asarray(y, dtype=float) + rate, shape, scale = self._jump_parameters(mu) + top = int(min(max(float(np.max(stats.poisson.ppf(1.0 - _JUMP_TAIL, rate))), 1.0) + 1.0, _MAX_JUMPS)) + jumps = np.arange(1, top + 1)[None, :] + weight = stats.poisson.pmf(jumps, rate[:, None]) + mass = special.gammainc(jumps * shape, np.maximum(y, 0.0)[:, None] / scale[:, None]) + return np.where(y < 0.0, 0.0, np.exp(-rate) + np.sum(weight * mass, axis=1)) + + def quantiles(self, mu: np.ndarray, taus: Sequence[float]) -> np.ndarray: + """Quantiles of the compound Poisson-gamma law by bisection on the CDF; ``0`` where ``tau <= P(y = 0)``.""" + mu = np.clip(np.asarray(mu, dtype=float), _TINY, None) + levels = _check_levels(taus) + zero = self.probability_of_zero(mu) + out = np.zeros((len(mu), len(levels))) + spread = np.sqrt(self.dispersion * mu ** self.power) + for column, tau in enumerate(levels): + rows = np.flatnonzero(tau > zero) + for start in range(0, len(rows), _CHUNK): + chunk = rows[start:start + _CHUNK] + low, high = np.zeros(len(chunk)), mu[chunk] + 12.0 * spread[chunk] + 1e-9 + while True: + short = self.cdf(high, mu[chunk]) < tau + if not short.any(): + break + high = np.where(short, high * 2.0, high) + for _ in range(_BISECTION_STEPS): + middle = 0.5 * (low + high) + below = self.cdf(middle, mu[chunk]) < tau + low, high = np.where(below, middle, low), np.where(below, high, middle) + out[chunk, column] = high + return out +#: + + +def poisson_quantiles(mu: np.ndarray, taus: Sequence[float]) -> np.ndarray: + """Discrete Poisson quantiles, shape ``(len(mu), len(taus))``.""" + mean = np.clip(np.asarray(mu, dtype=float), _TINY, None) + return np.stack([stats.poisson.ppf(tau, mean) for tau in _check_levels(taus)], axis=1) + + +def negative_binomial_family(dispersion: Optional[float] = None) -> NegativeBinomialFamily: + """Negative binomial family with log link; ``dispersion`` (alpha) is estimated when None.""" + return NegativeBinomialFamily(dispersion) + + +def tweedie_family(power: float = 1.5) -> TweedieFamily: + """Tweedie family (``1 < power < 2``) with log link.""" + return TweedieFamily(power) + + +def predict_count_quantiles(model, data: pd.DataFrame, taus: Sequence[float]) -> pd.DataFrame: + """Quantiles of a fitted Poisson, negative binomial or Tweedie ``StaticTAM`` on the rows of ``data``. + + Args: + model: A fitted StaticTAM whose loss has a ``quantiles`` method. + data: The rows to forecast. + taus: Quantile levels in (0, 1). + + Returns: + A DataFrame with one ``q`` column per level, non-negative, non-decreasing in the level. + """ + levels = _check_levels(taus) + predicted = model.predict(data) + mean = predicted[f"Estimated{model.target_col_}"].to_numpy(dtype=float) + strategy = model._reweighting_strategy_ + quantiles = strategy.quantiles(mean, levels) + return pd.DataFrame({f"q{tau}": quantiles[:, j] for j, tau in enumerate(levels)}, index=predicted.index) + + +# Poisson gets its quantiles here, without changing poisson_family() in _glm.py. +def poisson_family_with_quantiles() -> GLMFamily: + family = poisson_family() + family.quantiles = poisson_quantiles + return family diff --git a/src/tam/model/statistics/estimation/_distributional.py b/src/tam/model/statistics/estimation/_distributional.py index 40a7487..8c663fc 100644 --- a/src/tam/model/statistics/estimation/_distributional.py +++ b/src/tam/model/statistics/estimation/_distributional.py @@ -19,6 +19,7 @@ import numpy as np import pandas as pd from scipy import stats +from scipy.special import expit _TINY: float = 1e-12 # Mean of log(chi^2_1); de-biases an L2 fit on log(r^2) back to log(sigma^2). @@ -104,11 +105,148 @@ def init_config(model, formulas, loss, group_col, date_col, default_alpha_p, model.tail_family_ = None model.nu_ = None model._global_scale_variance_ = None + model._transform_params_ = {} # --- Helpers ---------------------------------------------------------------------------------------- +def require_positive_target(model, y: np.ndarray) -> None: + """Reject a training target with values <= 0 when the model works on log(y). + + Args: + model: A distributional or mixture StaticTAM. + y: The training target. + + Raises: + ValueError: If ``log_target`` is on and some finite values are <= 0 (they would be clipped to a tiny constant + before the log, i.e. turned into huge negative outliers). + """ + if not model._log_target_: + return + values = np.asarray(y, dtype=float) + finite = values[np.isfinite(values)] + n_bad = int(np.count_nonzero(finite <= 0)) + if n_bad: + raise ValueError( + f"{n_bad} of {finite.size} training target values are <= 0, but the model works on log(y) " + "(the default for a distributional or mixture StaticTAM). Pass dist_kwargs={'log_target': False} " + "to model the target on its own scale, or shift it to positive values." + ) + + +_TRANSFORMS = ("none", "log", "asinh", "logit") +# Offset from the bounds of a logit target, as a fraction of the interval: a value on a bound is nudged inside. +_LOGIT_EPS: float = 1e-6 +# Margin added on each side of the observed range when the bounds of a logit target are estimated. +_LOGIT_MARGIN: float = 0.05 + + +def resolve_target_transform(dist_kwargs: dict, mixture: bool = False) -> Tuple[str, bool, Optional[Tuple[float, float]]]: + """Read ``target_transform``, ``log_target`` and ``bounds`` of ``dist_kwargs``. + + Args: + dist_kwargs: The ``dist_kwargs`` of the model. + mixture: Whether the model is a mixture (only ``"none"`` and ``"log"`` are available there). + + Returns: + ``(transform, log_target, bounds)``: the transform name, ``transform == "log"``, and the given bounds or None. + + Raises: + ValueError: On an unknown transform, a contradiction with ``log_target``, invalid ``bounds``, or a mixture with + ``"asinh"`` or ``"logit"``. + """ + transform = dist_kwargs.get("target_transform") + if transform is None: + transform = "log" if bool(dist_kwargs.get("log_target", True)) else "none" + elif transform not in _TRANSFORMS: + raise ValueError(f"target_transform must be one of {_TRANSFORMS}; got {transform!r}") + elif "log_target" in dist_kwargs and bool(dist_kwargs["log_target"]) != (transform == "log"): + raise ValueError( + f"dist_kwargs target_transform={transform!r} contradicts log_target={dist_kwargs['log_target']!r}; give one of them." + ) + if mixture and transform in ("asinh", "logit"): + raise ValueError(f"target_transform={transform!r} is not available for a mixture; use 'log' or 'none'.") + bounds = dist_kwargs.get("bounds") + if bounds is not None: + if transform != "logit": + raise ValueError("bounds applies to target_transform='logit' only.") + low, high = float(bounds[0]), float(bounds[1]) + if not low < high: + raise ValueError(f"bounds must be (low, high) with low < high; got {bounds!r}") + bounds = (low, high) + return transform, transform == "log", bounds + + +def fit_target_transform(model, y: np.ndarray) -> None: + """Estimate and store the parameters of the target transform on the training target. + + ``"asinh"`` stores the median and a robust scale (normalised MAD, else the standard deviation, else 1); ``"logit"`` stores + the bounds (given, or the observed range widened by a margin). + + Args: + model: A distributional StaticTAM. + y: The training target. + + Raises: + ValueError: If a ``"logit"`` target has values outside the given bounds. + """ + values = np.asarray(y, dtype=float) + finite = values[np.isfinite(values)] + params: dict = {} + if model._target_transform_ == "asinh": + centre = float(np.median(finite)) + scale = float(1.4826 * np.median(np.abs(finite - centre))) + if not scale > 0.0: + scale = float(np.std(finite)) + params = {"centre": centre, "scale": scale if scale > 0.0 else 1.0} + elif model._target_transform_ == "logit": + bounds = model._transform_bounds_ + if bounds is None: + low, high = float(finite.min()), float(finite.max()) + margin = _LOGIT_MARGIN * (high - low if high > low else 1.0) + bounds = (low - margin, high + margin) + elif finite.min() < bounds[0] or finite.max() > bounds[1]: + n_out = int(np.count_nonzero((finite < bounds[0]) | (finite > bounds[1]))) + raise ValueError(f"{n_out} training target values lie outside the bounds {bounds} of the logit transform.") + params = {"bounds": bounds} + model._transform_params_ = params + + def _to_model_scale(model, y: np.ndarray) -> np.ndarray: - return np.log(np.clip(y, _TINY, None)) if model._log_target_ else np.asarray(y, dtype=float) + values = np.asarray(y, dtype=float) + transform = model._target_transform_ + if transform == "log": + return np.log(np.clip(values, _TINY, None)) + if transform == "asinh": + return np.arcsinh((values - model._transform_params_["centre"]) / model._transform_params_["scale"]) + if transform == "logit": + low, high = model._transform_params_["bounds"] + unit = np.clip((values - low) / (high - low), _LOGIT_EPS, 1.0 - _LOGIT_EPS) + return np.log(unit / (1.0 - unit)) + return values + + +def from_model_scale(model, z: np.ndarray, quantity: str = "the back-transformed values") -> np.ndarray: + """Back-transform model-scale values to the response scale (exact for a monotone transform, quantiles included). + + Args: + model: A distributional StaticTAM. + z: Values on the model scale. + quantity: Name used in the overflow warning of the log transform. + + Returns: + The values on the response scale. + """ + values = np.asarray(z, dtype=float) + transform = model._target_transform_ + if transform == "log": + return _to_response_scale(values, quantity) + if transform == "asinh": + with np.errstate(over="ignore"): + return model._transform_params_["centre"] + model._transform_params_["scale"] * np.sinh(values) + if transform == "logit": + low, high = model._transform_params_["bounds"] + return low + (high - low) * expit(values) + return values def _estimated(predicted: pd.DataFrame, internal_target: str) -> np.ndarray: @@ -157,6 +295,8 @@ def fit(model, data: pd.DataFrame, select: str = "fixed", smoothing stays fixed (Gaussian GCV is invalid for the Gamma GLM). """ working = data.copy() + require_positive_target(model, working[model.target_col_].to_numpy()) + fit_target_transform(model, working[model.target_col_].to_numpy()) working["__mu__"] = _to_model_scale(model, working[model.target_col_].to_numpy()) location_alpha = model._location_alpha_p_ if model._location_alpha_p_ is not None else model.default_alpha_p_ @@ -238,7 +378,7 @@ def _to_response_scale(model_scale: np.ndarray, quantity: str) -> np.ndarray: def predict_median(model, data: pd.DataFrame) -> np.ndarray: require_distributional(model, "predict_median") mu_hat, _ = mu_sigma(model, data) - return _to_response_scale(mu_hat, "the median (predict_median)") if model._log_target_ else mu_hat + return from_model_scale(model, mu_hat, "the median (predict_median)") def predict_quantile(model, data: pd.DataFrame, tau) -> np.ndarray: @@ -248,10 +388,7 @@ def predict_quantile(model, data: pd.DataFrame, tau) -> np.ndarray: columns = [] for level in tau_list: model_scale_quantile = mu_hat + sigma_hat * _standardized_ppf(model, float(level)) - columns.append( - _to_response_scale(model_scale_quantile, f"the quantile at level {float(level)} (predict_quantile)") - if model._log_target_ else model_scale_quantile - ) + columns.append(from_model_scale(model, model_scale_quantile, f"the quantile at level {float(level)} (predict_quantile)")) stacked = np.stack(columns, axis=1) return stacked[:, 0] if np.isscalar(tau) else stacked @@ -301,7 +438,7 @@ def anomaly_score(model, data: pd.DataFrame) -> pd.DataFrame: z_score = (observed - mu_hat) / sigma_hat upper_tail = _standardized_cdf(model, z_score) tail_pvalue = np.clip(2.0 * np.minimum(upper_tail, 1.0 - upper_tail), _TINY, 1.0) - median = _to_response_scale(mu_hat, "the median (anomaly_score)") if model._log_target_ else mu_hat + median = from_model_scale(model, mu_hat, "the median (anomaly_score)") return pd.DataFrame({ "predicted_median": median, "z_score": z_score, diff --git a/src/tam/model/statistics/estimation/_factory.py b/src/tam/model/statistics/estimation/_factory.py index 4391893..15d6f01 100644 --- a/src/tam/model/statistics/estimation/_factory.py +++ b/src/tam/model/statistics/estimation/_factory.py @@ -10,23 +10,34 @@ """ from __future__ import annotations +from typing import Optional + from ._base_strategy import ReweightingStrategy -from ._glm import gaussian_family, gamma_family, poisson_family, binomial_family +from ._glm import gaussian_family, gamma_family, binomial_family +from ._count_families import negative_binomial_family, poisson_family_with_quantiles, tweedie_family from ._expectile import ExpectileLoss from ._robust import HuberLoss, StudentTLoss _STRATEGY_ALIASES = { "l2": gaussian_family, "gaussian": gaussian_family, "normal": gaussian_family, - "gamma": gamma_family, "poisson": poisson_family, "binomial": binomial_family, + "gamma": gamma_family, "poisson": poisson_family_with_quantiles, "binomial": binomial_family, } #: -def build_strategy(loss: str, tau: float = 0.5, nu: float = 4.0, delta: float = 1.345) -> ReweightingStrategy: - """Map a loss name to a ReweightingStrategy.""" +def build_strategy(loss: str, tau: float = 0.5, nu: float = 4.0, delta: float = 1.345, + power: float = 1.5, dispersion: Optional[float] = None) -> ReweightingStrategy: + """Map a loss name to a ReweightingStrategy. + + ``power`` is the Tweedie power (1 < power < 2); ``dispersion`` fixes the negative binomial alpha (estimated when None). + """ key = loss.lower().replace("-", "_") if key in _STRATEGY_ALIASES: return _STRATEGY_ALIASES[key]() + if key in ("tweedie", "compound_poisson"): + return tweedie_family(power) + if key in ("negative_binomial", "negbin", "nb"): + return negative_binomial_family(dispersion) if key == "expectile": return ExpectileLoss(tau) if key == "huber": @@ -34,6 +45,6 @@ def build_strategy(loss: str, tau: float = 0.5, nu: float = 4.0, delta: float = if key in ("student_t", "studentt", "student"): return StudentTLoss(nu) raise ValueError( - f"Unknown loss '{loss}'. Supported: l2/gaussian, gamma, poisson, binomial, expectile, huber, student_t." + f"Unknown loss '{loss}'. Supported: l2/gaussian, gamma, poisson, negative_binomial, tweedie, binomial, expectile, huber, student_t." ) #: diff --git a/src/tam/model/statistics/estimation/_mixture.py b/src/tam/model/statistics/estimation/_mixture.py index 9fc7220..2c4cc65 100644 --- a/src/tam/model/statistics/estimation/_mixture.py +++ b/src/tam/model/statistics/estimation/_mixture.py @@ -23,6 +23,7 @@ from tam.common.utils import _ensure_dummies, _balance_groups from ..._base import BaseTAM +from ._distributional import require_positive_target _TINY: float = 1e-12 @@ -129,6 +130,7 @@ def _prepare_for_em(model, data, with_target: bool): def fit(model, data): """Fit a K-component Gaussian mixture by EM; the M-step is the responsibility-weighted atom.""" working = data.copy() + require_positive_target(model, working[model.target_col_].to_numpy()) working[model.target_col_] = _to_model_scale(model, working[model.target_col_].to_numpy()) # Establish normalization, effects and a baseline L2 fit on the (log-)target scale. diff --git a/src/tam/model/statistics/risk/conformal.py b/src/tam/model/statistics/risk/conformal.py index 6ebb647..304fcb9 100644 --- a/src/tam/model/statistics/risk/conformal.py +++ b/src/tam/model/statistics/risk/conformal.py @@ -11,7 +11,7 @@ calibration, and Adaptive Conformal Inference for streaming coverage under drift. It duck-types the model, needing only predict_quantile / anomaly_score / _mu_sigma / _to_model_scale / -target_col_ / _log_target_ - the surface any distributional StaticTAM exposes. The conformal mechanics (the +target_col_ / _from_model_scale - the surface any distributional StaticTAM exposes. The conformal mechanics (the finite-sample quantile, the p-value, the ACI level update) are delegated to the static SafetyTAM and the streaming aci module; this class only computes the model-specific nonconformity scores (a studentized residual and the CQR score) and assembles the outputs. One source of truth for conformal. @@ -117,8 +117,8 @@ def aci_intervals(self, ordered_data: pd.DataFrame, gamma: float = 0.02) -> pd.D lower = mu_hat[position] - sigma_hat[position] * radius upper = mu_hat[position] + sigma_hat[position] * radius inside = lower <= observed[position] <= upper - lowers.append(np.exp(lower) if self.model._log_target_ else lower) - uppers.append(np.exp(upper) if self.model._log_target_ else upper) + lowers.append(float(self.model._from_model_scale(np.array([lower]), "the lower bound")[0])) + uppers.append(float(self.model._from_model_scale(np.array([upper]), "the upper bound")[0])) levels.append(alpha_t) alpha_t = update_risk_level(alpha_t, 0.0 if inside else 1.0, self.alpha, gamma) return pd.DataFrame({"lower": lowers, "upper": uppers, "alpha_t": levels}, index=ordered_data.index) diff --git a/tests/regression/RESULTS.txt b/tests/regression/RESULTS.txt index 476f6dc..a9e65ac 100644 --- a/tests/regression/RESULTS.txt +++ b/tests/regression/RESULTS.txt @@ -1819,7 +1819,7 @@ 15_losses.student_t_nu10.q75 = 56478.5 15_losses.student_t_nu10.rmse = 1206.01 15_losses.student_t_nu10.std = 10248.6 -15_losses.unknown_loss.message = Unknown loss 'no_such_loss'. Supported: l2/gaussian, gamma, poisson, b +15_losses.unknown_loss.message = Unknown loss 'no_such_loss'. Supported: l2/gaussian, gamma, poisson, n 16_distributional.anomaly.anomaly_score.bias = 1.11479 16_distributional.anomaly.anomaly_score.max = 20.5691 16_distributional.anomaly.anomaly_score.mean = 1.11479 @@ -2725,7 +2725,7 @@ 29_errors.model.predict_before_fit = RuntimeError: Model must be fitted before prediction. 29_errors.model.quantiles_on_a_point_model = RuntimeError: predict_quantiles() is only available on a distributional StaticTAM (a dict formula such a 29_errors.model.unknown_group_at_predict = ValueError: TAM [Data Error]: group(s) [np.int64(99)] of column 'hour' were not seen in training; the -29_errors.model.unknown_loss = ValueError: Unknown loss 'no_such_loss'. Supported: l2/gaussian, gamma, poisson, binomial, expectile, +29_errors.model.unknown_loss = ValueError: Unknown loss 'no_such_loss'. Supported: l2/gaussian, gamma, poisson, negative_binomial, tw 29_errors.opera.no_expert = ValueError: No valid experts found. Ensure syntax uses 'l(expert_name)'. 29_errors.opera.unknown_algorithm = ValueError: Algorithm must be 'EWA' or 'MLPOL'. 30_warnings.adaptive_in_range.classes = TAMExtrapolationWarning diff --git a/tests/test_adaptive_quantiles.py b/tests/test_adaptive_quantiles.py new file mode 100644 index 0000000..b95a493 --- /dev/null +++ b/tests/test_adaptive_quantiles.py @@ -0,0 +1,137 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +r""" +``AdaptiveTAM.predict_quantiles``: online quantiles of a distributional base, whose scale is rescaled by the standardized residuals of +the adaptive location, using past residuals only. +""" +import numpy as np +import pandas as pd +import pytest + +import tam as ta +from tam.model import adaptative + +N, BREAK, BASE_ROWS, WINDOW, UPDATE = 2000, 600, 500, 60, 10 +_DIST = {"log_target": False, "tail_family": "normal"} + + +def _data(seed=0, n=N, break_at=BREAK, sd_after=3.0): + """y = 2 + 1.5 x + noise, with a noise standard deviation of 1 that becomes ``sd_after`` at ``break_at``.""" + rng = np.random.default_rng(seed) + x = rng.normal(0.0, 1.0, n) + sd = np.where(np.arange(n) < break_at, 1.0, sd_after) + return pd.DataFrame({"date": pd.date_range("2020-01-01", periods=n, freq="D"), "x": x, "y": 2.0 + 1.5 * x + sd * rng.normal(size=n)}) + + +def _base(df): + base = ta.StaticTAM({"mu": "y ~ l(x)", "sigma": "y ~ l(x)"}, loss={"mu": "l2", "sigma": "gamma"}, date_col="date", dist_kwargs=_DIST) + return base.fit(df.iloc[:BASE_ROWS]) + + +def _adaptive(df, horizon=1): + return ta.AdaptiveTAM(base_model=_base(df), adaptive_formula="Residualy ~ l(effect_x)", update_interval_periods=UPDATE, + training_window_periods=WINDOW, steps_per_period=1, horizon_steps=horizon, date_col="date", default_alpha_p=-6.0) + + +def test_a_plain_base_has_no_quantiles(): + df = _data() + plain = ta.StaticTAM("y ~ l(x)", date_col="date").fit(df.iloc[:BASE_ROWS]) + model = ta.AdaptiveTAM(base_model=plain, adaptive_formula="Residualy ~ l(effect_x)", update_interval_periods=UPDATE, + training_window_periods=WINDOW, steps_per_period=1, date_col="date", default_alpha_p=-6.0) + with pytest.raises(ValueError, match="distributional"): + model.predict_quantiles(df, (0.1, 0.9)) + + +def test_levels_outside_zero_one_are_rejected(): + df = _data() + with pytest.raises(ValueError, match="taus"): + _adaptive(df).predict_quantiles(df, (0.0, 0.5)) + + +def test_the_interval_follows_a_variance_break(): + df = _data() + quantiles = _adaptive(df).predict_quantiles(df, (0.1, 0.5, 0.9)) + late = (df["date"] >= df["date"].iloc[BREAK + 2 * WINDOW]).to_numpy() + y = df["y"].to_numpy() + inside = (y >= quantiles["q0.1"].to_numpy()) & (y <= quantiles["q0.9"].to_numpy()) + assert abs(float(inside[late].mean()) - 0.8) < 0.03 + + +def test_without_the_online_scale_the_interval_is_too_narrow_after_the_break(monkeypatch): + """The base alone (s = 1) learned a noise standard deviation of 1: after the break its 80% interval covers about a third.""" + monkeypatch.setattr(adaptative, "_online_scale", lambda standardized, *args: np.ones(len(standardized))) + df = _data() + quantiles = _adaptive(df).predict_quantiles(df, (0.1, 0.9)) + late = (df["date"] >= df["date"].iloc[BREAK + 2 * WINDOW]).to_numpy() + y = df["y"].to_numpy() + inside = (y >= quantiles["q0.1"].to_numpy()) & (y <= quantiles["q0.9"].to_numpy()) + assert float(inside[late].mean()) < 0.5 + + +def test_quantiles_do_not_decrease_with_the_level(): + df = _data() + quantiles = _adaptive(df).predict_quantiles(df, (0.05, 0.25, 0.5, 0.75, 0.95)) + values = quantiles.to_numpy() + assert (np.diff(values, axis=1) >= 0).all() and np.isfinite(values).all() + assert list(quantiles.columns) == ["q0.05", "q0.25", "q0.5", "q0.75", "q0.95"] + + +@pytest.mark.parametrize("horizon", [1, 3]) +def test_no_look_ahead_a_later_target_does_not_move_an_earlier_quantile(horizon): + df = _data() + before = _adaptive(df, horizon).predict_quantiles(df, (0.1, 0.9)) + changed = df.copy() + late = np.arange(len(df)) >= 1200 + changed.loc[late, "y"] = np.random.default_rng(1).normal(size=int(late.sum())) * 50.0 + after = _adaptive(df, horizon).predict_quantiles(changed, (0.1, 0.9)) + np.testing.assert_allclose(before.to_numpy()[:1200], after.to_numpy()[:1200], rtol=0, atol=1e-9) + assert not np.allclose(before.to_numpy()[1300:], after.to_numpy()[1300:]) + + +def test_with_a_scale_of_one_the_quantiles_are_the_base_sigma_around_the_adapted_location(monkeypatch): + monkeypatch.setattr(adaptative, "_online_scale", lambda standardized, *args: np.ones(len(standardized))) + df = _data() + model = _adaptive(df) + quantiles = model.predict_quantiles(df, (0.1, 0.5, 0.9)) + frame = model.predictions_ + location = frame["AdaptedEstimated__mu__"].to_numpy() + _, sigma = model.distributional_base_._mu_sigma(frame) + np.testing.assert_allclose(quantiles["q0.5"].to_numpy(), location, atol=1e-9) + z = 1.2815515655446004 + np.testing.assert_allclose(quantiles["q0.9"].to_numpy() - location, sigma * z, rtol=1e-9) + np.testing.assert_allclose(location - quantiles["q0.1"].to_numpy(), sigma * z, rtol=1e-9) + + +def test_online_scale_uses_only_residuals_known_at_the_origin(): + # window 10, interval 5, horizon 2: the first origin is row 11, then 16, 21, ...; an origin o reads the rows o-10 ... o-2 + residual = np.ones(40) + residual[15] = 100.0 # known from the origin 17 on, so first read at the origin 21 + scale = adaptative._online_scale(residual, np.zeros(40), np.arange(40), 10, 5, 2) + assert (scale[:11] == 1.0).all() # before the first origin: the base scale + assert (scale[11:21] == 1.0).all() # the spike at row 15 is not yet readable at the origins 11 and 16 + assert (scale[21:26] > 10.0).all() + + +def test_online_scale_keeps_the_base_scale_with_too_few_residuals(): + residual = np.full(30, np.nan) + residual[:3] = 2.0 + scale = adaptative._online_scale(residual, np.zeros(30), np.arange(30), 10, 5, 1) + assert (scale == 1.0).all() + + +def test_online_scale_works_group_by_group(): + residual = np.concatenate([np.full(30, 1.0), np.full(30, 3.0)]) + group = np.repeat([0, 1], 30) + scale = adaptative._online_scale(residual, group, np.tile(np.arange(30), 2), 10, 5, 1) + assert np.allclose(scale[20:30], 1.0) and np.allclose(scale[50:60], 3.0) + + +def test_a_plain_adaptive_model_is_unchanged_by_the_distributional_option(): + df = _data() + plain = ta.StaticTAM("y ~ l(x)", date_col="date").fit(df.iloc[:BASE_ROWS]) + model = ta.AdaptiveTAM(base_model=plain, adaptive_formula="Residualy ~ l(effect_x)", update_interval_periods=UPDATE, + training_window_periods=WINDOW, steps_per_period=1, date_col="date", default_alpha_p=-6.0) + assert model.distributional_base_ is None + assert "EstimatedResidualy" in model.predict_online(df).columns diff --git a/tests/test_count_families.py b/tests/test_count_families.py new file mode 100644 index 0000000..8332228 --- /dev/null +++ b/tests/test_count_families.py @@ -0,0 +1,152 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +r""" +Negative binomial and Tweedie families of :class:`StaticTAM` (``loss="negative_binomial"`` and ``loss="tweedie"``, log link), and +``predict_quantiles`` for the count and zero-inflated-positive models: discrete quantiles for Poisson and negative binomial, the +compound Poisson-gamma law for Tweedie. Quantiles are non-negative by construction. +""" +import numpy as np +import pandas as pd +import pytest +from scipy import stats + +import tam as ta +from tam.model.statistics.estimation import build_strategy + +N = 20000 + + +def _frame(y, x): + return pd.DataFrame({"date": pd.date_range("2020-01-01", periods=len(y), freq="h"), "x": x, "y": y}) + + +def _negative_binomial(seed=0, n=N, theta=2.0, intercept=4.0, slope=0.5): + rng = np.random.default_rng(seed) + x = rng.normal(size=n) + mu = np.exp(intercept + slope * x) + return _frame(rng.negative_binomial(theta, theta / (theta + mu)), x), mu + + +def _tweedie(seed=0, n=N, power=1.5, phi=2.0, intercept=-3.0, slope=1.0): + """Compound Poisson-gamma draws with mean mu = exp(intercept + slope x) and variance phi mu^power.""" + rng = np.random.default_rng(seed) + x = rng.normal(size=n) + mu = np.exp(intercept + slope * x) + lam = mu ** (2 - power) / (phi * (2 - power)) + shape = (2 - power) / (power - 1) + scale = phi * (power - 1) * mu ** (power - 1) + counts = rng.poisson(lam) + y = np.array([rng.gamma(shape * c, s) if c else 0.0 for c, s in zip(counts, scale)]) + return _frame(y, x), mu + + +def _model(loss, **loss_kwargs): + return ta.StaticTAM("y ~ l(x)", date_col="date", loss=loss, loss_kwargs=loss_kwargs or None) + + +def test_negative_binomial_recovers_the_dispersion(): + df, _ = _negative_binomial() + model = _model("negative_binomial").fit(df) + assert abs(model.dispersion_ - 0.5) / 0.5 < 0.10 # the dispersion alpha = 1 / theta + coef_ratio = model.predict(df)["Estimatedy"].mean() / df["y"].mean() + assert abs(coef_ratio - 1.0) < 0.02 + + +def test_negative_binomial_interval_covers(): + train, _ = _negative_binomial(seed=0) + test, _ = _negative_binomial(seed=1) + model = _model("negative_binomial").fit(train) + quantiles = model.predict_quantiles(test, (0.1, 0.5, 0.9)) + y = test["y"].to_numpy() + inside = (y >= quantiles["q0.1"].to_numpy()) & (y <= quantiles["q0.9"].to_numpy()) + assert abs(float(inside.mean()) - 0.8) < 0.03 + + +def test_a_given_dispersion_is_kept_and_not_estimated(): + df, _ = _negative_binomial() + model = _model("negative_binomial", dispersion=0.3).fit(df) + assert model.dispersion_ == 0.3 + + +def test_poisson_quantiles_are_the_discrete_poisson_quantiles(): + rng = np.random.default_rng(2) + x = rng.normal(size=5000) + df = _frame(rng.poisson(np.exp(2.0 + 0.4 * x)), x) + model = _model("poisson").fit(df) + mean = model.predict(df)["Estimatedy"].to_numpy() + quantiles = model.predict_quantiles(df, (0.1, 0.5, 0.9)) + for tau in (0.1, 0.5, 0.9): + np.testing.assert_array_equal(quantiles[f"q{tau}"].to_numpy(), stats.poisson.ppf(tau, mean)) + assert model.dispersion_ is None + + +def test_tweedie_recovers_the_mean_and_the_dispersion(): + df, mu = _tweedie() + model = _model("tweedie", power=1.5).fit(df) + assert abs(model.predict(df)["Estimatedy"].mean() / df["y"].mean() - 1.0) < 0.05 + assert abs(model.dispersion_ - 2.0) / 2.0 < 0.15 + assert (df["y"] == 0).mean() > 0.6 # a mostly-zero target + + +def test_tweedie_quantiles_are_non_negative_and_keep_the_mass_at_zero(): + train, _ = _tweedie(seed=0) + test, _ = _tweedie(seed=1, n=4000) + model = _model("tweedie", power=1.5).fit(train) + taus = (0.05, 0.25, 0.5, 0.75, 0.95) + quantiles = model.predict_quantiles(test, taus) + values = quantiles.to_numpy() + assert values.min() >= 0.0 and np.isfinite(values).all() + assert (np.diff(values, axis=1) >= -1e-9).all() + mean = model.predict(test)["Estimatedy"].to_numpy() + p_zero = np.exp(-(mean ** 0.5) / (model.dispersion_ * 0.5)) # lambda = mu^(2-p) / (phi (2-p)), p = 1.5 + sure_zero = p_zero > 0.5 + assert sure_zero.sum() > 100 + assert (quantiles["q0.5"].to_numpy()[sure_zero] == 0.0).all() + assert (quantiles["q0.5"].to_numpy()[p_zero < 0.45] > 0.0).all() + + +def test_tweedie_interval_covers(): + train, _ = _tweedie(seed=0) + test, _ = _tweedie(seed=3, n=4000) + model = _model("tweedie", power=1.5).fit(train) + quantiles = model.predict_quantiles(test, (0.1, 0.9)) + y = test["y"].to_numpy() + assert float(np.mean(y <= quantiles["q0.9"].to_numpy())) > 0.88 + assert float(np.mean(y <= quantiles["q0.9"].to_numpy())) < 0.94 + assert float(np.mean(y >= quantiles["q0.1"].to_numpy())) > 0.88 + + +@pytest.mark.parametrize("loss", ["negative_binomial", "tweedie"]) +def test_a_negative_target_is_rejected(loss): + df, _ = _negative_binomial(n=500) + df.loc[0, "y"] = -1 + with pytest.raises(ValueError, match="negative"): + _model(loss).fit(df) + + +@pytest.mark.parametrize("power", [1.0, 2.0, 0.5, 3.0]) +def test_the_tweedie_power_must_lie_in_one_two(power): + with pytest.raises(ValueError, match="power"): + build_strategy("tweedie", power=power) + + +def test_the_dispersion_must_be_positive(): + with pytest.raises(ValueError, match="dispersion"): + build_strategy("negative_binomial", dispersion=0.0) + + +def test_other_losses_have_no_discrete_quantiles_and_no_dispersion(): + df, _ = _negative_binomial(n=500) + model = ta.StaticTAM("y ~ l(x)", date_col="date").fit(df) + assert model.dispersion_ is None + with pytest.raises(RuntimeError, match="predict_quantiles"): + model.predict_quantiles(df, (0.5,)) + + +def test_the_mean_of_the_new_families_is_a_log_link_glm(): + df, mu = _negative_binomial(n=5000) + fitted = _model("negative_binomial").fit(df).predict(df)["Estimatedy"].to_numpy() + assert (fitted > 0).all() + assert np.corrcoef(np.log(fitted), np.log(mu))[0, 1] > 0.99 diff --git a/tests/test_cyclic_fixed_domain.py b/tests/test_cyclic_fixed_domain.py new file mode 100644 index 0000000..0ce07d6 --- /dev/null +++ b/tests/test_cyclic_fixed_domain.py @@ -0,0 +1,130 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +r""" +A Fourier term can declare the period of its input: ``f(toy, cyclic=True, period=(0, 1))``. + +Every feature is normalised to [-1, 1] with the minimum and maximum of its training rows, and a cyclic term takes [-1, 1] as one +period. With a time of year in [0, 1] and a history from January to June, the training range [0, 0.5] became one whole period: +January 1st and June 30th were the same point of the cycle and the forecast of July started the year again. A declared period is +the range the feature is normalised on instead. Without ``period`` nothing changes (the training range is the period, as before). +""" +import numpy as np +import pandas as pd +import pytest + +import tam as ta +from tam.common.utils import parse_formula_to_terms +from tam.model._data import _fit_normalization_params +from tam.model.spectrum import fixed_ranges +from tam.model.spectrum._factory import create_effects_from_parsed_terms + +AMPLITUDE = 3.0 + + +def _daily_yearly_sine(days: int, start: str = "2023-01-01", noise: float = 0.0, seed: int = 0) -> pd.DataFrame: + dates = pd.date_range(start, periods=days, freq="D") + toy = (dates.dayofyear.to_numpy() - 1) / 365.0 # time of year on [0, 1), the usual coding + rng = np.random.default_rng(seed) + y = 10.0 + AMPLITUDE * np.sin(2.0 * np.pi * toy) + noise * rng.standard_normal(days) + return pd.DataFrame({"date": dates, "toy": toy, "y": y}) + + +def _half_year_error(formula: str) -> float: + data = _daily_yearly_sine(365) + train, test = data.iloc[:182], data.iloc[182:] # January to June, then July to December + model = ta.StaticTAM(formula, date_col="date").fit(train) + return float(np.abs(model.predict(test)["Estimatedy"].to_numpy() - test["y"].to_numpy()).max()) + + +def test_declared_period_reads_half_a_year_as_half_a_cycle(): + assert _half_year_error("y ~ f(toy, m=1, s=1, cyclic=True, period=(0, 1))") <= 0.1 * AMPLITUDE + + +def test_without_period_half_a_year_is_taken_as_the_whole_cycle(): + # The default is unchanged: the training range is the period, which is wrong for a partial cycle (documented behaviour). + assert _half_year_error("y ~ f(toy, m=1, s=1, cyclic=True)") > AMPLITUDE + + +def test_values_outside_the_period_wrap_around(): + data = _daily_yearly_sine(365 * 2, noise=0.1) + model = ta.StaticTAM("y ~ f(toy, m=3, s=1, cyclic=True, period=(0, 1))", date_col="date").fit(data) + probe = data.iloc[:2].assign(toy=[0.02, 1.02]) # 1.02 is early January of the next cycle + inside, outside = model.predict(probe)["Estimatedy"].to_numpy() + assert inside == pytest.approx(outside, abs=1e-9) + + +def test_adaptive_windows_shorter_than_the_cycle_use_the_declared_period(): + # Each 60-day window was normalised on its own rows: two months were read as a whole year, even with years of data. + data = _daily_yearly_sine(365 * 2) + data["grp"] = "a" + model = ta.AdaptiveTAM(adaptive_formula="y ~ f(toy, m=1, s=1, cyclic=True, period=(0, 1))", update_interval_periods=7, + training_window_periods=60, steps_per_period=1, group_col="grp", date_col="date") + scored = model.predict_online(data).dropna(subset=["AdaptedEstimatedy"]) + error = np.abs(scored["AdaptedEstimatedy"].to_numpy() - scored["y"].to_numpy()) + assert np.median(error) <= 0.05 * AMPLITUDE + + +def _effects(formula: str): + return create_effects_from_parsed_terms(parse_formula_to_terms(formula)[1], token_values={}, default_alpha_p=-9.0, data_info={}) + + +def test_fixed_ranges_rules(): + assert fixed_ranges(_effects("y ~ f(toy, m=2, cyclic=True, period=(0, 1))")) == {"toy": (0.0, 1.0)} + assert fixed_ranges(_effects("y ~ f(toy, m=2, cyclic=True)")) == {} # no period: min/max as before + assert fixed_ranges(_effects("y ~ f(hour, m=2, cyclic=True, period=(0, 24)) + l(hour)")) == {"hour": (0.0, 24.0)} + assert fixed_ranges(_effects("y ~ te(f(toy, m=2, cyclic=True, period=(0, 1)), s(x))")) == {"toy": (0.0, 1.0)} + assert fixed_ranges(_effects("y ~ f(toy, m=2, cyclic=True, period=(0, 1)) + c(toy, n_cat=4)")) == {} + with pytest.raises(ValueError, match="Two periods"): + fixed_ranges(_effects("y ~ f(toy, m=2, period=(0, 1)) + f(toy, m=3, period=(0, 2))")) + with pytest.raises(ValueError, match="period"): + _effects("y ~ f(toy, m=2, period=(1, 0))") + + +def test_normalization_uses_the_declared_range(): + frame = pd.DataFrame({"g": 0, "toy": np.linspace(0.1, 0.4, 50), "x": np.linspace(5.0, 9.0, 50)}) + params, _ = _fit_normalization_params(frame, ["toy", "x"], "g", fixed_ranges={"toy": (0.0, 1.0)}) + assert params[0]["min"]["toy"] == 0.0 and params[0]["max"]["toy"] == 1.0 + assert params[0]["min"]["x"] == 5.0 and params[0]["max"]["x"] == 9.0 # other features unchanged + + +def test_integer_hours_need_the_period_to_tell_23h_from_0h(): + # Hours 0..23 normalised on their own range put 23 on +1 and 0 on -1, the same point of a cyclic basis. + hours = np.tile(np.arange(24), 200).astype(float) + profile = np.where(hours == 23, 8.0, np.where(hours == 0, 2.0, 5.0 + np.sin(2 * np.pi * hours / 24))) + data = pd.DataFrame({"hour": hours, "y": profile + 0.1 * np.random.default_rng(0).standard_normal(len(hours))}) + probe = pd.DataFrame({"hour": [0.0, 23.0], "y": [0.0, 0.0]}) + glued = ta.StaticTAM("y ~ f(hour, m=10, s=1, cyclic=True)").fit(data).predict(probe)["Estimatedy"].to_numpy() + apart = ta.StaticTAM("y ~ f(hour, m=10, s=1, cyclic=True, period=(0, 24))").fit(data).predict(probe)["Estimatedy"].to_numpy() + assert glued[0] == pytest.approx(glued[1], abs=1e-6) + assert apart[1] - apart[0] > 3.0 + + +# --- period="auto": (first value, last value + one step), the step read from the training rows ------------------------ +def test_auto_period_tells_23h_from_0h(): + hours = np.tile(np.arange(24), 200).astype(float) + profile = np.where(hours == 23, 8.0, np.where(hours == 0, 2.0, 5.0 + np.sin(2 * np.pi * hours / 24))) + data = pd.DataFrame({"hour": hours, "y": profile + 0.1 * np.random.default_rng(0).standard_normal(len(hours))}) + probe = pd.DataFrame({"hour": [0.0, 23.0], "y": [0.0, 0.0]}) + auto = ta.StaticTAM("y ~ f(hour, m=10, s=1, cyclic=True, period='auto')").fit(data).predict(probe)["Estimatedy"].to_numpy() + given = ta.StaticTAM("y ~ f(hour, m=10, s=1, cyclic=True, period=(0, 24))").fit(data).predict(probe)["Estimatedy"].to_numpy() + np.testing.assert_allclose(auto, given, atol=1e-9) + + +@pytest.mark.parametrize("values, expected", [ + (np.arange(24.0), (0.0, 24.0)), # hours + (np.arange(1.0, 13.0), (1.0, 13.0)), # months + (np.arange(365) / 365.0, (0.0, 1.0)), # time of year, one value per day +]) +def test_auto_period_is_first_value_to_last_plus_one_step(values, expected): + from tam.model._data import _resolve_fixed_ranges + frame = pd.DataFrame({"x": np.tile(values, 3)}) + low, high = _resolve_fixed_ranges(frame, {"x": "auto"})["x"] + assert low == pytest.approx(expected[0]) and high == pytest.approx(expected[1]) + + +def test_auto_period_in_fixed_ranges_and_conflicts(): + assert fixed_ranges(_effects("y ~ f(hour, m=2, cyclic=True, period='auto')")) == {"hour": "auto"} + with pytest.raises(ValueError, match="Two periods"): + fixed_ranges(_effects("y ~ f(hour, m=2, period='auto') + f(hour, m=3, period=(0, 24))")) diff --git a/tests/test_distributional.py b/tests/test_distributional.py index dbe94ff..3879a10 100644 --- a/tests/test_distributional.py +++ b/tests/test_distributional.py @@ -188,3 +188,42 @@ def test_scale_shrinkage_monotonically_shrinks_scale_dispersion(): dispersions.append(float(np.std(sigma_hat))) assert all(later <= earlier + 1e-9 for earlier, later in zip(dispersions, dispersions[1:])) assert dispersions[-1] < dispersions[0] + + + +# --- log_target on a target that is not positive ------------------------------------------------------------------- +def _signed_frame(rng, n=400): + x = rng.uniform(0.0, 1.0, n) + return pd.DataFrame({"x": x, "y": 2.0 * x - 1.0 + 0.3 * rng.standard_normal(n)}) # about half the values are negative + + +def test_log_target_rejects_non_positive_training_target(): + # The default log transform used to clip y <= 0 to a tiny constant before taking the log, so every such value became + # a huge negative outlier and the fitted location and scale were meaningless, without any message. + frame = _signed_frame(np.random.default_rng(0)) + n_bad = int((frame["y"] <= 0).sum()) + with pytest.raises(ValueError, match=rf"{n_bad} of {len(frame)} .*log_target"): + _distributional().fit(frame) + + +def test_log_target_false_fits_a_signed_target(): + frame = _signed_frame(np.random.default_rng(1)) + model = _distributional(dist_kwargs={"log_target": False}).fit(frame) + q = model.predict_quantiles(frame.head(5), taus=(0.1, 0.5, 0.9)) + assert np.isfinite(q.to_numpy()).all() + assert (q["q0.1"] < q["q0.9"]).all() + + +def test_mixture_log_target_rejects_non_positive_training_target(): + frame = _signed_frame(np.random.default_rng(2)) + with pytest.raises(ValueError, match="log_target"): + StaticTAM("y ~ s(x)", mixture_components=2).fit(frame) + + +def test_cdf_of_a_non_positive_observation_still_works(): + # Only the training target is checked: the CDF of a new observation y <= 0 under a log-scale law is ~0, a valid answer. + rng = np.random.default_rng(3) + frame, _, _ = _heteroscedastic_lognormal(rng, 2000) + model = _distributional().fit(frame) + probe = frame.head(3).assign(y=[-1.0, 0.0, 1e-9]) + assert np.allclose(model.cdf(probe), 0.0, atol=1e-6) diff --git a/tests/test_kalman_quantiles.py b/tests/test_kalman_quantiles.py new file mode 100644 index 0000000..85cbc6b --- /dev/null +++ b/tests/test_kalman_quantiles.py @@ -0,0 +1,116 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +r""" +``KalmanTAM.predict_quantiles``: Gaussian quantiles from the predictive variance of the filter, ``x' P x + R`` in the target scale, +for the one-step forecast and for the delayed state of a multi-step horizon. +""" +import numpy as np +import pandas as pd +import pytest +import torch + +import tam as ta + +N, SIGMA = 4000, 2.0 +CALIBRATION = 365 + + +def _frame(seed=0, n=N, groups=1): + rng = np.random.default_rng(seed) + frames = [] + for g in range(groups): + x = rng.normal(size=n) + frames.append(pd.DataFrame({"date": pd.date_range("2020-01-01", periods=n, freq="D"), "g": g, "x": x, + "y": 5.0 + 2.0 * x + SIGMA * rng.normal(size=n)})) + return pd.concat(frames, ignore_index=True) + + +def _model(df, observation_noise_var=None, horizon=1, groups=1): + """A tracker with a nearly constant coefficient; R is the true noise variance in the standardised target space.""" + kwargs = dict(kalman_formula="y ~ l(x)", group_col="g" if groups > 1 else None, date_col="date", horizon_steps=horizon, + process_noise_var=1e-8, calibration_steps=CALIBRATION) + if observation_noise_var is None: + probe = ta.KalmanTAM(**kwargs) + probe.prepare_data(df) + scale = float(next(iter(probe._reference_scale_.values()))) + observation_noise_var = SIGMA ** 2 / scale ** 2 + return ta.KalmanTAM(observation_noise_var=observation_noise_var, **kwargs) + + +def test_the_interval_covers_when_the_noise_is_known(): + df = _frame() + quantiles = _model(df).predict_quantiles(df, (0.1, 0.5, 0.9)) + y = df["y"].to_numpy()[CALIBRATION:] + inside = (y >= quantiles["q0.1"].to_numpy()[CALIBRATION:]) & (y <= quantiles["q0.9"].to_numpy()[CALIBRATION:]) + assert abs(float(inside.mean()) - 0.8) < 0.03 + + +def test_the_median_is_the_online_forecast(): + df = _frame() + model = _model(df) + quantiles = model.predict_quantiles(df, (0.1, 0.5, 0.9)) + forecast = model.predict_online(df)["KalmanAdapted_y"].to_numpy() + np.testing.assert_allclose(quantiles["q0.5"].to_numpy(), forecast, rtol=0, atol=1e-9) + + +def test_quantiles_do_not_decrease_with_the_level_and_the_columns_are_named_by_level(): + df = _frame(n=1000) + quantiles = _model(df).predict_quantiles(df, (0.05, 0.25, 0.5, 0.75, 0.95)) + assert list(quantiles.columns) == ["q0.05", "q0.25", "q0.5", "q0.75", "q0.95"] + assert (np.diff(quantiles.to_numpy(), axis=1) > 0).all() + + +def test_levels_outside_zero_one_are_rejected(): + df = _frame(n=600) + with pytest.raises(ValueError, match="taus"): + _model(df).predict_quantiles(df, (0.5, 1.0)) + + +def test_no_look_ahead_a_later_target_does_not_move_an_earlier_quantile(): + df = _frame() + model = _model(df) + before = model.predict_quantiles(df, (0.1, 0.9)).to_numpy() + changed = df.copy() + changed.loc[2000:, "y"] = np.random.default_rng(1).normal(size=len(df) - 2000) * 40.0 + after = _model(df).predict_quantiles(changed, (0.1, 0.9)).to_numpy() + np.testing.assert_allclose(before[:2000], after[:2000], rtol=0, atol=1e-9) + assert not np.allclose(before[2100:], after[2100:]) + + +def test_the_predictive_variance_does_not_depend_on_the_target(): + """The covariance of a linear Gaussian filter is a function of the design and the noise, not of the observations.""" + df = _frame(n=1500) + model = _model(df) + shifted = df.assign(y=df["y"] * 3.0 + 7.0) + other = _model(shifted, observation_noise_var=model.observation_noise_var_) + width = lambda m, frame: (lambda q: q["q0.9"] - q["q0.1"])(m.predict_quantiles(frame, (0.1, 0.9))).to_numpy() + # the target scale is 3 times larger, the standardised variance is the same: the widths differ by that factor only + np.testing.assert_allclose(width(other, shifted) / width(model, df), 3.0, rtol=1e-6) + + +def test_a_longer_horizon_gives_wider_intervals(): + df = _frame(n=1500) + one = _model(df, horizon=1).predict_quantiles(df, (0.1, 0.9)) + four = _model(df, horizon=4).predict_quantiles(df, (0.1, 0.9)) + width_one = (one["q0.9"] - one["q0.1"]).to_numpy()[CALIBRATION:] + width_four = (four["q0.9"] - four["q0.1"]).to_numpy()[CALIBRATION:] + assert np.isfinite(width_four).all() and (width_four >= width_one - 1e-12).all() + + +def test_the_groups_are_scored_separately(): + df = _frame(n=800, groups=2) + quantiles = _model(df, groups=2).predict_quantiles(df, (0.1, 0.9)) + assert len(quantiles) == len(df) + assert np.isfinite(quantiles.to_numpy()).all() + + +def test_asking_for_the_variance_does_not_change_the_filter(): + df = _frame(n=1000) + model = _model(df) + prepared = model.prepare_data(df) + plain = model._run_filter(prepared, model.P_init_diag_, model.observation_noise_var_, model.process_noise_var_).clone() + with_variance = model._run_filter(prepared, model.P_init_diag_, model.observation_noise_var_, model.process_noise_var_, + with_variance=True) + assert torch.equal(plain, with_variance) diff --git a/tests/test_opera_pinball.py b/tests/test_opera_pinball.py new file mode 100644 index 0000000..5bc52e5 --- /dev/null +++ b/tests/test_opera_pinball.py @@ -0,0 +1,61 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +r""" +Pinball loss for ``OperaTAM`` (``loss_type='pinball', tau=...``) and ``aggregate_quantiles`` (one aggregation per level). + +Aggregating quantile forecasts needs the loss of each level: the square and absolute losses aggregate means and medians. +""" +import numpy as np +import pandas as pd +import pytest + +import tam as ta +from tam.model.opera import aggregate_quantiles + + +def _two_biased_experts(n: int = 3000, seed: int = 0) -> pd.DataFrame: + rng = np.random.default_rng(seed) + y = rng.standard_normal(n) + return pd.DataFrame({"t": np.arange(n), "y": y, "high": np.full(n, 1.28), "low": np.full(n, -1.28)}) + + +def _final_weights(frame: pd.DataFrame, **kwargs) -> np.ndarray: + model = ta.OperaTAM(target_col="y", expert_cols=["high", "low"], date_col="t", **kwargs).fit(frame) + return next(iter(model.weights_history_.values()))[-1] + + +def test_pinball_at_the_median_gives_the_absolute_loss_weights(): + # The pinball subgradient at tau = 0.5 is half the absolute one; MLpol's adaptive rates cancel a constant factor, up to + # the fixed initial learning rate (2**-20), which does not rescale with the loss: equal to about 1e-4. + frame = _two_biased_experts() + frame["high"] = frame["y"] + 0.3 + 0.5 * np.random.default_rng(1).standard_normal(len(frame)) + absolute = _final_weights(frame, algorithm="MLPOL", loss_type="absolute") + pinball = _final_weights(frame, algorithm="MLPOL", loss_type="pinball", tau=0.5) + np.testing.assert_allclose(pinball, absolute, atol=1e-3) + + +@pytest.mark.parametrize("algorithm", ["MLPOL", "EWA"]) +def test_high_levels_favour_the_high_expert(algorithm): + # y ~ N(0, 1): the 0.9 quantile is +1.28 (the "high" expert), the 0.1 quantile is -1.28 (the "low" expert). + frame = _two_biased_experts() + upper = _final_weights(frame, algorithm=algorithm, loss_type="pinball", tau=0.9, eta=0.5) + lower = _final_weights(frame, algorithm=algorithm, loss_type="pinball", tau=0.1, eta=0.5) + assert upper[0] > 0.8 and lower[1] > 0.8 + + +def test_aggregate_quantiles_columns_and_order(): + frame = _two_biased_experts(500) + experts = {0.1: ["low", "high"], 0.5: ["low", "high"], 0.9: ["low", "high"]} + table = aggregate_quantiles(frame, target_col="y", expert_cols_by_level=experts, date_col="t") + assert list(table.columns) == ["q0.1", "q0.5", "q0.9"] + assert len(table) == len(frame) + assert (np.diff(table.to_numpy(), axis=1) >= 0).all() + + +def test_invalid_pinball_level_and_loss_are_rejected(): + with pytest.raises(ValueError, match="tau"): + ta.OperaTAM(target_col="y", expert_cols=["a"], loss_type="pinball", tau=1.0) + with pytest.raises(ValueError, match="loss_type"): + ta.OperaTAM(target_col="y", expert_cols=["a"], loss_type="huber") diff --git a/tests/test_safety.py b/tests/test_safety.py index 0c22653..3b4e831 100644 --- a/tests/test_safety.py +++ b/tests/test_safety.py @@ -88,3 +88,59 @@ def test_normalized_intervals_scale_with_sigma(): intervals = model.predict_intervals(np.zeros(5), scale=test_scale) ratio = intervals["Width"].to_numpy() / test_scale assert np.allclose(ratio, ratio[0]) + + + +# --- verbose and multi-level quantiles ------------------------------------------------------------------------------ +_LEVELS = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] + + +def test_calibrate_prints_by_default_and_verbose_false_is_silent(capsys): + y = np.random.default_rng(0).standard_normal(100) + SafetyTAM(alpha=0.1).calibrate(y, np.zeros(100)) + assert "calibrated on 100 samples" in capsys.readouterr().out + SafetyTAM(alpha=0.1).calibrate(y, np.zeros(100), verbose=False) + assert capsys.readouterr().out == "" + + +def test_symmetric_quantiles_match_the_intervals(): + rng = np.random.default_rng(1) + y_pred, scale = rng.normal(size=50), rng.uniform(0.5, 2.0, 50) + engine = SafetyTAM(alpha=0.2).calibrate(rng.standard_normal(500), np.zeros(500), verbose=False) + table = engine.predict_quantiles(y_pred, _LEVELS, scale=scale) + assert list(table.columns) == [f"q{t}" for t in _LEVELS] + band = engine.predict_intervals(y_pred, scale=scale) # alpha 0.2 -> levels 0.1 and 0.9 + np.testing.assert_allclose(table["q0.1"], band["Lower"]) + np.testing.assert_allclose(table["q0.9"], band["Upper"]) + np.testing.assert_allclose(table["q0.5"], y_pred) + + +def test_quantiles_are_non_decreasing_across_levels(): + rng = np.random.default_rng(2) + engine = SafetyTAM().calibrate(rng.lognormal(size=300), np.ones(300), verbose=False) + for signed in (False, True): + table = engine.predict_quantiles(rng.normal(size=40), _LEVELS, signed=signed).to_numpy() + assert (np.diff(table, axis=1) >= 0).all() + + +def test_signed_quantiles_cover_each_tail_of_skewed_errors(): + # Lognormal errors are skewed: the symmetric (absolute-score) band is too wide on the short side and too narrow on + # the long one; signed scores calibrate each tail on its own. + rng = np.random.default_rng(3) + errors = lambda n: rng.lognormal(0.0, 0.8, n) - np.exp(0.32) # noqa: E731 - centred at the mean + engine = SafetyTAM().calibrate(errors(4000), np.zeros(4000), verbose=False) + future = errors(20000) + for signed, tolerance in ((True, 0.02), (False, None)): + table = engine.predict_quantiles(np.zeros_like(future), [0.1, 0.9], signed=signed) + below_lower = float(np.mean(future < table["q0.1"])) + above_upper = float(np.mean(future > table["q0.9"])) + if signed: + assert abs(below_lower - 0.1) < tolerance and abs(above_upper - 0.1) < tolerance + else: + assert below_lower < 0.08 and above_upper > 0.12 # the symmetric band misplaces both tails (about 6.5 % and 13 %) + + +def test_signed_quantiles_need_signed_scores(): + engine = SafetyTAM().calibrate_scores(np.abs(np.random.default_rng(4).standard_normal(100))) + with pytest.raises(ValueError, match="signed"): + engine.predict_quantiles(np.zeros(3), _LEVELS, signed=True) diff --git a/tests/test_target_transform.py b/tests/test_target_transform.py new file mode 100644 index 0000000..39bda0f --- /dev/null +++ b/tests/test_target_transform.py @@ -0,0 +1,165 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +r""" +Tests for ``dist_kwargs={"target_transform": ...}`` of the distributional :class:`StaticTAM`: ``"asinh"`` for a skewed real-valued target +(negative values allowed) and ``"logit"`` for a target bounded in an interval. ``log_target`` keeps its meaning and its default. +""" +import numpy as np +import pandas as pd +import pytest +from scipy import stats + +from tam import StaticTAM + +_FORMULAS = {"mu": "y ~ s(x)", "sigma": "y ~ s(x)"} +_LOSS = {"mu": "l2", "sigma": "gamma"} +_TAUS = (0.05, 0.5, 0.95) + + +def _model(**dist_kwargs): + return StaticTAM(_FORMULAS, loss=_LOSS, dist_kwargs=dist_kwargs) + + +def _johnson_su(rng, n): + """A skewed, heavy-tailed, real-valued target with negatives: y = -1 + 4 sinh(mu(x) + 1.2 z), the asinh-normal (Johnson SU) family.""" + x = rng.uniform(0.0, 1.0, n) + y = -1.0 + 4.0 * np.sinh(1.5 * (x - 0.5) + 1.2 * rng.standard_normal(n)) + return pd.DataFrame({"x": x, "y": y}) + + +def _bounded(rng, n): + x = rng.uniform(0.0, 1.0, n) + mean = 0.15 + 0.7 * x + y = rng.beta(mean * 8.0, (1.0 - mean) * 8.0) + return pd.DataFrame({"x": x, "y": y}) + + +def _coverage(frame, quantiles, tau): + return float(np.mean(frame["y"].to_numpy() <= quantiles[f"q{tau}"].to_numpy())) + + +def test_asinh_recovers_the_tails_of_a_skewed_target_with_negatives(): + rng = np.random.default_rng(0) + train, test = _johnson_su(rng, 6000), _johnson_su(rng, 6000) + assert (train["y"] < 0).any() + asinh = _model(target_transform="asinh", tail_family="normal").fit(train) + plain = _model(log_target=False, tail_family="normal").fit(train) + quantiles = asinh.predict_quantiles(test, _TAUS) + assert abs(_coverage(test, quantiles, 0.05) - 0.05) < 0.02 + assert abs(_coverage(test, quantiles, 0.95) - 0.95) < 0.02 + # a symmetric law on the response scale misplaces the tails of a skewed target by more than two points + plain_quantiles = plain.predict_quantiles(test, _TAUS) + plain_errors = [abs(_coverage(test, plain_quantiles, tau) - tau) for tau in (0.05, 0.95)] + assert max(plain_errors) > 0.02 + + +def test_asinh_quantiles_do_not_cross_and_the_median_is_the_back_transformed_location(): + rng = np.random.default_rng(1) + train = _johnson_su(rng, 3000) + model = _model(target_transform="asinh").fit(train) + quantiles = model.predict_quantiles(train, (0.1, 0.5, 0.9)) + assert (quantiles["q0.1"] <= quantiles["q0.5"]).all() and (quantiles["q0.5"] <= quantiles["q0.9"]).all() + np.testing.assert_allclose(model.predict_median(train), quantiles["q0.5"].to_numpy(), rtol=1e-8, atol=1e-8) + + +def test_asinh_cdf_and_crps_and_anomaly_score_work_through_the_transform(): + rng = np.random.default_rng(2) + train = _johnson_su(rng, 3000) + model = _model(target_transform="asinh", tail_family="normal").fit(train) + pit = model.cdf(train) + assert 0.0 <= pit.min() and pit.max() <= 1.0 + assert abs(float(np.mean(pit)) - 0.5) < 0.03 + assert np.isfinite(model.crps(train)).all() and (model.crps(train) >= 0).all() + score = model.anomaly_score(train) + np.testing.assert_allclose(score["predicted_median"].to_numpy(), model.predict_median(train)) + + +def test_logit_quantiles_stay_inside_the_bounds_and_cover(): + rng = np.random.default_rng(3) + train, test = _bounded(rng, 5000), _bounded(rng, 5000) + model = _model(target_transform="logit", bounds=(0.0, 1.0), tail_family="normal").fit(train) + quantiles = model.predict_quantiles(test, (0.01, 0.05, 0.5, 0.95, 0.99)) + assert quantiles.to_numpy().min() > 0.0 and quantiles.to_numpy().max() < 1.0 + assert abs(_coverage(test, quantiles, 0.05) - 0.05) < 0.03 + assert abs(_coverage(test, quantiles, 0.95) - 0.95) < 0.03 + # a normal law on the response scale puts quantiles outside the bounds + plain = _model(log_target=False, tail_family="normal").fit(train).predict_quantiles(test, (0.01, 0.99)) + assert plain.to_numpy().min() < 0.0 or plain.to_numpy().max() > 1.0 + + +def test_logit_estimates_and_stores_bounds_when_none_are_given(): + rng = np.random.default_rng(4) + train = _bounded(rng, 3000) + model = _model(target_transform="logit", tail_family="normal").fit(train) + low, high = model._transform_params_["bounds"] + assert low < train["y"].min() and high > train["y"].max() + quantiles = model.predict_quantiles(train, (0.001, 0.999)) + assert quantiles.to_numpy().min() > low and quantiles.to_numpy().max() < high + + +def test_logit_values_at_the_bounds_are_nudged_not_rejected(): + rng = np.random.default_rng(5) + train = _bounded(rng, 2000) + train.loc[:20, "y"] = 0.0 + train.loc[21:40, "y"] = 1.0 + model = _model(target_transform="logit", bounds=(0.0, 1.0)).fit(train) + assert np.isfinite(model.predict_median(train)).all() + + +def test_logit_rejects_values_outside_the_given_bounds(): + rng = np.random.default_rng(6) + train = _bounded(rng, 500) + train.loc[0, "y"] = 1.5 + with pytest.raises(ValueError, match=r"outside the bounds"): + _model(target_transform="logit", bounds=(0.0, 1.0)).fit(train) + + +def test_an_unknown_transform_or_a_contradiction_with_log_target_is_rejected(): + with pytest.raises(ValueError, match="target_transform"): + _model(target_transform="sqrt") + with pytest.raises(ValueError, match="log_target"): + _model(target_transform="asinh", log_target=True) + with pytest.raises(ValueError, match="bounds"): + _model(target_transform="logit", bounds=(1.0, 0.0)) + + +def test_a_mixture_does_not_accept_a_target_transform(): + with pytest.raises(ValueError, match="mixture"): + StaticTAM("y ~ s(x)", mixture_components=2, dist_kwargs={"target_transform": "asinh"}) + + +def test_the_default_and_log_target_true_are_the_same_log_fit(): + rng = np.random.default_rng(7) + x = rng.uniform(0.0, 1.0, 2000) + train = pd.DataFrame({"x": x, "y": np.exp(1.0 + x + 0.3 * rng.standard_normal(2000))}) + default = _model().fit(train).predict_quantiles(train, _TAUS) + explicit = _model(log_target=True).fit(train).predict_quantiles(train, _TAUS) + named = _model(target_transform="log").fit(train).predict_quantiles(train, _TAUS) + pd.testing.assert_frame_equal(default, explicit) + pd.testing.assert_frame_equal(default, named) + assert default.to_numpy().min() > 0.0 + + +def test_the_conformal_intervals_cover_a_target_modelled_through_asinh(): + rng = np.random.default_rng(8) + train, calibration, test = (_johnson_su(rng, 3000) for _ in range(3)) + model = _model(target_transform="asinh", tail_family="normal").fit(train) + model.calibrate_conformal(calibration, alpha=0.1, studentized=True) + intervals = model.predict_intervals(test) + inside = (test["y"].to_numpy() >= intervals["Lower"].to_numpy()) & (test["y"].to_numpy() <= intervals["Upper"].to_numpy()) + assert abs(float(inside.mean()) - 0.9) < 0.03 + + +def test_the_back_transform_inverts_the_transform(): + rng = np.random.default_rng(9) + train = _johnson_su(rng, 2000) + for kwargs in ({"target_transform": "asinh"}, {"log_target": False}): + model = _model(**kwargs).fit(train) + back = model._from_model_scale(model._to_model_scale(train["y"].to_numpy()), "round trip") + np.testing.assert_allclose(back, train["y"].to_numpy(), rtol=1e-8, atol=1e-8) + bounded = _bounded(rng, 2000) + logit = _model(target_transform="logit", bounds=(0.0, 1.0)).fit(bounded) + back = logit._from_model_scale(logit._to_model_scale(bounded["y"].to_numpy()), "round trip") + np.testing.assert_allclose(back, bounded["y"].to_numpy(), rtol=1e-8, atol=1e-8) diff --git a/tests/test_thread_limits.py b/tests/test_thread_limits.py new file mode 100644 index 0000000..0860df8 --- /dev/null +++ b/tests/test_thread_limits.py @@ -0,0 +1,58 @@ +# SPDX-FileCopyrightText: 2026 EDF (Electricité De France) +# SPDX-FileContributor: Yann Allioux +# SPDX-License-Identifier: LGPL-3.0-or-later + +r""" +Limiting the CPU threads. The supported way is the environment (``OMP_NUM_THREADS`` and ``MKL_NUM_THREADS``, set before Python starts). + +Some PyTorch CPU builds (seen with torch 2.14.0+cpu on Windows, MKL 2026.1) make any batched LU factorisation of a matrix of a few hundred +rows (``torch.linalg.solve``, ``lu_factor``, ``inv``) print ``Intel oneMKL ERROR: Parameter 6 was incorrect on entry to DLASWP`` and +hang once ``torch.set_num_threads(n)`` has been called with ``n >= 2``, whatever ``n`` and without importing ``tam``. The test runs a +grouped fit large enough to need such a factorisation in a fresh process with the threads limited through the environment: it finishes, +uses the requested number of threads, prints no MKL error and takes about as long as without the limit. +""" +import os +import subprocess +import sys +import time + +import pytest + +_SCRIPT = r""" +import time, warnings +warnings.filterwarnings("ignore") +import numpy as np, pandas as pd, torch +import tam as ta +rng = np.random.default_rng(0) +days, groups = 400, 24 +n = days * groups +df = pd.DataFrame({"date": np.repeat(pd.date_range("2020-01-01", periods=days, freq="D"), groups), "h": np.tile(np.arange(groups), days), + "x": rng.normal(size=n), "z": rng.normal(size=n)}) +df["y"] = np.sin(df.x) + 0.5 * df.z + 0.1 * rng.normal(size=n) +start = time.time() +ta.StaticTAM("y ~ s(x, k=30) + s(z, k=30) + te(s(x, k=10), s(z, k=10))", group_col="h", date_col="date").fit(df) +print("THREADS", torch.get_num_threads(), "SECONDS", round(time.time() - start, 3)) +""" + + +def _run(threads=None, timeout=240): + env = {key: value for key, value in os.environ.items() if key not in ("OMP_NUM_THREADS", "MKL_NUM_THREADS")} + if threads is not None: + env.update(OMP_NUM_THREADS=str(threads), MKL_NUM_THREADS=str(threads)) + started = time.time() + done = subprocess.run([sys.executable, "-c", _SCRIPT], capture_output=True, text=True, timeout=timeout, env=env) + return done, time.time() - started + + +def _field(done, name): + return float(done.stdout.split(name)[1].split()[0]) + + +def test_the_environment_limits_the_threads_and_the_fit_finishes_without_mkl_errors(): + free, free_seconds = _run() + assert free.returncode == 0, free.stderr[-500:] + limited, limited_seconds = _run(threads=2) + assert limited.returncode == 0, limited.stderr[-500:] + assert _field(limited, "THREADS") == 2 + assert "DLASWP" not in limited.stdout + limited.stderr + assert limited_seconds < 10 * max(free_seconds, 5.0)