Quantile Regression

Beginner: ordinary least squares always answers the same question — "what is the average y I should expect for this x?" Sometimes that is not the question you actually have. A logistics company asking "how late could this delivery realistically run?" does not want the average delivery time — it wants a number that only 5% of deliveries exceed, so it can set an honest SLA. That is a question about a quantile of the delivery-time distribution — the 95th percentile, specifically — not about its mean. Quantile regression fits a model that directly targets whichever quantile you name, by swapping out squared error for a different loss function built exactly for that purpose.

Intermediate: recall from 2.1.2 (Statistical Decision Theory) that a loss function silently names the statistic your model is trying to learn — squared loss names the conditional mean, absolute loss names the conditional median. The pinball loss (also called the quantile loss), written ρ_τ for a chosen level τ ∈ (0, 1), generalizes that idea one step further: it is asymmetric, charging different penalties for under- and over-prediction depending on τ. Fit a linear model by minimizing the average pinball loss instead of squared error, and the function you get out is a genuine estimate of the conditional τ-quantile of y given x — not a mean shifted up or down by guesswork.

Advanced: the real payoff shows up under heteroscedasticity — when the spread of y | x itself changes with x, as it does for delivery times (short routes are predictable, long routes are erratic) or income (entry-level salaries cluster tightly, senior salaries scatter widely). In that regime, quantiles at different τ are not parallel shifts of the mean line — they fan out, and can even have different slopes. Fitting several quantile regressions at once, say τ = 0.1, 0.5, 0.9, recovers that fanning shape directly from data, which a single mean-plus-constant-margin model cannot do even in principle.

Worked intuition — the pinball loss, one residual at a time: pick τ = 0.9 and imagine two candidate predictions for a delivery that actually took 60 minutes. Predict q = 50 (an under-prediction, residual r = 10): the loss charged is τ·r = 0.9 × 10 = 9. Predict q = 70 instead (an over-prediction, residual r = -10): the loss is (τ-1)·r = -0.1 × -10 = 1. Same-sized miss, nine times the penalty for guessing too low versus too high — exactly the asymmetry a 90th-percentile estimate needs, since underestimating "how late it could get" is the expensive mistake there, and it is what pulls the fitted line up toward the top of the point cloud rather than through its middle.

The from-scratch tabs minimize average pinball loss directly via subgradient descent — the loss is piecewise linear, so it has no single derivative at a zero residual, only a subgradient, which is all that descent needs on a convex objective. The library tab solves the equivalent linear program statsmodels sets up internally and should land on essentially the same three fitted lines, confirming the manual implementation is minimizing the correct objective rather than something that merely looks similar.

#include <cmath>
#include <cstdio>
#include <vector>
#include <random>

struct Point { double x, y; };

double pinballLoss(double residual, double tau) {
    return residual >= 0.0 ? tau * residual : (tau - 1.0) * residual;
}

double avgPinballLoss(const std::vector<Point>& pts, double slope, double intercept, double tau) {
    double total = 0.0;
    for (const auto& p : pts) total += pinballLoss(p.y - (intercept + slope * p.x), tau);
    return total / pts.size();
}

// Subgradient of the average pinball loss w.r.t. (slope, intercept). Piecewise
// linear, so there is no single derivative at a zero residual -- the
// subgradient below picks a consistent branch there, which is all
// subgradient descent requires to converge on a convex objective.
void subgradient(const std::vector<Point>& pts, double slope, double intercept,
                  double tau, double& gSlope, double& gIntercept) {
    gSlope = 0.0;
    gIntercept = 0.0;
    for (const auto& p : pts) {
        double r = p.y - (intercept + slope * p.x);
        double sign = (r > 0.0) ? -tau : (1.0 - tau);
        gIntercept += sign;
        gSlope += sign * p.x;
    }
    gSlope /= pts.size();
    gIntercept /= pts.size();
}

void fitQuantileLine(const std::vector<Point>& pts, double tau, int steps, double lr,
                      double& slope, double& intercept) {
    slope = 0.0;
    intercept = 0.0;
    for (const auto& p : pts) intercept += p.y;
    intercept /= pts.size();
    for (int i = 0; i < steps; ++i) {
        double gSlope, gIntercept;
        subgradient(pts, slope, intercept, tau, gSlope, gIntercept);
        slope -= lr * gSlope;
        intercept -= lr * gIntercept;
    }
}

int main() {
    std::mt19937 rng(7);
    std::uniform_real_distribution<double> unif(0.0, 100.0);
    std::normal_distribution<double> stdNormal(0.0, 1.0);

    const double trueIntercept = 15.0, trueSlope = 0.5;
    std::vector<Point> points;
    for (int i = 0; i < 200; ++i) {
        double x = unif(rng);
        double sigma = 3.0 + 0.25 * x;
        double y = trueIntercept + trueSlope * x + stdNormal(rng) * sigma;
        points.push_back({x, y});
    }

    for (double tau : {0.1, 0.5, 0.9}) {
        double slope, intercept;
        fitQuantileLine(points, tau, 3000, 0.01, slope, intercept);
        double loss = avgPinballLoss(points, slope, intercept, tau);
        std::printf("tau=%.1f  y_hat = %.2f + %.3f*x   avg pinball loss=%.3f\n",
                     tau, intercept, slope, loss);
    }
    return 0;
}