Predicting Percentiles

A percentile is a quantile written as a percentage: the 25th percentile is the value below which a quarter of the outcomes fall, the 50th is the median, the 90th is the value below which nine tenths fall. Predicting percentiles means estimating those cut points of the outcome distribution, usually conditional on predictors.

Ordinary regression predicts the mean. Percentiles describe the rest of the distribution, and the reason to want them is that decisions are rarely made at the mean: a retailer sets stock from a high percentile of demand, a risk desk reads value-at-risk off a low percentile of returns, a hospital plans capacity from a high percentile of wait times.

Three routes produce them. Quantile regression models the conditional quantile directly, \(Q_\alpha(y \mid x) = x^{\top}\beta_\alpha\), fitted separately at each quantile level. Distributional forecasting fits a family and reads cut points off the CDF — if \(y \sim N(\mu, \sigma^2)\), the 95th percentile is \(\mu + 1.645\sigma\). Empirical methods — bootstrapping, posterior samples, ensembles — collect predictive samples and take percentiles of them.

A percentile is a promise about how often you are wrong

A model reporting a 10th percentile of 85, a median of 100 and a 90th of 120 is usually read as “about 100, as low as 85, as high as 120.” That reading loses the only thing the numbers actually pin down, which is a frequency.

import numpy as np
from scipy.stats import norm, lognorm

D = lognorm(s=0.30, scale=100.0)
"""Daily demand for one product. Percentiles of this are what a buffer
   is set from."""

print(f"{'percentile':>12}{'stock level':>13}{'expected days short':>22}{'extra vs the median':>21}")
for p in (0.50, 0.80, 0.90, 0.95, 0.99, 0.999):
    print(f"{p:>12.3f}{D.ppf(p):>13.2f}{365 * (1 - p):>22.1f}{D.ppf(p) - D.ppf(0.50):>21.2f}")
print(f"  going from 0.90 to 0.99 costs {D.ppf(0.99) - D.ppf(0.90):.2f} more units,"
      f" {(D.ppf(0.99) - D.ppf(0.90)) / (D.ppf(0.90) - D.ppf(0.50)):.2f} times")

rng = np.random.default_rng(0)
n, stores = 400_000, 10
z = norm.ppf(rng.random((n, stores)))
common = norm.ppf(rng.random((n, 1)))
print(f"\n  and percentiles do not add. {stores} stores, same demand each:")
print(f"{'correlation':>13}{'sum of each 90th':>19}{'90th of the total':>20}{'over-stock':>13}"
      f"{'pooled coverage':>18}")
for r in (0.0, 0.3, 0.7, 1.0):
    mix = np.sqrt(1 - r) * z + np.sqrt(r) * common
    d = lognorm(s=0.30, scale=100.0).ppf(norm.cdf(mix))
    each = np.quantile(d, 0.90, axis=0).sum()
    total = np.quantile(d.sum(axis=1), 0.90)
    print(f"{r:>13.1f}{each:>19.2f}{total:>20.2f}{each - total:>13.2f}"
          f"{np.mean(d.sum(axis=1) <= each):>18.4f}")

print(f"{'sample size':>13}{'sd of 50th':>13}{'sd of 90th':>13}{'sd of 99th':>13}{'99th vs 50th':>15}")
for m in (100, 400, 1600, 6400):
    s = D.ppf(rng.random((4000, m)))
    q = [np.quantile(s, p, axis=1).std(ddof=1) for p in (0.50, 0.90, 0.99)]
    print(f"{m:>13,}{q[0]:>13.4f}{q[1]:>13.4f}{q[2]:>13.4f}{q[2] / q[0]:>15.2f}")
#   percentile  stock level   expected days short  extra vs the median
#        0.500       100.00                 182.5                 0.00
#        0.800       128.72                  73.0                28.72
#        0.900       146.88                  36.5                46.88
#        0.950       163.80                  18.3                63.80
#        0.990       200.95                   3.7               100.95
#        0.999       252.71                   0.4               152.71
#   going from 0.90 to 0.99 costs 54.07 more units, 1.15 times
#
#   and percentiles do not add. 10 stores, same demand each:
#   correlation   sum of each 90th   90th of the total   over-stock   pooled coverage
#           0.0            1468.66             1178.67       289.98            0.9998
#           0.3            1468.00             1299.35       168.65            0.9745
#           0.7            1467.39             1404.49        62.90            0.9274
#           1.0            1468.24             1468.24        -0.00            0.9000
#   sample size   sd of 50th   sd of 90th   sd of 99th   99th vs 50th
#           100       3.7491       7.4115      18.3900           4.91
#           400       1.8649       3.7348      10.2234           5.48
#         1,600       0.9282       1.8703       5.4544           5.88
#         6,400       0.4787       0.9377       2.8280           5.91

Three things the reading misses

With calibrated continuous demand, stock at the 90th percentile has a 10% exceedance probability, or roughly 36.5 expected exceedance days in 365 comparable days. That is not a guarantee for a particular year. Discrete demand can have less than 10% exceedance at that cut point. In the illustrated distribution the 90th-to-99th increment costs more stock than the median-to-90th increment; other distributions need not behave that way.

Quantiles generally do not add. The table compares the sum of local stocks with total demand and obtains 99.98% total-demand coverage for independent stores. That reading requires freely pooled stock or timely transfer: it is not the probability that every store avoids a stockout. With fixed local stocks and independent continuous demands, each store has 90% coverage and all ten simultaneously have 0.9^10 = 34.87%. The nearly 290-unit saving concerns the pooled target. Quantiles add for comonotonic demands; the Gaussian example reaches that case at correlation 1, but comonotonicity is not generally equivalent to Pearson correlation 1.

Third, a cut point far out in the tail is estimated from few observations. At every sample size above, the 99th percentile is estimated about five times less precisely than the median, and at \(n = 100\) it sits between the two largest values that happened to be drawn. Quantile regression, a fitted distribution and a bootstrap all face this — fitting a family borrows strength from the middle to steady the tail, at the price of assuming the family is right there.

The median is a probability cut point, not generally the most likely value. The illustrated lognormal has mode below median below mean. Skewness alone does not guarantee that ordering for every distribution, nor does symmetry alone guarantee a unique mode at the mean.

Where percentile forecasts are the natural output

A prediction interval is two percentiles: a 90% interval is the 5th and 95th. Value-at-risk is a lower percentile of a return distribution. Supply chain, electricity load and weather extremes are all threshold decisions, and survival analysis reports percentile life expectancy. In each case the quantity asked for is a cut point rather than an average, and estimating it directly avoids assuming a shape that the middle of the data cannot check.

Two habits go with that. Score quantile forecasts with pinball loss at the levels you actually use, and check that the realized frequencies match the levels claimed — a 90th percentile that is exceeded 20% of the time is a broken promise, not a conservative one. And when separately fitted quantiles cross, which they can, the model has reported an impossible distribution and the crossing has to be resolved rather than ignored.

References

• SciPy documentation, scipy.stats.lognorm and scipy.stats.norm. docs.scipy.org
• NumPy documentation, numpy.quantile. numpy.org