Answer

Mean, median or quantile: which point prediction should a model give?

Inspired by a question on Cross Validated · October 30, 2026

regressiondistributions

The short answer

A fitted probabilistic model gives you a whole distribution for the outcome at each value of the predictors. A point prediction squeezes that distribution into one number, and which number is best depends entirely on how a wrong guess will be penalised. Decide the loss first, then report the summary that minimises it:

• Squared error (MSE, RMSE): predict the mean of the distribution.
• Absolute error (MAE): predict the median.
• Pinball (quantile) loss at level τ, which charges τ for every unit you under-predict and 1 − τ for every unit you over-predict: predict the τ quantile. For example, τ = .9 suits a stock level where running out costs nine times as much as holding one unit too many.
• 0–1 loss on a discrete outcome: predict the mode, the single most probable value.

When the distribution is symmetric, the mean and median coincide and the choice hardly matters. When it is skewed, as incomes, waiting times, costs and counts usually are, they can be far apart, and the wrong choice costs real accuracy.

Fitting the distribution and choosing the prediction are two steps

It is easy to merge these steps, because the most familiar likelihoods have a point-prediction twin. Maximising a normal likelihood for the errors is the same as least squares, which targets the conditional mean. Maximising a Laplace (double exponential) likelihood is the same as least absolute deviations, which targets the conditional median. But the likelihood is really an assumption about the shape of the noise, and its job is to estimate the distribution well.

So what happens if you fit with a Laplace likelihood and are then scored on squared error? Your model estimates the conditional median. If the true noise is symmetric, the median equals the mean, and nothing is wrong: you have only chosen a different estimator of the same line, which can even be more efficient when the noise has heavy tails. If the true noise is skewed, the median and mean differ, and a median-targeted prediction is biased for the quantity squared error rewards. The fix is not necessarily to change the likelihood. If the fitted model describes the full distribution well (scale and shape as well as location), you can compute its mean and report that instead.

A useful way to keep this straight: choose the likelihood to describe the data, and choose the summary to suit the decision. If you only ever need one kind of summary, you can skip the distribution and fit that summary directly with the matching loss (least squares for a mean, quantile regression for a quantile). If you need several, or need intervals too, model the distribution and read each summary off it.

See it in R and Python

Part 1 takes one right-skewed (lognormal) sample of 1,000 values and searches for the single number that minimises each loss. Part 2 fits a lognormal regression, which is an ordinary linear regression of log y on x and is the maximum likelihood fit for that model, to 500 training cases. From that one fitted distribution it makes three predictions: the median exp(μ̂), the mean exp(μ̂ + σ̂²/2) and the 90th percentile exp(μ̂ + 1.2816σ̂). Part 3 scores all three on 100,000 new cases under each loss.

set.seed(677464)

# 1. One skewed sample: which single number minimises each loss?
y <- rlnorm(1000, meanlog = 1, sdlog = 0.75)
sq   <- function(c) mean((y - c)^2)                      # squared error
ab   <- function(c) mean(abs(y - c))                     # absolute error
pin9 <- function(c) mean(pmax(0.9 * (y - c), -0.1 * (y - c)))   # pinball, tau = .9
best <- sapply(list(squared = sq, absolute = ab, pinball_90 = pin9),
               function(f) optimize(f, range(y), tol = 1e-8)$minimum)
round(rbind(minimiser = best,
            summary = c(mean(y), median(y), quantile(y, 0.9))), 3)

# 2. Regression: fit one distribution, then read off different predictions
sim <- function(n) {
  x <- runif(n, 0, 2)
  data.frame(x, y = exp(1 + 0.5 * x + rnorm(n, 0, 0.75)))   # lognormal noise
}
train <- sim(500); test <- sim(1e5)
fit <- lm(log(y) ~ x, data = train)    # MLE of a lognormal regression
s <- sigma(fit)
round(c(coef(fit), sigma = s), 3)
mu <- predict(fit, newdata = test)
pred <- list(median = exp(mu), mean = exp(mu + s^2 / 2),
             q90 = exp(mu + qnorm(0.9) * s))

# 3. Score every prediction under every loss on new data
u <- function(p) test$y - p
scores <- sapply(pred, function(p) c(
  MSE = mean(u(p)^2), MAE = mean(abs(u(p))),
  pinball_90 = mean(pmax(0.9 * u(p), -0.1 * u(p)))))
round(scores, 3)
round(sapply(pred, function(p) mean(test$y <= p)), 3)   # share of y below prediction

import numpy as np
from scipy import optimize, stats

rng = np.random.default_rng(677464)

# 1. One skewed sample: which single number minimises each loss?
y = rng.lognormal(mean=1, sigma=0.75, size=1000)
losses = {"squared": lambda c: np.mean((y - c)**2),
          "absolute": lambda c: np.mean(np.abs(y - c)),
          "pinball_90": lambda c: np.mean(np.maximum(0.9 * (y - c), -0.1 * (y - c)))}
best = {k: optimize.minimize_scalar(f, bounds=(y.min(), y.max()), method="bounded",
                                    options={"xatol": 1e-8}).x for k, f in losses.items()}
print("minimiser", {k: round(v, 3) for k, v in best.items()})
print("summary  ", np.round([y.mean(), np.median(y), np.quantile(y, 0.9)], 3))

# 2. Regression: fit one distribution, then read off different predictions
def sim(n):
    x = rng.uniform(0, 2, n)
    return x, np.exp(1 + 0.5 * x + rng.normal(0, 0.75, n))   # lognormal noise

x_tr, y_tr = sim(500)
x_te, y_te = sim(100_000)
X = np.column_stack([np.ones_like(x_tr), x_tr])
b, *_ = np.linalg.lstsq(X, np.log(y_tr), rcond=None)   # MLE of a lognormal regression
s = np.sqrt(np.sum((np.log(y_tr) - X @ b)**2) / (len(y_tr) - 2))
print("coef", np.round(b, 3), "sigma", round(s, 3))
mu = b[0] + b[1] * x_te
pred = {"median": np.exp(mu), "mean": np.exp(mu + s**2 / 2),
        "q90": np.exp(mu + stats.norm.ppf(0.9) * s)}

# 3. Score every prediction under every loss on new data
for name, p in pred.items():
    u = y_te - p
    print(name, "MSE", round(np.mean(u**2), 3), "MAE", round(np.mean(np.abs(u)), 3),
          "pinball_90", round(np.mean(np.maximum(0.9 * u, -0.1 * u)), 3),
          "share below", round(np.mean(y_te <= p), 3))

          squared absolute pinball_90
minimiser    3.62    2.509      7.308
summary      3.62    2.508      7.307
(Intercept)           x       sigma 
      1.038       0.414       0.768 
           median   mean    q90
MSE        34.536 31.074 63.187
MAE         3.379  3.532  6.911
pinball_90  2.402  1.875  1.222
median   mean    q90 
 0.473  0.631  0.895

All figures come from one seeded run of the R code. Python's random numbers differ from R's, so its figures differ slightly, but they show the same pattern: each loss is again minimised by its own summary, and each prediction wins under its own loss.

What the numbers show

• Each loss picks out its own summary. In the single sample, the number minimising squared error was 3.620, exactly the sample mean. The minimisers of absolute and pinball loss (2.509 and 7.308) matched the median (2.508) and the 90th percentile (7.307) to within 0.001, the small gap arising because those two losses are flat between neighbouring data points.
• On skewed data the summaries are far apart. The mean sat about 44% above the median, and the 90th percentile was nearly three times the median.
• The fit is close to the truth. The regression recovered the true intercept, slope and noise SD reasonably well (1.038, 0.414 and 0.768 against 1, 0.5 and 0.75).
• One fitted model, three different winners. Predicting its mean gave the lowest MSE (31.07, against 34.54 for its median). Predicting its median gave the lowest MAE (3.38, against 3.53 for its mean). Predicting its 90th percentile gave by far the lowest pinball loss (1.22) but its MSE (63.19) was roughly double that of the mean prediction.
• The coverage tells the same story. About 47% of new outcomes fell at or below the median prediction, 63% at or below the mean prediction and 89.5% at or below the 90th-percentile prediction. A mean forecast for a right-skewed outcome is exceeded less than half the time, which is correct for its purpose but surprising if you expected a fifty-fifty guess.

Choosing in practice

Work backwards from the use. If totals matter (a budget, total demand, total cost), the mean is the right target, because means add up and medians do not. If you want a typical case that half of outcomes fall below, use the median. If the costs of over- and under-predicting are unequal, use the quantile that balances them: with under-prediction costing a per unit and over-prediction costing b, the best quantile is a/(a + b).

Then evaluate with the matching loss. Comparing a median forecast by RMSE, or a mean forecast by MAE, rewards the wrong thing and can make a worse model look better. If you care about the whole distribution, score it as a distribution, for example by the log likelihood on held-out data. When the outcome is modelled on the log scale, remember that exp(fitted value) is a median, not a mean; this is a common source of biased totals, and our page on when to log transform data has more on it. For how RMSE relates to the spread of the data, see RMSE vs standard deviation, and for choosing an analysis in the first place try the test chooser or the DASS blog.

How to report MSE and MAE in APA style (7th edition)

In the method section, say which summary the model predicted and why, and name the loss used to evaluate it. In the results, give the error measures on held-out data, in the outcome's units. Using the example above:

"Outcomes were modelled with a lognormal regression (linear regression of log y on x, b = 0.41). Because accuracy was evaluated by mean squared error (MSE), predictions were the conditional means implied by the model. On 100,000 held-out cases, the mean predictions gave MSE = 31.07, compared with 34.54 for the conditional medians, whereas the medians gave the lower mean absolute error (MAE = 3.38 vs. 3.53)."

Spell out each abbreviation at first use, and when several models or predictors are compared, put their error measures side by side in a table with a note naming the evaluation sample.

Related tools and guides

• APA 7 formatter for regression
• RMSE vs standard deviation: what is the difference?
• When should you log transform data?
• How to interpret regression output in R
• Which statistical test should I use?
• Quantile regression (Wikipedia)

More answered questions

• Skewed data with zeros: is log(x + 1) OK, and what else works?
• Why does a t test use the t distribution and ANOVA the F distribution?
• Where should histogram bin edges go when your data are rounded?

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.