Calibrated Multivariate Regression with Localized PIT Mappings

In the univariate case, where $\cal{Y}\subseteq\mathbb{R}$, several notions of calibration exist within the literature (e.g., Gneiting and Resin; 2023). Here, we focus on marginal calibration and probability calibration, which are two choices commonly considered in practice.

$\widehat{F}$ is said to be marginally calibrated (Gneiting et al.; 2007) if

That is, the average predictive CDF $\frac{1}{n_{\text{val}}}\sum_{i=1}^{n_{\text{val}}}\widehat{F}(y\mid\bm{x}_{\text{val}}^{(i)})$ matches with the empirical CDF of the observations $\frac{1}{n_{\text{val}}}\sum_{i=1}^{n_{\text{val}}}\mathds{1}_{\left\{y_{\text{val}}^{(i)}\leq y\right\}}$ asymptotically for all $y\in\cal Y$. Hence, Gneiting et al. (2007) suggest plotting the average predictive CDF versus the empirical CDF to graphically assess marginal calibration.

The random variable

where $\operatorname{\mathcal{V}}\sim\operatorname{U}(0,1)$ and $\widehat{F}(Y^{-}\mid\bm{X})$ is the left-handed limit of $\widehat{F}(y\mid\bm{X})$ as $y$ approaches $Y$ from below, is the randomized PIT value (Czado et al.; 2009). Note that $P$ depends both on $F(Y,\bm{X})$ and the model $\widehat{F}(Y\mid\bm{X})$. If $\widehat{F}$ is continuous, $P$ is not randomized

$\widehat{F}$ is said to be probability calibrated if $P\sim\operatorname{U}(0,1)$. For $i=1,\dots,n_{\text{val}}$ let $\nu^{(i)}$ be independent uniform random variables on $[0,1]$ and write

$p^{(i)}$ is an empirical evaluation of (1) across the validation set. Thus, probability calibration can be graphically checked by plotting a histogram of $\{p^{(i)},i=1,\dots,n_{\text{val}}\}$. Formal tests for uniformity of the PIT values based on the Wasserstein distance (Zhou et al.; 2021; Zhao et al.; 2020) and the Cramér-von Mises distance (Kuleshov et al.; 2018) are popular alternatives to graphical checks. Gneiting et al. (2007) show that probability calibration is under mild conditions equivalent to quantile calibration (Kuleshov et al.; 2018), which requires

where $\widehat{F}^{-1}(\cdot\mid\bm{x}_{\text{val}}^{(i)})$ denotes the generalized inverse of $\widehat{F}(\cdot\mid\bm{x}_{\text{val}}^{(i)})$. This perspective has the nice interpretation that prediction intervals derived from $\widehat{F}$ have the correct coverage.

Histograms of the PIT values can also be used for model criticism as they indicate the type of miscalibration at hand. For example, U-shaped histograms indicate overconfidence, while triangular shapes indicate a biased model (Gneiting et al.; 2007).

Several techniques to recalibrate a potentially miscalibrated model $\widehat{F}$ in a post-hoc step exist in the literature. Here, we describe a simple method for doing this due to Kuleshov et al. (2018). Write $G(p)$ for the distribution function of the PIT values (1), where dependence on $F$ and $\widehat{F}$ is left implicit in the notation. It is easy to check that $G(\widehat{F}(y\mid\bm{x}))$ is a distribution function for every $\bm{x}\in\cal X$ and probability calibrated with respect to $F(Y,\bm{X})$. In practice, the distribution function $G(p)$ is not known, and it must be estimated from $\mathcal{D}_{\text{val}}$. Kuleshov et al. (2018) suggest using a method based on isotonic regression. Dheur and Ben Taieb (2023) extend this idea and consider KDEs. Among other choices, they propose to use

which is the empirical CDF from the PIT values over the validation set $\cal{D}_{\text{val}}$. An alternative approach to recalibration for regression models is given by Song et al. (2019). Recently, Torres et al. (2024) proposed nonparametric local recalibration for neural networks. Their approach uses a fast KNN algorithm to localize the recalibration and can be used in any layer of the neural network scaling to potentially high-dimensional feature spaces $\cal X$.

Extending the different notions of calibration from the univariate to the multivariate case, $\mathcal{Y}\subseteq\mathbb{R}^{d}$, is not straightforward. One reason for this is that the multivariate integral transformation

is, in contrast to the univariate case, generally not uniformly distributed (e.g., Genest and Rivest; 2001), but follows the so-called Kendall distribution of $F$. The Kendall distribution depends only on the copula of the multivariate probability measure, and thus summarizes the dependence structure of $F$. Based on this observation, Ziegel and Gneiting (2014) introduce copula probability integral transform (CopPIT) values as analogous to the univariate PIT values described in (1). The CopPIT values are given as

where $\Upsilon\sim\operatorname{U}(0,1)$, $(\bm{Y},\bm{X})\sim F(\bm{Y},\bm{X})$, and $\mathcal{K}_{\bm{X}}$ denotes the Kendall distribution of $\widehat{F}(\bm{Y}\mid\bm{X})$. $\widehat{F}$ is said to be copula calibrated if the CopPIT values are uniformly distributed on the unit interval (Ziegel and Gneiting; 2014). In this way, copula calibration can be seen as a multivariate extension to probability calibration. In particular for $d=1$, $\mathcal{K}_{\bm{X}}$ is the uniform distribution on $[0,1]$, so that (3) is equal to $\eqref{eq:pit}$. As in the univariate case, let

denote the empirical CopPIT values from the validation set, $i=1,\dots,n_{\text{val}}$, where $\upsilon^{(i)}$ are independent uniform variates on $[0,1]$. Again, copula calibration can be assessed by checking uniformity of $\{u^{(i)},i=1,\dots,n_{\text{val}}\}$. However, interpretation of the CopPIT histograms is more challenging than in the univariate case, as they not only summarize potential miscalibration of the dependence structure, but also of the marginal distributions. However, in the special case that all margins of $\widehat{F}$ are uniformly probability calibrated, the CopPIT values summarize miscalibration of the copula of $\widehat{F}$ only (Ziegel and Gneiting; 2014). Thus, in practice it is sensible to assess multivariate calibration by checking for univariate calibration of each marginal in terms of the marginal PIT values (1) and copula calibration in terms of the CopPIT values (3).

Ziegel and Gneiting (2014) also introduce Kendall calibration

as the multivariate analogue to marginal calibration. Kendall calibration can be assessed by a so called Kendall diagram, which is a scatter plot of the empirical left hand side versus the empirical right hand side of (4) for different values of $\omega$.

Both copula calibration and Kendall calibration necessitate the derivation of the Kendall distributions $\mathcal{K}_{\bm{x}}$. The Kendall distribution can be calculated in closed form only for a few special cases (e.g., Genest and Rivest; 2001). So, in practice, $\mathcal{K}_{\bm{x}}$ in (3) and (4) is replaced by an approximation given as the empirical CDF of the pseudo observations (Barbe et al.; 1996)

where $\bm{y}_{j}=(y_{j1},\dots,y_{jd})\preceq\bm{y}_{k}=(y_{k1},\dots,y_{kd})$ if $y_{jl}\leq y_{kl}$ for all $l=1,\dots,d$ and $\bm{y}_{1},\dots,\bm{y}_{m}$ is a large sample from $\widehat{F}(\bm{y}\mid\bm{x})$.

Suppose that we have a mapping on the feature space, $h:\mathcal{X}\to\mathbb{R}^{d}$. The purpose of the function $h$ is to reduce the dimension of $\bm{x}$ and we use $\lVert h(\bm{x})-h(\bm{x}^{\prime})\rVert$ to measure the similarity of the feature vectors $\bm{x}$ and $\bm{x}^{\prime}$. Let

If $N_{k}(\bm{x})$ is a sufficiently small neighbourhood around $\bm{x}$, $\{\bm{p}^{(i)},i\in N_{k}(\bm{x})\}$ approximates a sample from $\bm{P}=(P_{1},\dots,P_{d})$, where $P_{l}$ is the PIT value for the $l$-th response of $\widehat{F}(\bm{Y}\mid\bm{x})$ as given in (1). Let $G_{\bm{x}}$ denote the joint distribution of $\bm{P}$ given $\bm{x}$ with marginal distributions $G_{\bm{x},l}$, $l=1,\dots,d$. Theoretical properties of $G_{\bm{x}}$ were studied in Rodrigues et al. (2018). In particular, $G_{\bm{x}}(\widehat{F}(\bm{y}\mid\bm{x}))=\left(G_{\bm{x},1}(\widehat{F}_{1}(y_{1}\mid\bm{x})),\dots,G_{\bm{x},d}(\widehat{F}_{d}(y_{d}\mid\bm{x}))\right)$ has probability calibrated marginals following the same arguments as for the univariate recalibration techniques described in Section 2.1. Note that the use of $G_{\bm{x},l}$ instead of the global unconditional distribution $G_{l}$ as considered in Kuleshov et al. (2018) and Dheur and Ben Taieb (2023) gives a stronger form of calibration, as the resulting model is locally, that is conditional on $\bm{x}$, probability calibrated. In addition to the marginal information, $G_{\bm{x}}$ also matches the dependence structure of $\widehat{F}(\bm{Y}\mid\bm{x})$ under $F(\bm{Y}\mid\bm{x})$. For a given $\bm{x}\in\cal X$ and continuous marginals $\widehat{F}_{l}(y\mid\bm{x})$, $p_{l}$ is a non-random transformation of $y_{l}$ and, in particular, $\bm{p}$ is an invertible transformation of $\bm{y}$. The Kendall distribution is invariant under such transformations and thus the CopPIT value $u\mid\bm{x}$ could be calculated purely on $\bm{p}\mid\bm{x}$, without access to $\bm{y}\mid\bm{x}$. Also, $\bm{P}\mid\bm{x}$ has copula $C_{\bm{x}}$, which is the copula of $F(\bm{Y}\mid\bm{x})$ and, from the arguments above, $G_{\bm{x}}(\widehat{F}(\bm{Y}\mid\bm{x}))$ has probability calibrated marginals. Following Ziegel and Gneiting (2014) this implies copula calibration.

Thus, for a given $\bm{x}\in\cal X$, a sample of size $k$ from an approximately calibrated predictive distribution $\widetilde{F}(\bm{Y}\mid\bm{x})$ can be generated as

Here, $\bm{p}^{(i)}=(p_{1}^{(i)},\dots,p_{d}^{(i)})$ with $p^{(i)}_{l}$ the empirical PIT value for the $l$-th marginal distribution $\widehat{F}_{l}(Y_{l}\mid\bm{X})$ evaluated on the $i$-th entry of the validation set. If $\widehat{F}^{-1}_{l}(\cdot\mid\bm{x})$ is not available in closed form it can be easily approximated using a sample from $\widehat{F}_{l}(\bm{Y}\mid\bm{x})$ making our approach model free.

However, (5) is only an approximation to $G_{\bm{x}}(\widehat{F}(\bm{y}\mid\bm{x}))$ and we will illustrate how well this works in practice on a number of simulated and real data examples in Sections 4 and 5.

We can think of the nearest neighbour approach introduced in Section 3.1 as obtaining an approximate sample from $\bm{P}\mid\bm{X}=\bm{x}$ for a target feature vector $\bm{x}$, and then transforming back to the original space of the responses to obtain approximate samples of $\bm{Y}\mid\bm{X}=\bm{x}$. In the nearest neighbour approach, no explicit expression for $\widetilde{F}(\bm{Y}\mid\bm{x})$ is constructed, and the number of potential draws from $\widetilde{F}(\bm{Y}\mid\bm{x})$ is restricted by $k$ heavily depending on $n_{\text{val}}$. However, some applications require calculating complex summary statistics from the potentially intricate distribution $\widetilde{F}$, which necessitates the ability to draw arbitrary large samples from the recalibrated model. Thus, we propose a similar method to the KNN approach using normalizing flows to draw approximate samples from $\bm{P}\mid\bm{x}$.

The basic idea is as follows. Let $\rho(\bm{z})$ be a reference density with respect to the Lebesgue measure on $\mathbb{R}^{d}$, which we take to be the standard normal density. We consider a bijective transformation $T_{\bm{\zeta}}(\bm{z}\mid\bm{x})$, and transform $\bm{Z}\sim\rho(\bm{z})$ to a random vector $\bm{P}\mid\bm{x}$, where $\bm{\zeta}$ is a set of learnable parameters. The density of $\bm{P}\mid\bm{x}$ is thus approximated as

where $T^{-1}_{\bm{\zeta}}(\bm{p}\mid\bm{x})$ is the inverse of $T_{\bm{\zeta}}(\bm{z}\mid\bm{x})$, and $J_{T^{-1}_{\bm{\zeta}}}(\bm{p}\mid\bm{x})$ is its Jacobian matrix. Based on this, the parameter ${\bm{\zeta}}$ is learned using observations $(\bm{p}^{(i)},\bm{x}^{(i)}_{\text{val}})$, where $\bm{p}^{(i)}$ denotes the vector of PIT values for the marginal distributions evaluated on the validation set. To avoid boundary effects, we consider normalized PIT values $\bm{p}_{N}=\Phi^{-1}(\bm{p})=\left(\Phi^{-1}(p_{1}),\dots,\Phi^{-1}(p_{d})\right)\in\mathbb{R}^{d}$, where $\Phi(\cdot)$ is the CDF of the standard Gaussian distribution, instead of the usual PIT values on $[0,1]$. As for the KNN approach, $\bm{x}$ can be replaced with a lower dimensional representation $h(\bm{x})$ in the construction of $T_{\bm{\zeta}}(\bm{z}\mid\bm{x})$. Having learned ${\bm{\zeta}}$ as $\hat{{\bm{\zeta}}}$, samples from $\bm{P}\mid\bm{x}$ can be generated by sampling $\bm{z}\sim\rho(\bm{z})$ and setting $\bm{p}=T_{\hat{{\bm{\zeta}}}}(\bm{z}\mid\bm{x})$. These samples can be then used similar to (5) to generate samples from the approximately recalibrated predictive distribution $\widetilde{F}(\bm{Y}\mid\bm{x})$. Pseudo code for the full approach is given in Appendix A.

There are many ways to construct suitable transformations $T_{{\bm{\zeta}}}(\bm{z}\mid\bm{x})$ in the literature on transport maps (e.g., Marzouk et al.; 2016) and normalizing flows (e.g., Rippel and Adams; 2013; Yao et al.; 2023) as well as standard software for conditional density estimation. Here, we consider the real-valued non-volume preserving (real-NVP) approach by Dinh et al. (2017), where $T_{{\bm{\zeta}}}(\bm{z}\mid\bm{x})$ is given as a stack of affine coupling layers. We found this to be a satisfactory choice in all examples considered, but alternative approaches might be better suited depending on the structure of the recalibration task at hand.

The normalizing flow can be interpreted as a conditional density estimate for $\bm{P}\mid\bm{x}$. If $d=1$, our approach can therefore be considered a localized version of the KDE based method by Dheur and Ben Taieb (2023).

Foreign currencies constitute a popular class of assets among investors. In so-called Forex trading, traders exchange currencies with the goal of making a profit. The ability to make reliable predictions for currency exchange rates is crucial for a successful trading strategy. To this end, we analyze five time series of daily exchange rates for five currencies relative to the US dollar: the Australian Dollar (AUD), the Chinese Yuan (CNY), the Euro (EUR), the Pound Sterling (GBP), and the Singapore Dollar (SGD). These data span five years from August 01, 2019, to August 01, 2024, and were sourced from Yahoo Finance. Given the high correlation among currency exchange rates, multivariate calibration is a desirable feature for any currency exchange forecasting model.

Ending all forms of malnutrition is one of the sustainable development goals of the United Nations (United Nations; 2015). Here, we consider a sample from the Demographic and Health surveys (www.measuredhs.com) containing $n=24,286$ observations on several factors of undernutrition in India. Following previous analyses of the data (Klein et al.; 2020; Briseño Sanchez et al.; 2024) we consider two responses. The continuous indicator wasting reports weight for height as a z-score and the binary response fever indicates fever within the two weeks prior to the interview. Following Klein et al. (2020) we consider the covariables csex indicating the sex of the child, cage the age of the child in months, breastfeeding the duration of breastfeeding in months, mbmi the body mass index of the mother, and dist the district in India the child lives in.