Quantile regression via the Extended Log-F (ELF) pseudo-family.
QuantileFamily(
tau=0.5,
sigma=1.0,
)
QuantileFamily fits a single conditional quantile 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 or 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 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 per tau and compare/combine the resulting models, or see 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 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:
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
|
Whether the dispersion (scale) parameter is fixed rather than estimated.
|
|
sigma
|
Bandwidth parameter sigma of the ELF loss approximation.
|
|
tau
|
Quantile level tau being estimated.
|
scale_known
Whether the dispersion (scale) parameter is fixed rather than estimated.
Always True for QuantileFamily, since the ELF loss has no free scale parameter to estimate.
sigma
Bandwidth parameter sigma of the ELF loss approximation.
tau
Quantile level tau being estimated.
Methods
|
Name
|
Description
|
|
deviance()
|
Total deviance, twice the summed ELF loss over observations.
|
|
initialize()
|
Initialize the mean mu from the observed response y.
|
|
irls_update()
|
Compute the working response and weight for one custom IRLS step.
|
|
link()
|
Identity link: map the mean mu directly to the linear predictor eta.
|
|
link_derivative()
|
Derivative of the identity link with respect to mu.
|
|
link_inverse()
|
Identity inverse link: map the linear predictor eta back to mu.
|
|
log_likelihood()
|
Pseudo log-likelihood implied by the ELF loss.
|
|
simulate()
|
Draw random samples from the Asymmetric Laplace distribution.
|
|
unit_deviance()
|
Per-observation deviance contribution, twice the ELF loss.
|
|
variance()
|
Variance function for QuantileFamily, identically equal to 1.
|
deviance()
Total deviance, twice the summed ELF loss over observations.
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.
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.
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.
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.
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.
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.
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.
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 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.
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 — 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.
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()
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.