Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Quantile Regression › Implementation in NumPy, SciPy, and scikit-learn

Michael BrenndoerferOctober 26, 202610 min read

Machine Learning from ScratchAll chapters

Part 5 of 15Linear and Generalized Linear Models

Chapter 47Quantile Regression › Implementation in NumPy, SciPy, and scikit-learn

Reading section 6 of 10 · Implementation in NumPy, SciPy, and scikit-learn · 10 minSection 6 of 10 · 10 min

Section 6 of 10 · 10 min read

Quantile Regression › Implementation in NumPy, SciPy, and scikit-learnSection view · Sections in this chapter6 / 10 ↓

Sections in this chapter10

Reading tools Definitions, links & notes

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Notes & comments

Notes & commentsLeave notes and join the discussionSelection toolbarShow comment actions when selecting text. Turn off for the browser menu only.Comment indicatorsShow small counts beside annotated passagesPassage highlightsGently shade annotated text selections

We build the linear program from first principles with NumPy and SciPy, because the point of this chapter is to see the mechanism rather than to call it. Only after our own solver reproduces the reference fit do we bring in scikit-learn's QuantileRegressor, which solves an equivalent unpenalized problem with the HiGHS backend when configured with the settings used below.

We start with a synthetic dataset that is deliberately heteroscedastic: the response is linear in $x$ , and the noise scales up with $x$ . This is a best case for quantile regression, since the conditional quantiles really are linear, so any coverage shortfall we observe later comes from estimation noise rather than misspecification.

PythonIn [11]
CollapseExpandrng = default_rng(20240517)

n_total = 900
x_all = rng.uniform(0.0, 10.0, size=n_total)
sigma_all = 0.6 + 0.45 * x_all  # spread grows with x
y_all = 2.0 + 1.1 * x_all + sigma_all * rng.standard_normal(n_total)

perm = rng.permutation(n_total)
train_rows, test_rows = perm[:600], perm[600:]


def add_intercept(x):
    return np.column_stack([np.ones_like(x), x])


X_tr, y_tr = add_intercept(x_all[train_rows]), y_all[train_rows]
X_te, y_te = add_intercept(x_all[test_rows]), y_all[test_rows]

ConsoleOut [12]
CollapseExpandtraining rows: 600, test rows: 300, parameters: 2
train noise scale ranges from 0.60 to 5.07

Two functions carry the whole method: the pinball loss for scoring, and the linear program for fitting. We write the constraint matrix in the block form $\begin{bmatrix} X & I_n & -I_n\end{bmatrix}$ derived above, and let SciPy's HiGHS interface choose its supported LP algorithm.

PythonIn [13]
CollapseExpanddef pinball_loss(y, yhat, tau):
    """Average sample pinball loss for supplied predictions; its population counterpart targets the conditional tau-quantile."""
    resid = y - yhat
    return float(
        np.mean(np.where(resid >= 0.0, tau * resid, (tau - 1.0) * resid))
    )


def fit_quantile_lp(X, y, tau, return_result=False):
    """Solve min sum(tau*u + (1-tau)*v) s.t. X beta + u - v = y, u, v >= 0."""
    n, p = X.shape
    cost = np.concatenate(
        [np.zeros(p), tau * np.ones(n), (1.0 - tau) * np.ones(n)]
    )
    A_eq = np.hstack([X, np.eye(n), -np.eye(n)])
    bounds = [(None, None)] * p + [(0.0, None)] * (2 * n)
    result = linprog(cost, A_eq=A_eq, b_eq=y, bounds=bounds, method="highs")
    if not result.success:
        raise RuntimeError(f"linear program failed: {result.message}")
    beta = result.x[:p]
    return (beta, result) if return_result else beta

We fit seven quantile levels on the training data. Each call solves one linear program. Running more levels adds their individual solve costs; those solve costs can differ across levels.

PythonIn [14]
CollapseExpandtaus = np.array([0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95])
betas = {tau: fit_quantile_lp(X_tr, y_tr, tau) for tau in taus}

coef_table = pd.DataFrame(
    {
        "intercept": [betas[t][0] for t in taus],
        "slope": [betas[t][1] for t in taus],
    },
    index=pd.Index([f"{t:g}" for t in taus], name="tau"),
).round(3)
coef_table

Out[15]:

	intercept	slope
tau
0.05	1.126	0.362
0.1	1.106	0.570
0.25	1.504	0.831
0.5	2.074	1.058
0.75	2.485	1.385
0.9	2.443	1.794
0.95	2.815	1.906

Read the table as a description of a stretching distribution. The true data-generating quantiles are $Q_\tau(x)=2+1.1x+(0.6+0.45x)\Phi^{-1}(\tau)$ , so their slopes are $1.1+0.45\Phi^{-1}(\tau)$ . At levels 0.05 and 0.95 these are approximately 0.360 and 1.840. The fitted slopes are estimates and their intercepts need not be monotone; we check the resulting predictions for crossings rather than inferring order from their coefficient table. One fitted sample also does not estimate the sampling variability of its coefficients.

The coefficient table can be read more directly by comparing its slope column with the true slope function. The next figure does that at the fitted levels, using the already computed coefficient vectors.

PythonIn [16]
CollapseExpandtau_slope_grid = np.linspace(0.05, 0.95, 201)
true_slope_curve = 1.1 + 0.45 * norm.ppf(tau_slope_grid)
fitted_slopes = np.array([betas[t][1] for t in taus])
true_slopes_at_fits = 1.1 + 0.45 * norm.ppf(taus)

VisualizationOut [17]
CollapseExpand

Estimated slope points against tau lie close to a smooth increasing true slope curve. — Fitted quantile slopes at seven levels compared with the true slope function implied by the heteroscedastic Gaussian generator. The slopes rise with tau, and the fitted points lie close to the true curve, showing that separate level fits recover the growing conditional spread as the noise scales up with the predictor. Each level is estimated by its own linear program, yet the estimates follow the same underlying pattern.

The helper normally returns only coefficients. Setting return_result=True also returns SciPy's successful LP result so we can inspect its actual equality-constraint dual multipliers. We check the counting inequality with a stated numerical zero tolerance, the full stationarity condition, the dual bounds, and the primal-dual objective gap.

PythonIn [18]
CollapseExpandcheck_rows = []
for tau in (0.10, 0.50, 0.90):
    beta, lp_result = fit_quantile_lp(X_tr, y_tr, tau, return_result=True)
    resid = y_tr - X_tr @ beta
    zero_tol = 1e-7 * max(1.0, float(np.max(np.abs(y_tr))))
    n_negative = int(np.sum(resid < -zero_tol))
    n_zero = int(np.sum(np.abs(resid) <= zero_tol))
    dual = lp_result.eqlin.marginals
    stationarity = float(np.max(np.abs(X_tr.T @ dual)))
    dual_bounds_ok = bool(
        np.all(dual >= tau - 1 - 1e-8) and np.all(dual <= tau + 1e-8)
    )
    count_ok = bool(n_negative <= len(y_tr) * tau <= n_negative + n_zero)
    objective_gap = float(abs(lp_result.fun - y_tr @ dual) / len(y_tr))
    assert (
        count_ok
        and dual_bounds_ok
        and stationarity < 1e-7
        and objective_gap < 1e-7
    )
    check_rows.append(
        {
            "tau": tau,
            "negative residuals": n_negative,
            "near-zero residuals": n_zero,
            "count inequality": count_ok,
            "max stationarity residual": stationarity,
            "mean primal-dual gap": objective_gap,
        }
    )
pd.DataFrame(check_rows)

Out[19]:

	tau	negative residuals	near-zero residuals	count inequality	max stationarity residual	mean primal-dual gap
0	0.1	59	2	True	5.224987e-13	1.042129e-15
1	0.5	299	2	True	4.574119e-14	3.789561e-16
2	0.9	539	2	True	4.218847e-15	8.526513e-16

The zero tolerance classifies residuals near a fitted hyperplane; it is not a license to assign them an arbitrary gradient. The returned dual values satisfy the permitted intervals and balance all design columns. Primal feasibility from the LP, dual feasibility, and agreement of their objectives provide an optimality check to the stated solver precision, while leaving coefficient uniqueness as a separate question.

Now the reference implementation. scikit-learn's QuantileRegressor minimizes the same objective with an optional $\ell_1$ penalty on the coefficients; setting alpha=0 removes the penalty, and fit_intercept=False prevents it from adding a second intercept column to the one we already built.

PythonIn [20]
CollapseExpandreference = QuantileRegressor(
    quantile=0.5, alpha=0.0, fit_intercept=False, solver="highs"
)
reference.fit(X_tr, y_tr)

gap = np.max(np.abs(reference.coef_ - betas[0.5]))

ConsoleOut [21]
CollapseExpandour intercept, slope:      2.073922,  1.057594
scikit-learn intercept, slope:  2.073922,  1.057594
largest absolute difference: 8.88e-16

The two fits agree to solver precision for this sample. They minimize the same unpenalized objective. When its optimum is non-unique, different solvers can return different optimal coefficient vectors even without numerical error.

With the model in hand, we evaluate it on the held-out 300 observations. We report the pinball loss at each level and, more importantly for this chapter, the empirical coverage of two nominal prediction bands. The test-share-below-fit column is a per-level marginal calibration check; the two nominal bands aggregate that check over $\tau=0.05$ and $0.95$ and over $\tau=0.25$ and $0.75$ respectively.

PythonIn [22]
CollapseExpandgrid_x = np.linspace(0.0, 10.0, 201)
G = add_intercept(grid_x)

pred_te = {tau: X_te @ betas[tau] for tau in taus}
grid_pred = {tau: G @ betas[tau] for tau in taus}

eval_table = pd.DataFrame(
    {
        "tau": taus,
        "test pinball loss": [pinball_loss(y_te, pred_te[t], t) for t in taus],
        "test share below fit": [
            float(np.mean(y_te <= pred_te[t])) for t in taus
        ],
    }
).round(4)
eval_table

Out[23]:

	tau	test pinball loss	test share below fit
0	0.05	0.3182	0.0767
1	0.10	0.5520	0.1200
2	0.25	0.9933	0.2767
3	0.50	1.2191	0.4733
4	0.75	0.9508	0.7567
5	0.90	0.5269	0.9167
6	0.95	0.3153	0.9467

The third column is a per-level marginal calibration check: for accurate quantiles in a continuous model its population target is $\tau$ . It is not a histogram of a fully estimated probability integral transform. Deviations are sampling noise plus any estimation error in $\beta(\tau)$ . Our linear form is correctly specified here, so there is no model-form approximation error; the deviations happen to be small in this sample. Turning those numbers into band coverage:

PythonIn [24]
CollapseExpandcover_90 = float(np.mean((y_te >= pred_te[0.05]) & (y_te <= pred_te[0.95])))
cover_50 = float(np.mean((y_te >= pred_te[0.25]) & (y_te <= pred_te[0.75])))

ols_beta, *_ = np.linalg.lstsq(X_tr, y_tr, rcond=None)
ols_sigma = float(np.sqrt(np.mean((y_tr - X_tr @ ols_beta) ** 2)))
ols_halfwidth = float(norm.ppf(0.95) * ols_sigma)

ConsoleOut [25]
CollapseExpandnominal 90% band  -> empirical coverage 0.870
nominal 50% band  -> empirical coverage 0.480
OLS constant 90% half-width (normal assumption): 5.25

The 90 percent band covers close to 90 percent of held-out responses and the 50 percent band close to 50 percent. This is a good outcome for this sample. The model form matches the data-generating process, and the independent holdout provides an honest check; neither fact guarantees that a fitted band will attain nominal coverage in a finite sample. The last line foreshadows the contrast we will draw visually: the OLS band has one constant half-width, while ours varies with $x$ .

Key ParametersLink Copied

The quantile level $\tau$ is the parameter that defines the model, not a tuning knob in the usual sense. It selects which feature of the conditional distribution we estimate, and different values of $\tau$ produce different coefficient vectors that are not required to be ordered. Choosing $\tau$ is a decision-theoretic act: pick the level that appears in your cost structure or your reporting requirement, such as 0.95 for a service-level target or 0.05 for a lower reference limit.

For fit_intercept=True, QuantileRegressor adds $\alpha\sum_{j=1}^p|\beta_j|$ to the mean pinball loss and leaves its separately fitted intercept unpenalized. If we instead supply a constant column with fit_intercept=False, that column is treated as an ordinary coefficient and is penalized when alpha>0. Our comparison uses alpha=0, so no coefficient is penalized. This is the $\ell_1$ regularization studied in lasso regularization. It shrinks coefficients toward zero and can zero out predictors entirely, which is useful when the number of predictors is large relative to the sample. Because the penalty acts on coefficient magnitudes, features must be on comparable scales for it to be meaningful; standardize first, as described in standardization.

The solver argument selects the underlying optimizer. "highs" is the default; practical solve time depends on the row and feature counts, sparsity, conditioning, and available hardware. The "highs-ds" and "highs-ipm" methods choose the dual simplex and interior-point variants, which can behave differently on degenerate designs. The fit_intercept flag should be set to False whenever we supply our own intercept column, to avoid fitting a redundant second intercept.

Key MethodsLink Copied

fit(X, y) builds and solves the linear program and stores the coefficient vector in coef_, with an intercept in intercept_ when fit_intercept=True. predict(X) evaluates the fitted linear quantile function at new rows; there is no link function and no inverse link, because the fitted conditional quantile is evaluated directly on the response scale used in training.

We compute the pinball loss ourselves rather than relying on a built-in score. The pinball loss is not a proportion of variance explained, and interpreting it requires knowing $\tau$ . Reporting an $R^2$ -like number at an arbitrary quantile invites confusion. The loss is in the units of $y$ , so a value of $0.30$ at $\tau = 0.1$ means the average weighted residual magnitude is $0.30$ response units, not $30$ percent of anything. For inference on the coefficients (standard errors, confidence intervals, joint tests across $\tau$ ), the usual practice is a bootstrap over observations, refitting every $\tau$ within each replicate so that the correlation between levels is preserved.

Comments

Add to the conversation

A question, a correction, or an example from your own work.

No comments yet. Be the first to share your thoughts!

Key ParametersLink Copied

Key MethodsLink Copied

Comments

Add to the conversation

Stay up to date