Conformal intervals under a log transform: the blow-up isn't a back-transform bug

Someone log-transforms their target, fits a model, wraps it in conformal prediction, inverse-transforms the intervals back with expm1, and the upper bound comes out at 350x the point forecast. The natural reaction is "the back-transform is broken." I ran into exactly this framing on a real bug report recently, and the interesting part is that the back-transform is correct — the interval is doing precisely what it should. What looks like a bug is two honest effects stacking. Here's the mental model, with a runnable check.

The setup

You have skewed, positive data (demand, counts, prices). The standard move is to fit on log1p(y) and report on the original scale with expm1. You also want distribution-free intervals, so you add split conformal prediction on top. The conformal machinery runs in the model's scale — here, log space — and hands you interval endpoints lo and hi in log space. You expm1 those two columns and get your original-scale interval.

Then hi explodes. point ≈ 400, lo ≈ 200, hi ≈ 145,000. Surely you transformed something wrong.

The one fact that resolves it

Coverage is invariant under a monotone transform.

If g is strictly increasing (and expm1 is), then for any interval [lo, hi] in log space:

y ∈ [lo, hi]   ⟺   g(y) ∈ [g(lo), g(hi)]

The two events are the same event. So if your log-space interval covers the truth 90% of the time, the expm1-transformed interval [expm1(lo), expm1(hi)] covers the (back-transformed) truth 90% of the time too — exactly, not approximately. The back-transform adds zero error.

The crucial detail: you transform the endpoints, not the scores. expm1(lo) and expm1(hi), never expm1(score) added to a point forecast. Transforming an additive score through a nonlinear map is what actually breaks — but returning endpoint columns and mapping those is the correct operation, and it's usually what libraries already hand you.

So why is hi enormous?

Two things, both honest:

1. Lognormal skew under a convex map. expm1 is convex, so a roughly symmetric interval in log space is supposed to become very asymmetric and wide on the top in the original scale. A right-skewed quantity has a genuinely long upper tail. Some of that 350x is simply correct.

2. An inflated log-space score, then exponentiated. Conformal width is set by an upper quantile of the calibration residuals. On a short, badly-specified, or thin calibration set, that quantile is large and unstable — and exp of a large number is a very large number. This is small-sample variance amplified by the transform, not miscoverage.

Show me

Self-contained, no forecasting library needed — just the conformal mechanism on lognormal data:

import numpy as np
rng = np.random.default_rng(0)

def split_conformal_logspace(cal_true, cal_pred, mean_log, alpha):
    scores = np.abs(cal_true - cal_pred)            # |residual| in log space
    n = scores.size
    k = int(np.ceil((n + 1) * (1 - alpha)))         # finite-sample rank
    qhat = np.inf if k > n else np.sort(scores)[k - 1]
    return mean_log - qhat, mean_log + qhat, qhat

def coverage(y, lo, hi):
    return float(np.mean((y >= lo) & (y <= hi)))

alpha, mu, sigma = 0.10, 6.0, 0.5
cal   = rng.normal(mu, sigma, 2000)
test  = rng.normal(mu, sigma, 20000)

lo, hi, qhat = split_conformal_logspace(cal, np.full_like(cal, mu), mu, alpha)

cov_log  = coverage(test, lo, hi)                          # log space
cov_orig = coverage(np.expm1(test), np.expm1(lo), np.expm1(hi))  # after expm1

print(f"log-space coverage : {cov_log:.4f}")
print(f"orig-space coverage: {cov_orig:.4f}")
print(f"identical          : {cov_log == cov_orig}")
print(f"hi/point ratio     : {np.expm1(hi)/np.expm1(mu):.1f}x")

Output:

log-space coverage : 0.8999
orig-space coverage: 0.8999
identical          : True
hi/point ratio     : 2.3x

Coverage is bit-for-bit identical across scales. Now shrink and corrupt the calibration set the way a misspecified model on a short series would — few effective residuals, inflated variance:

for n_cal, extra in [(2000, 0.0), (30, 0.0), (30, 1.5)]:
    c   = rng.normal(mu, sigma, n_cal)
    pred = np.full(n_cal, mu) + rng.normal(0, extra, n_cal)   # misspecification
    lo, hi, q = split_conformal_logspace(c, pred, mu, alpha)
    print(f"n={n_cal:>4} extra={extra}  qhat={q:5.2f}  hi/point={np.expm1(hi)/np.expm1(mu):7.1f}x")

n=2000 extra=0.0  qhat= 0.82  hi/point=    2.3x
n=  30 extra=0.0  qhat= 0.92  hi/point=    2.5x
n=  30 extra=1.5  qhat= 2.99  hi/point=   19.9x

The upper endpoint balloons with qhat — driven by the calibration set, not the back-transform — and coverage stays valid the whole way.

What to actually do

• Keep transforming endpoints, not scores. The expm1(lo), expm1(hi) step is correct. Don't "fix" it.
• Read a giant upper bound as a diagnostic. It usually means your log-space score quantile is unstable: more calibration history, more windows, or a better-specified model. Inspecting intervals in the modeling scale makes this obvious.
• Prefer an internal transform when the library supports it. If the transform lives inside the model (e.g. a working Box-Cox path), conformity scores are computed and back-transformed consistently and you never hand-roll log1p/expm1 — which sidesteps the whole class of confusion.
• Expect asymmetry, and treat it as honest. For right-skewed targets a wide, one-sided upper interval is the correct answer, not a rendering glitch.

The through-line: a monotone transform is one of the few things conformal prediction handles for free. The coverage guarantee rides along untouched. When an interval looks absurd after expm1, the transform is rarely the culprit — the calibration is where to look.

Disclaimer: This article was drafted with AI assistance and reviewed and edited by the author. The technical design and opinions are my own.