# QuantileFamily Quantile regression via the Extended Log-F (ELF) pseudo-family. Usage ``` python QuantileFamily( tau=0.5, sigma=1.0, ) ``` [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily) fits a single conditional quantile [tau](QuantileFamily.md#whittaker.QuantileFamily.tau) of the response -- for example the median (`tau=0.5`) or the 90th percentile (`tau=0.9`) -- rather than the conditional mean targeted by families such as [Gaussian](Gaussian.md#whittaker.Gaussian) or [Gamma](Gamma.md#whittaker.Gamma). This is useful whenever the object of interest is not the average behavior of the response but a specific point in its distribution, e.g. modeling the upper tail of a skewed cost distribution, or building prediction bands by fitting several quantiles ([tau](QuantileFamily.md#whittaker.QuantileFamily.tau) values) side by side. Rather than the non-smooth pinball ("check") loss used by classical quantile regression, Whittaker uses the smooth Extended Log-F (ELF) approximation of Fasiolo et al. (2021), which is differentiable and therefore fits within the standard P-IRLS loop via a custom `irls_update`. To fit several quantiles jointly with a shared smoothness structure, fit one [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily) per [tau](QuantileFamily.md#whittaker.QuantileFamily.tau) and compare/combine the resulting models, or see [ConformalPredictor](ConformalPredictor.md#whittaker.ConformalPredictor) for distribution-free coverage guarantees around a fitted mean model. ## Parameters `tau: float = ``0.5` Quantile level to estimate, in `(0, 1)`. `tau=0.5` corresponds to median regression; smaller values target lower quantiles and larger values target upper quantiles. `sigma: float = ``1.0` Bandwidth (smoothing) parameter controlling how closely the ELF loss approximates the non-smooth pinball loss. Smaller values give a sharper, more faithful approximation to the pinball loss (and to the check-function optimum) but a less smooth optimization surface; larger values give a smoother but more biased approximation. `sigma` may be adjusted after construction via the `sigma` property, e.g. to anneal it across fitting iterations. ## Notes [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily) has no meaningful link or variance function in the usual GLM sense -- `link` and `link_inverse` are the identity on `eta` -- because fitting instead minimizes the ELF loss directly. For a residual `r = y - mu`, the ELF loss is \rho\_{\tau,\sigma}(r) = \tau r + \sigma \log\\\left(1 + e^{-r/\sigma}\right), which converges to the pinball loss \rho\_\tau(r) = \tau r \\ \mathbb{1}\[r \ge 0\] - (1-\tau) r \\ \mathbb{1}\[r \< 0\] as \sigma \to 0. Its first and second derivatives with respect to `mu`, \frac{\partial \rho}{\partial \mu} = -\left\[\tau - 1 + \operatorname{expit}(r/\sigma)\right\], \qquad \frac{\partial^2 \rho}{\partial \mu^2} = \frac{1}{\sigma}\\ s (1 - s), \quad s = \operatorname{expit}(r/\sigma), supply the working response `z` and working weight `W` used by the custom `irls_update`. The reported "deviance" is `2 * sum(ELF loss)`, so that it reduces to twice the usual pinball loss in the limit `sigma -> 0`. ## Examples Fit the 10th, 50th, and 90th percentile curves of a heteroscedastic response: ``` python import numpy as np import whittaker as wk rng = np.random.default_rng(0) n = 300 x = np.linspace(0, 2 * np.pi, n) mu = np.sin(x) noise_scale = 0.2 + 0.3 * np.abs(np.cos(x)) y = mu + rng.normal(0, noise_scale, n) data = {"x": x, "y": y} for tau in (0.1, 0.5, 0.9): model = wk.GAM("y ~ s(x)", family=wk.QuantileFamily(tau=tau)) model.fit(data, method="REML") print(f"tau={tau}:") print(model.summary()) ``` tau=0.1: GAM fit summary ============================================================ Formula: y ~ s(x) Family: Quantile(tau=0.1, sigma=1.0) Inference: REML Observations: 300 Coefficients: 10 Parametric coefficients: Term Estimate Std.Err z value p-value ------------------------ ---------- ---------- ---------- ---------- (Intercept) -2.3565 0.1956 -12.045 < 1e-16 Approximate significance of smooth terms: Term EDF Ref.df Chi.sq p-value ------------------------ ------ ------ ---------- ---------- s(x) 1.00 2 8.706 0.01287 Total EDF: 2.00 Scale est: 1.000000 Deviance: 204.5488 Null dev: 465.6825 Dev. expl: 56.1% GCV score: 0.691012 AIC: 208.55 BIC: 215.96 tau=0.5: GAM fit summary ============================================================ Formula: y ~ s(x) Family: Quantile(tau=0.5, sigma=1.0) Inference: REML Observations: 300 Coefficients: 10 Parametric coefficients: Term Estimate Std.Err z value p-value ------------------------ ---------- ---------- ---------- ---------- (Intercept) -0.0196 0.1177 -0.166 0.8678 Approximate significance of smooth terms: Term EDF Ref.df Chi.sq p-value ------------------------ ------ ------ ---------- ---------- s(x) 3.66 4 36.020 2.866e-07 Total EDF: 4.66 Scale est: 1.000000 Deviance: 428.0705 Null dev: 465.6825 Dev. expl: 8.1% GCV score: 1.472271 AIC: 437.39 BIC: 454.64 tau=0.9: GAM fit summary ============================================================ Formula: y ~ s(x) Family: Quantile(tau=0.9, sigma=1.0) Inference: REML Observations: 300 Coefficients: 10 Parametric coefficients: Term Estimate Std.Err z value p-value ------------------------ ---------- ---------- ---------- ---------- (Intercept) 2.3146 0.1955 11.842 < 1e-16 Approximate significance of smooth terms: Term EDF Ref.df Chi.sq p-value ------------------------ ------ ------ ---------- ---------- s(x) 1.00 2 8.887 0.01175 Total EDF: 2.00 Scale est: 1.000000 Deviance: 204.3576 Null dev: 465.6825 Dev. expl: 56.1% GCV score: 0.690366 AIC: 208.36 BIC: 215.77 ## Attributes | Name | Description | |----|----| | [scale_known](#scale_known) | Whether the dispersion (scale) parameter is fixed rather than estimated. | | [sigma](#sigma) | Bandwidth parameter `sigma` of the ELF loss approximation. | | [tau](#tau) | Quantile level [tau](QuantileFamily.md#whittaker.QuantileFamily.tau) being estimated. | ------------------------------------------------------------------------ ### scale_known Whether the dispersion (scale) parameter is fixed rather than estimated. `scale_known: bool` Always True for [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily), since the ELF loss has no free scale parameter to estimate. ------------------------------------------------------------------------ ### sigma Bandwidth parameter `sigma` of the ELF loss approximation. `sigma: float` ------------------------------------------------------------------------ ### tau Quantile level [tau](QuantileFamily.md#whittaker.QuantileFamily.tau) being estimated. `tau: float` ## Methods | Name | Description | |----|----| | [deviance()](#deviance) | Total deviance, twice the summed ELF loss over observations. | | [initialize()](#initialize) | Initialize the mean `mu` from the observed response `y`. | | [irls_update()](#irls_update) | Compute the working response and weight for one custom IRLS step. | | [link()](#link) | Identity link: map the mean `mu` directly to the linear predictor `eta`. | | [link_derivative()](#link_derivative) | Derivative of the identity link with respect to `mu`. | | [link_inverse()](#link_inverse) | Identity inverse link: map the linear predictor `eta` back to `mu`. | | [log_likelihood()](#log_likelihood) | Pseudo log-likelihood implied by the ELF loss. | | [simulate()](#simulate) | Draw random samples from the Asymmetric Laplace distribution. | | [unit_deviance()](#unit_deviance) | Per-observation deviance contribution, twice the ELF loss. | | [variance()](#variance) | Variance function for [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily), identically equal to 1. | ------------------------------------------------------------------------ ### deviance() Total deviance, twice the summed ELF loss over observations. Usage ``` python deviance( y, mu, *, weights=None, ) ``` Implements the family-specific deviance as D = 2 \sum_i w_i\\ \rho\_{\tau,\sigma}(y_i - \mu_i), where `w_i` are optional observation weights (defaulting to 1) and `rho` is the ELF loss. As `sigma -> 0` this converges to twice the usual pinball loss, matching the deviance convention used by other families. #### Parameters `y: NDArray` Observed response values. `mu: NDArray` Fitted mean values. `weights: NDArray or None = None` Optional observation weights. If None, all weights are treated as 1. #### Returns `float` The total (weighted) deviance. ------------------------------------------------------------------------ ### initialize() Initialize the mean `mu` from the observed response `y`. Usage ``` python initialize(y) ``` Implements the family-specific starting values for P-IRLS as a direct copy of the observed response, `mu = y`, which is a reasonable starting point for quantile fitting regardless of [tau](QuantileFamily.md#whittaker.QuantileFamily.tau). #### Parameters `y: NDArray` Observed response values. #### Returns `NDArray` Copy of `y`, used as the initial mean estimate. ------------------------------------------------------------------------ ### irls_update() Compute the working response and weight for one custom IRLS step. Usage ``` python irls_update( y, mu, eta, ) ``` Implements the family-specific IRLS update used in place of the standard GLM working response/weight pair, since the ELF loss does not come from an exponential-family model. For residual `r = y - mu`, the update uses the ELF score and curvature u = \tau - 1 + \operatorname{expit}(r/\sigma), \qquad W = \max\\\left(\frac{1}{\sigma}\\ s(1-s),\\ \epsilon\right), \quad s = \operatorname{expit}(r/\sigma), where `epsilon` is a small floor that prevents zero weights, and forms the working response `z = eta + u / W`. #### Parameters `y: NDArray` Observed response values. `mu: NDArray` Current fitted mean values. `eta: NDArray` Current linear predictor values. #### Returns `tuple[NDArray, NDArray]` The working response `z` and working weight `W`, both with the same shape as `y`. ------------------------------------------------------------------------ ### link() Identity link: map the mean `mu` directly to the linear predictor `eta`. Usage ``` python link(mu) ``` [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily) does not use a link function in the usual GLM sense because it minimizes the ELF loss directly rather than modeling a mean-variance relationship, so this implements the family-specific link as the identity, `eta = mu`. #### Parameters `mu: NDArray` Mean values on the response scale. #### Returns `NDArray` Same values, unchanged, interpreted as the linear predictor `eta`. ------------------------------------------------------------------------ ### link_derivative() Derivative of the identity link with respect to `mu`. Usage ``` python link_derivative(mu) ``` Since `link` is the identity, `d(eta)/d(mu) = 1` everywhere, so this implements the family-specific derivative as an array of ones with the same shape as `mu`. #### Parameters `mu: NDArray` Mean values on the response scale. #### Returns `NDArray` Array of ones, same shape as `mu`. ------------------------------------------------------------------------ ### link_inverse() Identity inverse link: map the linear predictor `eta` back to `mu`. Usage ``` python link_inverse(eta) ``` Implements the family-specific inverse link as the identity, `mu = eta`, matching `link` above. #### Parameters `eta: NDArray` Linear predictor values. #### Returns `NDArray` Same values, unchanged, interpreted as the mean `mu`. ------------------------------------------------------------------------ ### log_likelihood() Pseudo log-likelihood implied by the ELF loss. Usage ``` python log_likelihood( y, mu, scale, *, weights=None, ) ``` Implements the family-specific log-likelihood as the negative summed ELF loss, \ell = -\sum_i w_i\\ \rho\_{\tau,\sigma}(y_i - \mu_i), with optional observation weights `w_i`. This is a pseudo-likelihood used for model comparison (e.g. AIC) rather than a true likelihood, since [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily) does not correspond to a proper probability model; the `scale` argument is accepted for interface compatibility but unused because `scale_known` is always True. #### Parameters `y: NDArray` Observed response values. `mu: NDArray` Fitted mean values. `scale: float` Dispersion parameter; unused since the ELF loss has no free scale. `weights: NDArray or None = None` Optional observation weights. If None, all weights are treated as 1. #### Returns `float` The (weighted) pseudo log-likelihood. ------------------------------------------------------------------------ ### simulate() Draw random samples from the Asymmetric Laplace distribution. Usage ``` python simulate( mu, scale, rng, ) ``` Implements the family-specific sampler by inverting the CDF of the Asymmetric Laplace distribution with location `mu`, scale `sigma`, and asymmetry [tau](QuantileFamily.md#whittaker.QuantileFamily.tau) -- the distribution whose negative log-density is proportional to the pinball loss that the ELF loss approximates. Draws `u ~ Uniform(0, 1)` and returns y = \begin{cases} \mu + \sigma \log(u / \tau), & u \< \tau, \\ \mu - \sigma \log\\\left(\dfrac{1-u}{1-\tau}\right), & u \ge \tau. \end{cases} #### Parameters `mu: NDArray` Location (fitted mean) values. `scale: float` Dispersion parameter; unused since sampling instead uses `sigma`. `rng: object` Random number generator exposing a `uniform(size=...)` method. #### Returns `NDArray` Simulated response values, same shape as `mu`. ------------------------------------------------------------------------ ### unit_deviance() Per-observation deviance contribution, twice the ELF loss. Usage ``` python unit_deviance( y, mu, ) ``` Implements the family-specific unit deviance as d_i = 2 \rho\_{\tau,\sigma}(y_i - \mu_i), the unweighted, un-summed term whose sum (optionally weighted) gives `deviance`. #### Parameters `y: NDArray` Observed response values. `mu: NDArray` Fitted mean values. #### Returns `NDArray` Per-observation deviance contributions, same shape as `y`. ------------------------------------------------------------------------ ### variance() Variance function for [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily), identically equal to 1. Usage ``` python variance(mu) ``` [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily) does not model a mean-variance relationship since fitting minimizes the ELF loss directly, so the family-specific variance function returns a constant array of ones regardless of `mu`; the working weights used during fitting instead come from `irls_update`. #### Parameters `mu: NDArray` Mean values on the response scale (unused beyond determining shape). #### Returns `NDArray` Array of ones, same shape as `mu`.