Heteroscedasticity and Variance Stabilization

The Exam Grading Analogy

Imagine comparing test scores from two different classes. Class A had a straightforward exam: scores range from 85 to 100, and most students cluster around 93. Class B had a much harder exam: scores range from 40 to 90, and most students cluster around 65.

Now suppose a student from each class scores 85. In Class A, that is an unremarkable result — right around the middle. In Class B, it places the student well above the class average. The raw number 85 is hard to interpret on its own; you need to know the exam's difficulty to give it context. To compare students fairly, you adjust for the difficulty: divide each score by its class's standard deviation, or convert to a percentile within the class.

The same issue arises with model errors. Some predictions are inherently harder (more variable) than others, and a residual of a given size means different things depending on the local difficulty. A residual of 10 at an easy-to-predict point suggests the model performed poorly. But a residual of 10 at a hard-to-predict point may be well within the expected range of variation — not a sign of model failure at all. Without knowing the local difficulty, we cannot tell whether a particular residual is large or small in a meaningful sense.

Variance stabilization is essentially "grading on a curve" for each region of feature space. It adjusts each residual by the local difficulty so that all residuals become directly comparable. Once this adjustment is made, a residual of 10 means the same thing everywhere — whether the prediction was easy or hard.

Variance stabilization: divide each residual by the local difficulty to make all residuals comparable.

Why Leverage Creates Unequal Difficulty

Even when the noise in the data is perfectly constant — every observation has the same amount of random error — the model's prediction errors still have different variances at different points. This is the leverage effect from Part 6. Points far from the training data center have high leverage, and the model is more uncertain about its predictions there. Points near the center have low leverage, and the model is more confident.

The practical consequence is that raw residuals naturally fan out as leverage increases: high-leverage points produce larger residuals on average, not because the data is noisier there, but because the model is less certain. If you plotted residuals against leverage, you would see a funnel shape widening to the right. After dividing by the difficulty factor $\sqrt{1+h}$, this funnel collapses into a horizontal band where all residuals have the same spread. The stabilized residuals can then be compared on a level footing.

Raw residuals fan out with leverage (left). After dividing by the difficulty factor, they collapse into a uniform band (right).

Homoscedasticity vs. Heteroscedasticity

In regression, homoscedasticity means the error variance is constant across all observations:

Heteroscedasticity means the error variance depends on the features:

The function $g$ captures the local difficulty: where $g(x) > 1$, errors are larger than average; where $g(x) < 1$, errors are smaller. In most real data, some form of heteroscedasticity is the rule rather than the exception.

Scale Families

A particularly clean form of heteroscedasticity is the scale family. In a scale family, the shape of the error distribution is the same everywhere — only the scale changes:

where $g$ is a known or estimable function, and the $\eta_i$ are identically distributed with mean 0 and variance 1. The key property is that $\varepsilon_i / \sqrt{g(X_i)}$ has the same distribution for all $i$ — not just the same variance, but the same shape and the same quantiles. This is a stronger condition than equal variance alone, and it is what makes scale families particularly well-suited to conformal prediction.

Variance-Stabilizing Transformations

Given a scale family with known $g$, we can stabilize the variance by dividing out the scale factor. If $Z_i$ has variance proportional to $g(X_i)$, define:

The transformed variables $\tilde{Z}_i$ have constant variance. More importantly, if the $Z_i$ form a scale family, the $\tilde{Z}_i$ are identically distributed — exchangeable in the sense conformal prediction requires. To give a concrete example: if the prediction error at a point with leverage $h = 0.3$ has standard deviation $\sigma\sqrt{1.3}$, dividing by $\sqrt{1.3}$ rescales it back to standard deviation $\sigma$. A point with leverage $h = 0.05$ has standard deviation $\sigma\sqrt{1.05} \approx \sigma$, so the division barely changes it. The transformation applies a large correction where the variance is large and a small correction where the variance is already close to the baseline.

VST in Conformal Prediction

Apply this to conformal prediction. Instead of raw absolute residuals $S_i = |Y_i - \hat{f}(X_i)|$, use the stabilized scores:

If the errors form a scale family, these stabilized scores are identically distributed across the calibration set. Computing the conformal quantile on these stabilized scores yields a threshold that reflects the common scale rather than a mixture of different scales. The resulting prediction interval then adapts to local difficulty when we invert the transformation:

The interval is wider where $g(x)$ is large (harder predictions) and narrower where $g(x)$ is small (easier predictions).

The variance stabilization pipeline: raw scores are multiplied by the weight function, a quantile is computed, and the interval width is scaled by 1/w(x).

The Weight Function Perspective

An equivalent way to express variance stabilization is through a weight function $w$:

The stabilized score is $\tilde{S}_i = |Y_i - \hat{f}(X_i)| \cdot w(X_i)$, and the prediction interval width at $x$ is $2\hat{q} / w(x)$. Different weight functions correspond to different assumptions about the variance structure:

Two common weight functions and what they assume about the variance structure.

Comparing the Two Approaches

The General Variance-Stabilizing Principle

The general approach to variance stabilization is straightforward: if a random variable $Z$ has variance that depends on a parameter $\theta$, specifically $\text{Var}(Z) = f(\theta)$, then the transformation

stabilizes the variance to a constant. This is a fundamental tool in statistics dating back to the early 20th century, with several classical instantiations:

• Square root transformation for Poisson data: If $Z \sim \text{Poisson}(\lambda)$, then $\text{Var}(Z) = \lambda$, so $Z / \sqrt{\lambda}$ stabilizes. In practice, $\sqrt{Z}$ is used since $\text{Var}(\sqrt{Z}) \approx 1/4$ for large $\lambda$. For example, if one city reports 400 events and another reports 25, the raw counts suggest the first is 16 times larger — but $\sqrt{400} = 20$ and $\sqrt{25} = 5$ have similar relative uncertainty (standard deviation $\approx 0.5$ for both).
• Arcsine transformation for binomial proportions: If $\hat{p} = Z/n$ with $Z \sim \text{Binomial}(n,p)$, then $\text{Var}(\hat{p}) = p(1-p)/n$. The variance peaks at $p = 0.5$ and shrinks toward zero at $p = 0$ or $p = 1$. The transformation $\arcsin(\sqrt{\hat{p}})$ stabilizes the variance to approximately $1/(4n)$, making proportions near the extremes and near the center equally amenable to standard statistical comparisons.
• Log transformation for proportional variance: If $\text{Var}(Z) \propto [\mathbb{E}(Z)]^2$, then $\log(Z)$ stabilizes the variance (the coefficient of variation is constant). This arises naturally with financial returns, gene expression levels, and other quantities that are multiplicative in nature. Taking logarithms converts multiplicative variation into additive variation with constant spread.

Optimal Weight Under Leverage Heteroscedasticity

In the linear model setting, the prediction error variance is exactly $\sigma^2(1+h(x))$, so the variance function is $g(h) = 1+h$. The optimal weight follows immediately:

This gives stabilized scores $\tilde{S}_i = |Y_i - \hat{f}(X_i)| / \sqrt{1 + h(X_i)}$, which have constant variance $\sigma^2$ regardless of leverage.

When Noise Is Also Heteroscedastic

The analysis requires more care when the noise itself is heteroscedastic. Suppose:

where $g$ depends on leverage. Then the total prediction error variance combines the noise variance and the estimation error variance:

The optimal weight for stabilizing the total prediction error is therefore:

The optimal weight depends on whether the noise itself varies with leverage. Two special cases are shown.

Special Cases

Two special cases illustrate the range of behavior:

Case 1: Homoscedastic noise ($g = 1$). The noise variance is constant at $\sigma^2$. The total prediction variance is $\sigma^2(1 + h)$, and the optimal weight is $(1+h)^{-1/2}$. This is the simplest and most common setting. The weight corrects only for the estimation error term, and it is exact at any sample size because the formula $\sigma^2(1+h)$ is an exact identity, not an approximation.

Case 2: Noise proportional to $(1+h)$ ($g = 1+h$). The noise itself grows with leverage — points that are geometrically unusual also happen to be noisier. The total variance is:

and the optimal weight is $(1+2h)^{-1/2}$. The weight must correct for both the noise heteroscedasticity and the estimation error, leading to a steeper decay with $h$.

When Does Variance Stabilization Help?

Variance stabilization improves conformal intervals when two conditions hold simultaneously:

Condition 1: The variance actually varies. If leverage scores are nearly constant across the feature space (all points are roughly equidistant from the training centroid and the intrinsic noise is constant), there is nothing to stabilize. The diagnostic $\hat{\eta} = \text{std}(h) / \text{mean}(h)$ measures leverage heterogeneity: when $\hat{\eta} > 1$, leverage variation is substantial and leverage weighting is likely to help; when $\hat{\eta} < 0.5$, leverage variation is small and vanilla CP performs comparably.

Condition 2: The variance function is known or well-estimated. Dividing by the wrong scale factor can make things worse — it introduces a systematic distortion rather than removing one. The leverage-based weight $w(h) = (1+h)^{-1/2}$ satisfies this condition by construction because it uses the exact variance formula from the OLS prediction error decomposition, not an estimate. No auxiliary model is needed, and there is no estimation error to worry about. Studentized CP, which estimates $g$ from training residuals, is vulnerable to the sign-flip problem (Part 6): the estimator learns $\sigma^2(1-h)$ instead of $\sigma^2(1+h)$, producing corrections that point in the wrong direction.

Decision diagram: variance stabilization helps when leverage varies substantially and the variance function is correctly specified.