JAPAN: Joint Adaptive Prediction Areas with Normalising Flow
Authors: Eshant English
Christoph Lippert
Venue: ICLR 2026 Poster (2026)
Keywords: Uncertainty Quantification
Normalising Flows
Joint Prediction Areas
Primary area: generative models
TL;DR: An application of Normalising Flows to obtain compact disjoint prediction regions without any geometrical assumption

Abstract: Conformal prediction provides a model-agnostic framework for uncertainty quantification with finite-sample validity guarantees, making it an attractive tool for constructing reliable prediction sets. However, existing approaches commonly rely on residual-based conformity scores, which impose geometric constraints and struggle when the underlying distribution is multimodal. In particular, they tend to produce overly conservative prediction areas centred around the mean, often failing to capture the true shape of complex predictive distributions. In this work, we introduce JAPAN (Joint Adaptive Prediction Areas with Normalising-Flows), a flow-based framework that uses density estimates for several conformal scores. By leveraging flow-based models, JAPAN estimates the (predictive) density and constructs prediction areas by thresholding on the estimated density scores, enabling compact, potentially disjoint, and context-adaptive regions that retain finite-sample coverage guarantees. We theoretically motivate the efficiency of JAPAN and empirically validate it across multivariate regression and forecasting tasks, demonstrating good calibration and tighter prediction areas compared to existing baselines. Furthermore, several density-based conformity scores showcase the flexibility of our proposed framework.

Published as a conference paper at ICLR 2026
JAPAN: JOINT ADAPTIVE PREDICTION AREAS WITH FORMALISING-FLOWS
Eshant English 1.23; Christoph Lippert! 1.4
Hasso Plattner Institute for Digital Engineering, Germany University of Tokyo, Tokyo, Japan
3 Centre for Advanced Intelligence Project, RIKEN, Tokyo, Japan 'Hasso Plattner Institute for Digital Health at Mount Sinai, USA
ABSTRACT
Conformal prediction provides a model-agnostic framework for uncertainty quan- tification with finite-sample validity guarantees, making it attractive tool for constructing reliable prediction sets. However, existing approaches commonly rely residual-based conformity scores, which impose geometric constraints and struggle when the underlying distribution is multimodal. In particular, they tend produce overly conservative prediction areas centred around the mean, often failing to capture the true shape of complex predictive distributions. In this work, we introduce JAPAN (Joint Adaptive Prediction Areas with Kormalising-Floww flow-based framework that uses density estimates for several conformal scores. By leveraging flow-based models, JAPAN estimates the (predictive) density and con- structs prediction areas by thresholding on the estimated density scores, enabling compact, potentially disjoint, and context-adaptive regions that retain finite-samin sample coverage guarantees. We theoretically motivate the efficiency of JAPAN and empir- ically validate it across multivariate regression and forecasting tasks, demonstrating good calibration and tighter prediction areas compared to existing baselines. Fur- thermore, the several density-based conformity scores showcase the ficxibility our proposed framework.
INTRODUCTION
Conformal Prediction (CP) (Gammerman et al. 1998 Shafer & Vovk 2008) offers a model-agnostic
framework for uncertainty quantification with finite-smple coverage guarantees. For user-specified significance level e E (0, 1), CP constructs a prediction region I'. C y (the label space for prediction y) satisfying (2 E (x)) 2 1- i.e. the true output lies in the prediction set with high probability given by the significance level c. This makes CP an appealing choice for enhancing the reliability of predictive models, particularly in safety-critical or high-stakes applications.
Traditionally, CP has been implemented using residual based conformity scores, such as the absolute error ly - 5) While intuitive and effective for univariate response y, it does not remain so in higher dimensions as the residuals become multivariate. The residuals can be mapped to univariate conformity score, but this requires a choice of geometry-e.e -nomm construct spherical regions, and -norms give rise to rectangular regions. However, imposing such geometric constraints on the prediction area may not reflect the true nature of the underlying uncertainty. Additionally, residual- based conformity scores tend to over-concentrate around single mode, resulting in overly broad sets when the predictive distribution is multimodal.
"contact email: sshantenglish@@ppidd. wtw.::/wwhheeleeeeeeeeeeeeiiiiaa cshanteglisttggga..........

---

Published as a conference paper at ICLR 2026
True Density
JAPAN region
CONTRA
PCP Elliptical region
Rectangular
Figure I: For the spiral density, only JAPAN produces compact regions that align with high-density areas, as the prediction area closely follows the spiral, whereas other baselines fail to capture this.
Although CP guarantees validity, practical usefulness also requires efficiency: the size of the prediction region. Ideally, we seek the smallest possible region that maintains coverage, tightly capturing likely outcomes while excluding unnecessary low-density areas. It is already known (Lei et al. 2013) that conformal prediction arcas obtained by thresholding the conditional density ply x) are compact, and can even be the smallest possible when the true threshold (T₂) is not function of x. While the conditional density p(y z) can yield efficient conformal areas in theory, the density is typically unknown in practice. Fortunately, it be approximated using estimators such as flow-based models Paaaaaakarios et al.] 2021; Dinh et al.] 2017), which are flexible generative models that learn a mapping either explicitly or implicitly from base distribution tractable likelihood and can be used to approximate p(y x). Additionally, as we will see, their latent space sampling provides an efficient mechanism for computing the regions, making them attractive in comparison to other density estimation methods such as Kernel Density Estimation (KDE). When the estimated density p(y 2²) is sufficiently accurate, the resulting conformal prediction sets approach the size under the true density threshold while preserving validity. This insight motivates a shift away from distance-based conformity scores and towards density-based ones.
target distribution. They have
In this work, we propose JAPAN (Joint Adaptive Prediction Areas for Normalising- flowss method that uses the (conditional) density estimates as the conformity score computed by flow-based models. Largely, motivate IAPAN-consttuuttd prediction areas by thresholding
densities from calibration set, such as: r (x) == {y E y log p(y .) 2 Te) where т, is the (1 e)-quantile of calibration scores. These areas are geometry-free, avoiding assumptions like ellipsoids or rectangles; potentially disjoint, adapting to the shape of the predictive distribution by only aligning with areas of high density, especially when it is multimodal; and context-adaptive, adjusting with the input æ, particularly in structured settings like time series.
We summarise our main contributions follows:
We propose JAPAN, a flow-based method that uses log-densities as conformity scores, enabling efficient, adaptive, and potentially multimodal prediction areas.
We theoretically show that the regions formed are close to the ones formed by the true unconditional threshold under mild ranking errors.
We empirically evaluate JAPAN across diverse tasks multivariate regression and time series forecasting (while also proposing a new deriresturredddmmmmsrrrtttnn good coverage and tighter areas compared to other baselines.
We highlight the versatility of JAPAN by several density-based conformity scores, including unconditional, conditional, boster, latent-based, and T, adaptive modelling setups.
Throughout the paper, we use the term area as measure of prediction sets, though in higher dimensions this naturally generalises volume. Similarly, use Normalising Flows (NFs) as an umbrella term for the broader class of flow-based generative models.
PRELIMINARIES
2.1 NORMALISING FLOWS
Normalising Flows (NFs) are class of generative models that transform simple base distribution into complex target distribution either by sequence of bijective mappings or by integrating continuous dynamics (e.g.. via learned velocity fields score functions). This includes discrete flows as well as continuous- time flows such as score-based diffusion models (bijective mapping exists through the probability flow ODE) (Song et al. 2021 2020), continuous normalising flows (CNFs)
2

---

where ((,,2) 2) encodes the change in log-volume change due to the underlying flow mapping, e-g., log-detenn deg-deeerminatt Jacobian for discrete flows (Dinh et al., 2017). or a time-integrated divergence term for continuous- time models (Chen et al.] 2018). For the unconditional case, simply remove the Published as a conference paper at ICLR 2026
(Grathwohl et al. 2018) trained with without Flow-Matching interpolants objectives (Albergo & Vanden- Eijjnenn 2022).
Lipman et al. 2022 stochastic
Formally, a NF realise bijective mapping h P2 + RP,, such that h(y) pz (2), where Pz is a simple base distribution, typically standard Gaussian. The mapping can be explicitly learned by sequence of invertible transformations Dinh et al. (017)) or obtained by integrating ODE
defined by a learned velocity field (Lipman et al. 2022) or the score function (Song et al. 2020). For NF the density of y can be computed by the change-of- varriabllss formula:
log p(u x) = log (v, x)) + (1)
context & and have h, only depend on y, and thus yield p(y) with a similar log-likelihood equation.
2.2 INDUCTIVE CONFORMAL PREDICTION
Inductive Conformal Prediction (ICP) PPapadopoulos al. 2002 Vovk et al. 2005 is scalable variant of conformal prediction that enables valid prediction sets by splitting the data into training and calibration sets. Let D ( yt ]== be dataset of i.i.d. samples from unknown distribution PP. This dataset is partitioned into a training set D { 1:)) and a calibration ser ceal { y4 y=))==++++ where n m +
Typically, a regression model h X' - y is fitted on Duain, and a conformity score ai is computed on each calibration point (zii, y:). For test input x" and candidate output y. the test conformity score a" is compared against the calibration scores via p-value:
T()== - 1+1 1
1
1
saa)
:=1
which quantifies how plausible y is under the model and calibration set. The prediction region is then defined as:
Te(x") = fyr E y m(y)) c)-
In the classical univariate regression setting, it is natural to work with nonconformity score, i.e. testing the implausibility of a label instead, and is defined as the absolute residual: as = ly p(xx)) with μ the predictor. In this case, the prediction set (x()) can be constructed directly symmetric interval centred at p(="). ). This interval can be obtained via a simple inversion of the above p-value condition for the residual score such that (x*) := {u(x*) q1-.) where ęı- is the
(1 - e)-quantile of the nonconformity scores {a.)/= based on the absolute residual error.
3 JAPAN: JOINT ADAPTIVE PREDICTION AREAS VIA NORMALISING-FLOWW
We apply JAPAN to two settings: multivariate regression and exchangeable trajectories of time series. Despite differing in structure. from a conformal inference perspective, both share the same core challenge: designing conformity scores that yield valid and efficient prediction areas. A particular difficulty in multivariate outputs, especially for time series, is the presence of inter-dimensional (temporal) dependencies. These make simple geometric constructions like L2 balls or Bonferroni- corrected rectangular areas inefficient, as they do not incorporate inter-dimensional dependencies. Since JAPAN works directly with the distribution of P(Y by leveraging normalising flows, we work with both settings with minimal changes. Note that JAPAN does not require a Normalising Flow to be trained conditionally on X, as we will see later in Section 5 or Section A Normalising Flows can be trained unconditionally, or conditioned on Y, or like a generative model (-IY), however, we use the notion P(Y 1)) first for simplicity. Also, we have slight abuse notation where ær, y, z are variables and the subscript versions specify specific instances of the variables, i.e. I
3.1 GENERAL SETTING: MULTIVARIATE REGRESSION
Let R2= denote the covariate (conditioning) variable and y Px the multivariate response. Our goal is to construct a prediction set I', (xr C P, such that
TP(y T. (x)))) 6,
3

---

Published as a conference paper at ICLR 2026
with (2) being compact as possible. JAPAN models the conditional distribution p(y [)) using Normalising Flow (NF), which realises the following bijective mapping:
h(y, 2), y = h -=
pz(zz),
where pz a base distribution (typically standard Gaussian).
We compute the log- degsity from this approximation as the conformity score as defined in Eq. (I. After training the NF on the training set, we compute conformity scores on the held-out calibration set. For any candidate test point (d"))), we define the conformal prediction set as:
I re(x") = E e
log x 2 Te),
where Te is the (1 qu-untile of the log-densities observed on the calibration set. Section 3.3 shows the full algorithm. This region satisfies the CP coverage guarantee by construction (Vovk et al. 2023) and, under assumptions such as preservation of rank (Proposition [T) and approximate rank preservation (Proposition 2 the prediction areas closely approximate the prediction areas realised by thresholding on the true unknown density. Proposition[ showws how the inaccuracies in normalising flows estimation impact the efficiency of the computed area.
Proposition (Optimality under Rank Preservation). Let pey x) be the true conditional density, and let fe yy x) be the prodel-estimated density. If there exists a strictly increasing function R - R such that (y xr) = g(ply z)) for all y, then the conformal prediction region defined by fo the same under the true density:
Area(T, (x)) = (r))(
Proposition 2 (Approximate Optimality under Approximate Rank Preservation). Assume foly a) approximates a monotonic transformation of ply x) within a uniform error 8, i.e., ISo(y x) g(p(y x)) < 6 for strictly increasing Then the area difference between the predicted and optimal regions bounded:
(rea(T((x)) (T)) 5(())
where C(ő) + 0 as 0.
3.2 SPECIAL CASE: MULTI-STEP TIME SERIES TRAJECTORIES
JAPAN naturally extends to uncertainty quantification for predicting future of time series trajectories by treating the conditioning variable as a temporal context window. Here, we have access collection of independent time-series trajectories, where each entire trajectory is treated as one data point. All trajectories are assumed to be sampled from the same underlying distribution, meaning that the resulting conformity scores are exchangeable and do not suffer from serial dependence. Note that this setting is different from one where we have access to a single time series and each time index is a data point and exchangeability is violated. Note that this setting yields exchangeable conformity scores and by construction the prediction areas satisfics the coverage guarantees (validity) (English eet al. 2024 Vovk et al. 2023).
Let (f1:++H1 12 be collection of multivariate time series trajectories, where 21.7 RTxD is the observed history and (8)
(4)
+T+1.7+ E HD is the prediction target. We summarise the temporal
context into representation cf computed using recurrent attention iidention-bbsedd encoder such an RNN or Transformer Vaswani et al. (2017) as done in Rasul et al. (2021) Feng et al. (2023). ie.
(t)
C RNN(2 (z,-1.
(8)
(2) (0) 222)), where, of context encoder(z))))
The NF then models the conditional distribution:
i = h( M1THH 2.),
TTITTH = - (2 (0) 20 ~ p2(2),
60
where h is bijective transformation realised by the flow mapping and conditioned on the context vector e The conditional density p(-7)1-T 1:7 | x can evaluated similar to (z; 1: pairs in the multivariate regression case in Eq. (1) and used the conformity score.
4

---

Published as a conference paper at ICLR 2026
3.3 ESTIMATING PREDICTION AREA
Formally, we want to estimate the area of the conformal set
(x)x) == {u E RD logp(u 1)) 2 Te)
where 7, is the density threshold obtained from the calibration scores at level 1 -.
While evaluating coverage in JAPAN is straightforward- requiriin only log-likelihood evaluation and thresholding- the computation of the prediction area is not. One straightforward way can be with methods like Monte Carlo or Quasi-Monte Carlo, but they can be computationally expensive, especially if the support of the label space Y is unknown. Fortunately, due to the bijectíve structure of normalising flows, we can shift the area estimation to the latent space, which allows for efficient computation as in Proposition (see Table for run time comparison against naive MC sampling).
Proposition 3 (Area/Volume computation). Let z pz be samples from the base distribution and let y (z, x) be the realised inverse flow transformation conditioned x, and exp (1,(6,z,))) be area/volume correction under the inverse mapping. Then the area of the prediction set T' (2) can be computed as
exp(o(z,,
Area(T.(x)) N
Σ 1 (2,
Te)
pz(2:)
This estimator accounts for the transformation from latent to data space and includes only those latent samples whose transformed outputs fall within the conformal region. Since the base distribution pz is typically standard Gaussian, can be sampled from efficiently, making this approach practical in moderately high-dimensional settings, See Section B for the proofs. Additionally, the estimator is linked importance sampling (see Section B. which also shows why the estimator is efficient and has low variance.
Algorithm 1 JAPAN: Joint Adaptive Prediction Areas with Normalising -aaaiisingg--wws
Require: Training set D Pirain {Z 9:)] calibration set D {1z,, 8;) 8;) significance level
CE 10)
Ensure: Prediction set I (x) satisfying marginal coverage 2 -
1: Train a conditional normalising flow to estimate p(y x) on Durain:
2: for each calibration example (Tr Mj) E D do
3: Compute conformity score: aj = log p(y; æz) based on Eq.
4: end for
5: Compute threshold т. as the [c (1 - 1 1/m) m mJjth smallest value in {az m
6: At test time: For new input x, define the prediction using Proposition such that:
(x) {u ) log p(v x) 2 T.
4 EXPERIMENTS
We evaluate JAPAN on various datasets (multivariate and time series) and against various baselines with normalising flows used base predictor. In particular, we use CONTRA (Fang et al. 2024
PCP (Wang et al. 2023), RCP (Johnstone & Cox, 2021), NLE (Messoudi et al. 2020). MCQR
(Fang et al. 2024). SPLIT (Izbicki al.J 2020), CQR (Romano al. al.] 2019). COPULACPTS
(Sun & Yu 2024). VANILA
(Messoudi et al. 2021). CRNN (Stankeviciute et al., 2021).
JANET (English et al. 2024) as baselines. All the baselines used either official implementation
public implementations closcly following the set up from the original works. See Table for
qualitative comparison of all the baselines.
In our experiments, we exclusively use discrete normalising flows (Dinh et al. 2017 that is.
flows composed of a finite sequence of explicit invertible transformations such as afhnn (Kingma & Dhariwal 2018: English et al. 2023), spline (Durkan al. 2019). autoregressive Paaamakarios ct al. (2018). These models allow for efficient sampling and exact log- likelihoo evaluation without requiring the numerical integration of ODEs. Moreover, recent works such TaRFl.ow (Zhai
5

---

Published as a conference paper at ICLR 2026
Property
JAPAN
CONTRA
PCP
JANET
Dist-Split
NLE
CQR
MCQR
RCP
CopolaCPTS
VanillaCopula
CFRNN
Disoint
x
X
K
Geometry free
Context-busee
Table I: Comparisons of different baselines, where / Yes; A partially/soft geometric constraint; X = No. JAPAN is the only method that can provide disjoint regions, is context-adaptive, and has no constraint the geometry, whereas PCP and CONTRA have soft constraints on the geometry.
et al. 2024) demonstrated that with sufficiently expressive architectures, their generation quality and likelihood estimation rival those of continuous-time flow models and score-based diffusion models. For discrete models, the conditional density is computed via the change-of-variable formula (Winkler
al. 2019)
(Oh(y, x))
log ply z) = log pz (h(y, z)) + log det
Oy
Time series modelling: While any normalising flow architecture can be used in principle, during modelling it is important incorporate the causal structure inherent time series data. In practice, architectures that preserve the temporal relationship tend to perform better. To this end, we adapt the TARFLOW architecture (Zhai et al. 2024), originally for image data, for time series modelling. Further architectural details are provided in Section C.1
4.1 TOY DENSITIES
To illustrate the behaviour of different conformal prediction methods, we visualise prediction areas on a few synthetic 2D distributions. These datasets have sharp boundaries and disconnected high-density regions, serving as a test-bed for evaluating how well different methods adapt to complex geometries and avoid low- den-iinsity areas.
Figure shows a visualisation of different regions on a spiral density, JAPAN constructs the most compact prediction area, closely following the spiral density and effectively avoiding low-density regions. This reflects its ability to threshold directly on the estimated density in the data space, without imposing geometric constraints. CONTRA, while geometry- spherical regions in the latent space. As a result, low-density regions in the data space that correspond high-density latent samples mistakenly included in the conformal region. PCP also allows for disjoint regions and avoids global geometric constraints, but its local ball-based area construction introduces soft constraints that fail on densities with ring-like modes. In this case, co-eccentric balls cover both dense and non-dense regions indiscriminately. Elliptical and Rectangular conformal methods impose hard geometric assumptions, leading to highly inefficient prediction areas that encapsulate large regions of near-zero density. Computed area along with visualisation of other toy densities can be found in Section
in data space, constructs
4.2 MULTIVARIATE REGRESSION
We evaluate four multivariate regression benchmarks: Energy T'anas & Xifara 2012). RF2D (Spyromitros- Xioufis et al. 2016). RF4D (Spyromitros-Xionn Xixuffs et al. 2016), and SCM20D (Spyromitros-Xion
al. (2016). These tasks span domains from energy efficiency modelling, river flow forecasting, to structured output prediction under supply chain simulation. We repeated the experiments over 25 random data splits. Details the datasets and experiments can be found in Section D and Section
C.2
As shown in Table 2 most methods-including including CQR, MCQR, RCP, DIST-SPLIT, NLE, CFRNN. and VANILLACOPUL. Coheee-achieve comparable marginal coverage across datasets. Figure 2 shows area plots for five methods, with the smallest areas with JAPAN consistently yielding the smallest prediction areas, demonstrating strong volume efficiency without compromising coverage. This advantage is particularly pronounced that JANET and CopulaCPTS perform poorly in some settings, leading to high standard deviations. This is due to these methods using half of the calibration set to train an auxiliary copula or residual predictor model, and their final prediction intervals when the secondary model produces near-zero outputs used as denominators for score rescaling, the resulting prediction regions may explode, leading to highly inflated areas in some runs.
relatively high-dimensional setting RF4D. We also note
become numerically unstable. In particular,

---

Published as a conference paper at ICLR 2026
Methed
Cov. (Knnegyy
Area
Cov. (RF2D)
Area
Cov. (RF4D)
Area
Cov. (SCM)
Area (X10))
CONTRA
028 + 0.03
18.81 7
0.91 +0001
0.35+
0.89 +0001
59.27 + 18.60
0
61.
5.70
PCP
0.88 0.03
16.58 3.87
0.91 +.001
7.39 + 1.02
0.91 +.001
111.63 18.93
0.89 0.01
68.
7.75
NLE
0.57 0.03
9
0.91 +.00
1
0.90 0.00
2732.15
1598.94
0.01
102.1 21.19
30.79
23,86
4175.90
0.01
61.69
MOQR
0.83 0.03
28.99
11.83
0.91 1000
12.07 0.97
0.91 10000
1031 893189911 2279.10
0.91 0.00
70.93 6.85
Dist-ppiit
0.86 a 0.00
35.36 1 12.64
0.91 0.00
19.22 x 1:60
0.91 0.00
1595.92
x: 6355 70
0.91 0:00
87 87.10 9.23
OQR
0.88 + 0.02
$1,12 12.05
0.91 +.00
12.50 1.01
0.91 +.001
1180,5 9536.00
0.91 D.00
84.48 8.74
CogulaCPSS
0.90 0.03
67.20
+000
25.98 4.19
0.91 +001
1591.26 6631
0.89
71.57 +9.09
Copila
0.56 a: 0.03
39.50 30.91
0.91 0.0
25.68
0.89 000
1363.76
5296.00
1
67.84 6.38
CERNN
0 + 11:03
56.22 + 44.29
0.91 +000
27.15 3.31
0.92 +001
$322,47-
14397 00
0.91 0.01
83.60 7.81
JANET
0.88 +000
42.69 + 45.52
0.88 0.02
21,87 13:30
0.85 0.10
1002.3 + 3237.00
0.00000001
64.30 21,06
JAPAN
0.80
1
0.01
5.06
0.91 0.00
24.11 2.33
0.90 0.01
61.28 5.72
Table 2: Coverage and prediction areas for the Energy, RF2D, RF4D, and SCM datasets computed over 25 random data splits. Areas for SCM are shown in units of x10 Almost all methods achieve the desired coverage level of 0.9. We bold the methods for each datasets that have mean coverage of 0.90+0.01 and has the lowest area for the dataset. Thus, JAPAN bolded in all datasets except Drone, but has a smaller volume of all with only a 0.02 difference with the bolded CFRNN method.
Merhade Apea lite Acress Mutthvrrate Datwents
30-
30-
-005
-205
10
NCOE
Figure 2: Prediction areas per dataset (Energy, RF2D, RF4D, SCM) for the top methods computed over 25 random data splits. Each bar's length signifies the standard deviations, with methods sorted by average area for each dataset. JAPAN consistently has the lowest areas and hence is more informative.
4.3 TIME SERIES DATA
We also evaluated on four time series datasets: COVID-19 (Stankeviciute et al. 2021), Particle-1 (Kipf et al. 2018). Drone (Sakai al., 2018), and Pedestrian (van den Berg et al., 2008). each representing different forecasting challenges in scientific or physical systems. These datasets include multi-step prediction tasks where temporal context plays a central role. We repeated the experiment over 25 random data splits. Further details on the dataset and experiments, such as input windows, horizons, and preprocessing, are provided in Section D and Section
As shown in Table 3 all methods achieve similar marginal coverage across datasets. Figure BigureBJshrr shows area plots for five methods, with the smallest areas with JAPAN consistently producing the most compact prediction regions, outperforming all the baselines. Notably, JAPAN yields up to an order- of-magnitude smaller volumes in COVID-19 while maintaining reliable coverage. On the other hand, methods such as CQR, MCQR, and CFRNN tend to generate significantly larger regions. We also observe RCP produce unreliable results in some cases. RCP suffers from numerical instability and undefined scores in some settings, as reflected by missing volume (due to infinity) for COVID- 19 data. Extended Plots for Figure 2 and Figure 3 and the calibration curve aggregated over all datasets for all methods can be found in A
EXTENSIONS: APPLICATIONS OF OTHER DENSITY-BASEL SCORES
Now, we present several extensions that demonstrate the flexibility of using normalising flows for constructing prediction areas. In our main experiments, we used JAPAN as both a predictor and an estimator of log-densities. However, there are practical situations where one may already have a trained predictive model or may wish to use a different base model due to task-specific reasons. In such cases, JAPAN can still be applied on top of the base model and thus decouple the predictive mechanism from the density estimation. This can be viewed using a fixed context encoder, whose representation is learned through an external or unrelated training objective, or could be trained under single objective.
7

---

Published as a conference paper at ICLR 2026
Method
(COVIDDI11
Area
(Particle-1)
(Drone)
Cov. (Pedestrian)
Area
CONTRA
0.92
563.64 98,31
0.87 0.01
2
0.89 = 0.02
1.51 0.11
0.89 1 0.02
0.55 0.03
PCP
0
610.56 1 359,73
1.8.100.0
1.71 1.04
0.89 0.00
1.211
0.88 1 0.00
1.60
NLE
0.90
416.13 57.55
0.87 s:00
2.31 a 0:00
0.88 a 0.02
3.12 0.08
0.88 0.03
2.23
RCP
0.57 0.03
0.88 0.02
2 0.01
0.88 0000
3.16
0.87 0.02
2.00
MCQR
2.03
1276.50 222.80
2.89
0.86
0.12
0.86 0.02
Sisi-Splii
0.85 0.03
652.66 80.66
0.96 0.0]
0.93 0.0
4,47 0.228
0.9 0.00
2.20 0.07
CQR
0.9 0.02
1048 60 175.6%
0.96 0.01
3.33++0.06
0.98 + 0.03
5.12 0.74
0.96 0.02
2.29 +.15
CopulaCPTS
0,85 = 0.077
627 72 167.06
0.80+0000
4.12+0.12
0.92 = 0.05
4,46 0.73
0.91 +.00
2.19 D.26
ConilaCopula
0.04
539.09 1: 79.93
0.87 0.00
2.91 0.05
0.86
3.48 0.08
0.87 0.02
1,91 0.05
CERNN
0.92=0003
112,44
0.09
0.99
0.76
JANET
0.88
669 + 364.47 47
0:92 0.00
2.78 0.09
0.90
0.26
0.92 0.00
2.20
JAPAN
0.91 0.00
400.94 74.39
0.91 0.83
0.89 + 0.02
0.88 = 0.02
1.47+00111
0.89 + 0.02
0.50 0.04
Table Coverage and prediction areas on the COVID-19, Particle-1, Drone, and Pedestrian time series datasets computed over 25 random data splits. Almost all methods achieve the desired coverage level of 0.9. We bold the methods for each datasets that have a mean area of 0.90+0.01 and has the lowest area for the dataset. Thus, JAPAN is bolded in all datasets except Drone, but has a smaller volume of all with only a 0.01 difference with the bolded CONTRA method.
Notrodd femtili
COVC-43
900
14-
A99
112
530
2.0
Figure 3: Prediction areas per dataset (COVID, Particle-1, Drone, Pedestrian) for the top 5 methods computed over 25 random data splits. Each bar's length signifies the standard deviations, with methods sorted by average area for each dataset. JAPAN consistently has the lowest areas and hence more informative.
We explore several such applications, namely unconditional modelling (Section A. where the flow is trained independently of input features or the context to learn p(t); conditional modelling (Section A.3), where the flow is conditioned on representation from a pretrained base model to learn p(ylū): posterior modelling (Section A.4). where the flow is trained to learn the posterior predictive distribution (p(ğly) or more generally, p(rlyy, instead of p(x[y)); latent modelling (Section A.5). where we use density of the transformed sample in the latent space, thus making CONTRA (Fang et al. 2024) special case of JAPAN; adaptive т. modelling (Section AA6), where we discuss adapting based on the condition and effectively getting closer to conditional coverage under mild assumptions. Here, we briefly introduce conditional modelling and posterior modelling and defer to Section A for detailed discussion on all the extensions, and to Section A.7 for discussion on the assumption violations with underfitted models and the corresponding ablation experiments.
5.1 CONDITIONAL MODELLING
In this setting, we assume access base predictive model jı that outputs point predictions و = p(z). Then, we train a conditional normalising flow to model P(Y f). This can be seen as obtaining a fixed representation extracted from fixed base model and training a flow separately from the base predictor, in contrast to standard end-to-end JAPAN.
This setting creates connection between the predictions of the base model and the true targets y. allowing the resulting density to capture input-dependent uncertainty. While the addition of another model introduces potential estimation error, we find that the regions constructed via log- likelihood thresholding remain both compact and valid as the flow model tries to correct for the base predictor's 's distribution estimation error. Consequently, this form of conditional modelling achieves more compact regions than those obtained from other baselines. Please refer to Section A.3 for experimental analysis.
8

---

Published as a conference paper at ICLR 2026
5.2 POSTERIOR MODELLING
In this setup, we aim to model the posterior distribution over input conditions given an output, i.e., P(X Y) or as a special case P(Y Y), using normalising flow. We then use the resulting log-pon log-posterior density as a conformity score. While the general framework remains similar, this formulation allows to test different possible input conditions for given output Y, and select those that are most plausible under the learned posterior.
Posterior modelling has several practical advantages. For instance, in settings where the response variable is categorical or discrete, it is often difficult or ill- ill-posed to model P(Y | XX) directly using normalising flows. Modelling P(X Y), however, is tractable and aligns with standard generative modelling pipelines- for example, generating images conditioned on class labels. In such cases, one can perform conformal-based calibration posterior log-densities to produce valid prediction sets. This approach can be used both to assess how plausible a generated output is and to identify whether an observed outcome falls within calibrated prediction set.
While it is known that using the conditional density p(y x:) with the threshold T, adaptive on the condition can yield the smallest conformal region under ideal conditions, it is less clear how tight the regions are when the conformity score is based on the posterior density p(z y). Nevertheless, we find that this approach still produces compact and valid regions in practice, albeit typically less efficient than those obtained via conditional likelihood modelling. Please refer to Section A.4 for experimental analysis.
6 RELATED WORK
The literature on conformal prediction (CP) (Gammerman et al. 1998 Vovk et al. 2005; Angelopou
los & Bates| 2022) is extensive, with works on new (non- )conformity scores (Sun & Yu 2024
English et al. 2024), going beyond exchangeability (Chermozhukov et al., 2018: Barber al. 2023
etc. In the context of non-univariate response variables, CP research primarily branches into three directions: (i) multivariate regression with i.i.d. data, (ii) forecasting in single long time series with long- tinghhorioon or multivariate time steps (where temporal dependence violates exchangeability of
conformity scores) (Gibbs & Candes 202T Chernozhukov et al. 2018 Lee et al. 2025), and (iii) settings where exchangeable set of time series trajectories are available allowing CP to be trivially applied.
In the multivariate regression setting, RCP (Johnstone & Cox, 2021 adapts to multivariate depen- dencies using covariance structure to form elliptical areas. Similar RCP, NLE (Messoudi et al.. 2020) creates elliptical areas but is context adaptive. PCP (Wang et al. (023) leverages generative models to sample new conditional outputs and then defines prediction sets via unions of local L2 balls, with (Plassier et al.) 2024) extending to get approximate conditional validity. While these avoid global shape constraints, still relies on local spherical balls based on residual distances, which refer to soft geometric constraint. CONTRA (Fang et al. 2024). more similar to our work, employs normalising flows but constructs conformal balls in the latent space that we again term soft constraint. This allows some flexibility but risks including low-density areas in the output space. Another work (Izbicki et al., 2021) provides a disjoint region but requires a lot of hyperparameters and can have stability issues. Similarly, ST- (Feldman al., (023) uses conditional variational autoencoders to learn low-dimensional latent representations but does not work directly with densities. Other related methods originally for univariate response but adapted for multi-variate data by Fang et al. (2024) include DIST-SPLIT [Izbicki et al. 2020) (which estimates CDFs to construct areas), LOCALLY ADAPTIVE SPLIT CP (which introduces local calibration via splitting), and quantile-based methods like CQR (Romano et al., 2019) and MCQR (Fang et al. 2024]. which are well-calibrated but limited in their ability to handle multimodality. Another work, Colombo (2024) also uses Normalising Flows in the context of conformal prediction, but they used it to train the calibration process for univariate response, i.c. on the distribution of the residuals, only flow model to construct a region.
whereas
For time series trajectories, many conformal approaches focus on structured prediction over multiple horizons. BJ-RNN (Alaa & van der Schaar 2020) applies the jackknife+ method Barber al. (2020) to recurrent networks, but its theoretical coverage guarantee is 1 2e, which is looser than the standard CN-RNN (Stankeviciute al.) 2021) assumes conditionally i.i.d. residuals and

---

Published as a conference paper at ICLR 2026
applies Bonferroni correction (unnn 1961), which become overly conservative horizon
length increases. (Cleaveland et al., 2024) proposed an LP-based approach for multi-step error
modelling, but it requires an additional dataset and extensive parameter tuning. COPULACPTS (Sun & Yu 2024) builds on Messoudi et al. (2021) to use copulas for temporal dependency modelling calibration sets and gradient-based training,
in multivariate time series. However, it requires introducing inefficiency. Recent methods like CAFHT (Zhou et al. 2024) and JANET (English al. 2024) construct rectangular areas that are adaptive the input context but cannot produce disjoint or multimodal areas. By contrast, JAPAN's use of normalising flows enables areas that are both context context-adaptive and geometrically unconstrained, yielding compact, expressive prediction areas across both multivariate and time series settings.
recent work (Lee et al., 2025 uses flow-matching but works in the setting where only a single time series is provided and tackles the problem of exchangeability violation, whereas JAPAN tackles a completely different problem with different approach. CCN-JPR (English & Lippert (024) introduced discrete normalising flows for time-series with empirical coverage guarantees, but without area computations. JAPAN extends CCN-JPR by (i) providing a comprehensive analysis across multivariate regression and time-series settings, (ii) establishing theoretical results, (iii) proposing novel conformity scores, and (iv) developing an efficient method for area computation.
7 DISCUSSION AND LIMITATIONS
Our experiments demonstrate that JAPAN flexible and effective framework for constructing prediction regions using normalising flows. It can be applied in a wide range of settings, including multivariate regression and time series forecasting, and can be seamlessly integrated with existing predictive models. Across tasks, we observe that JAPAN consistently achieves valid coverage while producing significantly smaller prediction regions than baseline methods, particularly in high-dimensional settings.
Despite these promising results, several limitations remain. In particular, while we discuss conditional coverage when using adaptive Te extension (see Section A.6).. JAPAN-like all methods based on marginal conformal prediction does not guarantee conditional coverage. That is, while the prediction set satisfies IP(y E T2 (x)) 2 may vary significantly between different subpopulations or contexts. Finally, while our approach can be extended to use predictions from another model, posterior densities, or adaptive threshold, each step introduces additional assumptions and potential estimation errors. Future work may investigate the theoretical benefits these alternatives offer.
in expectation over the distribution of æ. the coverage
Reproducibility Statement: All experimental results reported in this paper can be reproduced using the implementation details provided in Section C Dataset descriptions and preprocessing steps are detailed in Section Complete proofs of all theorems and propositions are included in Section B
This work does not raise any specific ethical concerns. All datasets used are publicly available, and no human or animal subjects were involved. We have adhered to the ICLR Code of Ethics throughout the preparation of this paper.
Ethics Statement:
ACKNOWLEDGMENTS
We extend our gratitude to Prof Vladimir Vovk for helpful discussions and to Thomas Gaertner, Juliana Schneider, Noel Danz, and Wei-Cheng Lai for their feedback on an earlier draft. This research was supported by the HPI Research School for Digital Health. We gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - project number 459422098.
REFERENCES
Ahmed M. Alaa and Mihaela van der Schaar. Frequentist uncertainty in recurrent neural networks via blockwise influence functions, 2020.
10

---

Published as a conference paper at ICLR 2026
Michael Samuel Albergo and Eric Vanden- Building Normalizing Flows with Stochastic Inter- polants. September 2022. URL httpps://ppenreviieii net 1tt.://..nntt/toooomidd-d---
Anastasios N. Angelopoulos and Stephen Bates. A gentle introduction to conformal prediction and distribution-free uncertainty quantification, 2022.
Rina Foygel Barber, Emmanuel Candes, Aaditya Ramdas, and Ryan J. Tibshirani. Predictive inference with the jackknife+,, 2020.
Rina Foygel Barber, Emmanuel J. Candes, Aaditya Ramdas, and Ryan J. Tibshirani. Conformal prediction beyond exchangeability. 2023.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud. Neu- ral ordinary differential equations. Advances in neural information processing sys- tems, 31. 2018. URL
neurips cC ct/ cctpapper//01/haas.. (a8.8ffffff666222226333359596-33
Victor Chemoozuuko, Kaspar Wüthrich, and Zhu Yinchu. Exact and robust conformal inference methods for predictive machine learning with dependent data. In Sébastien Bubeck. Vianney Perchet, and Philippe Rigollet (eds.), Proceedings of the 31 Conference On Learning Theory, PML.RMR.. PMLR, 06-09 Jul 2018. mlr wts.
volume 75 of Proceedings of Machine Learning Research, pp. 732 749.
URL
Matthew Cleaveland, Insup Lee, George J. Pappas, and Lars Lindemann. Conformal prediction regions for time series using linear complementarity programming. 2024.
Nicolo Colombo. Normalizing flows for conformal regression, 2024. URL https: //arxivv
org/abs/ 2406. 03346
Laurent Dinh, Jascha Sohl- Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp, 2017.
Olive Jean Dunn. Multiple Comparisons among Means. Journal of the American Statistical Asso- ciation, 56(293):52-66 64, 1961. ISSN 0162-1459, 1537-274X. doi: 10. 00.180000121455.1199 10482090. URL http:/ http://wwwww.dd..ooo../eo Landfoninne com/ dtt//bs/1/.a.ooeooooo 100.0000000121000
1961 10482090
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows, 2019.
Eshant English and Christoph Lippert. Conformalised conditional normalising flows for joint prediction regions in time series. In NeurIPS 2024 Workshop on Bayesian Decision-making and Uncertainty, 2024.
Eshant English, Matthias Kirchler, and Christoph Lippert. Kernelised normalising flows, 2023.
Eshant English, Eliot Wong- Toi, Matteo Fontana, Stephan Mandt, Padhraic Smyth, and Christoph Lippert. Janet: Joint adaptive prediction-region estimation for time-series, 2024. URL https
Zhenhan Fang, Aixin Tan, and Jian Huang. CONTRA: Conformal Prediction Region via Normal- izing Flow Transformation. October 2024. URL httpps//pppeeeeveeeweeerrrorr
p009cqLq70
Shai Feldman, Stephen Bates, and Yaniv Romano. Calibrated Multiple-Output Quantile Regression with Representation Learning. Journal of Machine Learning Research, 24(24):1-48 2023. ISSN
1533-7928. TRRhhttp:// jmlr org/ wwtppps/ww//1-12222222 html
Shibo Feng, Chunyan Miao, Ke Xu, Jiaxiang Wu, Pengcheng Wu, Yang Zhang, and Peilin Zhao. Multi-scale attention flow for probabilistic time series forecasting, 2023.
A. Gammerman, V. Vovk, and V. Vapnik. Learning by transduction. In Proceedings of the Fourteenth Conference on Uncertainty in Artificial Intelligence, UAI'98, pp. 148-155, San Francisco, CA. USA, 1998. Morgan Kaufmann Publishers Inc. ISBN 155860555X.
Isaac Gibbs and Emmanuel Candès. Adaptive conformal inference under distribution shift, 2021.

---

Published as a conference paper at ICLR 2026
Will Grathwohl, Ricky Т. Q. Chen, Jesse Bettencourt, llya Sutskever, and David Duvenaud. FFJORD: Free-Form Contimuous Dynamics for Scalable Reversible Generative Models. September 2018.
URL httptppss//ppeeeveee...oo http:t//oooooonodddddddaooooaaa
Rafael Izbicki, Gilson Shimizu, and Rafael Stern. Flexible distribution-free conditional predictive bands using density estimators. In Proceedings of the Twenty Third International Conference
on Artificial Intelligence and Statistics, pp. 3068-3077. PMLR, June 2020. URL https:/7
proccedinns.mlr... http://w..........eeee///////a.a
Rafael Izbicki, Gilson Shimizu, and Rafael B. Stern. Cd-split and hpd-split: efficient conformal
regions in high dimensions, 2021.
Chancellor Johnstone and Bruce Cox. Conformal uncertainty sets for robust optimization. In Proceedings of the Tenth Symposium Conformal and Probabilisticc Prediction and Applications, pp. 72-90. PMLR, September 2021. URL http:::/ppoceeedinn/ .twp:eeeiiiiiiirm.pppo..e.m...ee mlr p11essssssee
Johnstone2la, html
Diederik P. Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1xl convolutions, 2018.
Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems, 2018.
Junghwan Lee, Chen Xu, and Yao Xie. Flow-based conformal prediction for multi-dimensional time series, 2025. URL https: arxiv. www.ooggg//b....... 25w. 05709
Jing Lei, James Robins, and Larry Wasserman. Distribution Free Prediction Sets. Journal of the American Statistical Association, 108(501):288-287, 2013, ISSN 0162-1459. doi: 1080/
Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow
Matching for Generative Modeling. September 2022. URL
Forum? BerrummLd--ogWWAcCCT9
Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Conformal multi-target regression using neural networks. Alexander Gammerman, Vladimir Vovk, Zhiyuan Luo, Evgueni Smirnov, and Giovanni Cherubin (eds.), Proceedings of the Ninth Symposium on Conformal and Proba- bilistic Prediction and Applications, volume 128 of Proceedings of Machine Learning Research, wtt.ooss.......rseeeem.m
pp. 65-83, PMLR, 09-11 Sep 2020. URL htttts://proooeeeo
messsodd20aa html
Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Copula-based conformal predic- tion for multi-target regression, 2021.
Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. Machine learning: ECML 2002: 13th European conference on machine learning Helsinki, Finland, August 19-23, 2002 proceedings 13, pp. 345-356. Springer, 2002.
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation, 2018.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference, 2021.
Vincent Plassier, Alexander Fishkov, Mohsen Guizani, Maxim Panoy, and Eric Moulines. Probabilistic Conformal Prediction with Approximate Conditional Validity. October 2024. URL https
Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann, and Roland Vollgraf. Multi- variate probabilistic time series forecasting via conditioned normalizing flows, 2021.
Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized Quantile Regres- sion. In Advances in Neural Information Processing Systems, volume 32. Curran Asso- ciates, Inc., 2019. URL https:: http:tppppppppcceeedigieeeeeeeee wwwpuuupss.cpppppe.eeppp html
12

---

Published as a conference paper at ICLR 2026
Atsushi Sakai, Daniel Ingram, Joseph Dinius, Karan Chawla, Antonin Raffin, and Alexis Paques. Pythonrobotics: a python code collection of robotics algorithms, 2018.
Glenn Shafer and Vladimir Vovk. A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3), 2008.
Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Based Generative Modeling through Stochastic Differential Equations. October 2020.
URL/httppss/ppperevviie wwwp://wwwwww..eeooooooooaaaaaae http:////waaageeeeoooocooooooooa
Yang Song, Conor Durkan, lain Murray, and Stefano Ermon.
Maximum likelihood training of score-based diffusion models. Advances in neural information processing systems, 34:
1415-1428, 2021. URL
0afaabb11ebccoorrrr0008222222222 .hmml
Eleftherios Spyromitros-Xions Grigorios Tsoumakas, William Groves, and loannis Vlahavas. Multi-Target Regression via Input Space Expansion: Treating Targets as Inputs. Machine Learning, 104(1):55-98, July 2016. ISSN 0885-6125, 1573-0565. doi: 10. 10000//11994--1--5555666 URL
http://arxiv> org/ abs/ a211 6581
arXiv:1211.658) [cs].
Kamile Stankeviciute, Ahmed M. Alaa, and Mihaela van der Schaar. Conformal Time-series Fore- casting. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P. S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 6216-6228. Curran Asso- ciates, Inc., 2021. URL http:ttt/pocceeeeeeeee wtw.e/iiiinnsueeuueeeeeeeee nwwrips.up.... wwt.:eewwwwrr....eeeeeeeeeeeeeee cc/paper files/paper/ 2011/fLee/12111ae222223eeee
Paper.pdi
Sophia Huiwen Sun and Rose Yu. Copula conformal prediction for multi-step time series prediction.
In The Twelfth International Conference on Learning Representations, 2024. URL https
[
net./ (aarrummd==ojjIJaa
Athanasios Tsanas and Angeliki Xifara. Accurate quantitative estimation of energy performance of
residential buildings using statistical machine learning tools. Energy and Buildings, 49:560-567,
June 2012. ISSN 03787788 doi: 10.1016/j enbuild..0012.033000. URL
elsevier com http:////eeeeereeeeea/////131110 ptt.eeaeee..cccc//aaaaaaaaaaaaaa
Jur van den Berg, Ming Lin, and Dinesh Manocha. Reciprocal Velocity Obstacles for real-time multi-agent navigation. In 2008 IEEE International Conference Robotics and Automation, pp. 1928-1935, May 2008. doi: 10. 10000//00....000055533000 URL
leee org/ wow.meett/4443444
ISSN: 1050-4729.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, E ukasz Kaiser, and Illia Polosukhin. Attention is All you Need. In Ad- vances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL htttps::/prooceeeiiii. .ew.eeeeeeddnns.eeenn.s.... cttpcccpppeee ppwper//01//// Abst html
V. Vovk, A. Gammerman, and G. Shafer. Algorithmic Learning in Random World. Springer International Publishing, 2023. ISBN 9783031066511. URL hhttps::/bbbttt googe google
Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning in random world, volume 29, Springer, 2005.
Zhendong Wang, Ruijiang Gao, Mingzhang Yin, Mingyuan Zhou, and David Blei. Probabilistic Con- formal Prediction Using Conditional Random Samples. In Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, pp. 8814-8836. PMLR, April 2023. URL
https::/proceeei.ngg mlr. www.mmm...........
press v20 / wwn.///ag222 2aanng2nn.. 2n. html
Christina Winkler, Daniel Worrall, Emiel Hoogeboom, and Max Welling. Learning Likelihoods with Conditional Normalizing Flows. September 2019. URL https: httttp::///peeeeeiew.. net/ forum? d
13

---

Published as a conference paper at ICLR 2026
Shuangfei Zhai, Ruixiang Zhang, Prectum Nakkiran, David Berthelot, Jiatao Gu, Huangjie Zheng, Tanrong Chen, Miguel Angel Bautista, Navdeep Jaitly, and Josh Susskind. Normalizing Flows are Capable Generative Models, December 2024. http:tts:: //ui adsabs hwwvard. edu/ abs 3bas/2022arrxiv22222662222
ADS Bibcode: 2024419912412066299
Yanfei Zhou, Lars Lindemann, and Matteo Sesia. Conformalized adaptive forecasting of heteroge- neous trajectories, 2024.
14

---

Published as a conference paper at ICLR 2026
A EXTENSIONS AND ADDITIONAL EXPERIMENTS
A.1 TOY DENSITIES
To illustrate the behaviour of different conformal prediction methods, we evaluate them on a set of synthetic 2D distributions: Spiral, Checkerboard, Circles, and Moons. These toy densities are multimodal and feature sharp discontinuities or disconnected high-density regions, making them strong testbed for evaluating geometric flexibility of conformal methods.
Figure 4 and Figure show visualisations of the conformal prediction regions produced by different methods on these datasets. JAPAN consistently produces compact regions that accurately track the support of the underlying distribution, avoiding low-density areas and capturing disconnected modes. By contrast, hard-constrained methods like Elliptical and Rectangular conformal prediction yield inefficient regions that fail to capture the data geometry. CONTRA, while geometry- free in the data space, forms spherical regions in the latent space; these often map to big continuous shapes in the data space that may include large, low-density regions. PCP introduces soft geometric constraints by creating unions of balls around generated samples. While this setup allows for disconnected regions, vacuous space when the density is highly non-uniform..
the local ball structure still
Quantitative metrics including coverage and area at c = 0.1 are reported in Table 4 and area comparison plotted in Figure JAPAN consistently achieves the lowest areas while maintaining valid coverage, Additionally, calibration curves across a range of e E [0.05, 0.95] are shown in Figures 10 and confirming that JAPAN remains well-calibrated across different significance levels.
Density
JARAN region
rogion
region
region
Rectangular reglon
JAPAN region
region
region
region
Rectanguler reglon
region
region
region
region
Rectenguler region
Figure 4: Visualisation of different densities (Moons, Rings, and Checkerboard) at 0.1. Only JAPAN captures the modes and avoids low-density regions. Hard constrained areas, such as Elliptical and Rectangular, naturally fail capturing the true shape. CONTRA encapsulates low-density regions, completely ignoring discontinuities in the density. PCP while captures discontinuity includes some low-density regions due to the soft constraints.
A.2 UNCONDITIONAL MODELLING
In the unconditional modelling setting, we assume access base predictive model μ that produces point predictions = pzz). Rather than modelling the conditional density p(y x) directly, we instead fit normalising flow to estimate the marginal distribution IP(Y') over predicted outputs on the training set. Once the flow trained, we apply JAPAN with a slight modification: for each calibration example (Ir y), evaluate the log-density log Py (y) under the learned marginal model (rather than conditioning on æ) and use this as the conformity score (see Algorithm 2). This change ensures that calibration is still performed with respect the true output y-
15

---

Published as a conference paper at ICLR 2026
Method
(Checkerbourd)
(Moens)
(Cireles)
Cov. (Spiral)
Area
CONTRA
0.91 0.02
41.01 0.455
0.91 1 0.01
1
0.91 1 0.02
34.06 0.10
0.91 1 1.00
25.21 0.01
PCP
0.89.00.00
32.62.10.33
0.89 0.01
2
0,90 0.01
18.43 2,09
0 0.00
14.28 0.01
4.83
RCP
0.91 1:001
57,50 0.61
0.89 0.00
1.84
0.90 0.02
3
0.50 0.01
26.38 2.64
MCQR
0.90
57.48 0.66
0.89
4.52
0.90 0.01
35.53
0.90 0.02
31151
Dist-Split
0.91
57.67 0.45
4.70 0.09
0.90 0.01
a.44
0
32.64
CQR
57.77 0.21
0.01
32.
JANET
0.90
57.59 0.31
0.89 = 0.02
4.57
0.90 0.01
35.57
0.90 0.01
31.57 3.16
Capila
+.90
57.17+0.33 0.336
0.89 + 0.01
4.37 0.00
0.90 + 0.01
35.39
0.90 0.02
32.11 2.90 2.96
CERNN
0.01
39.11 0.52
0.02
4.92 0.19
0.01
37.49
0.90 0.01
34.24 3.22
JAPAN
0.90
27.95 0.04
0.90 0.01
2.06 0.13
0.91
14.76 0.67
0.90
11.65 1.17
CopulaCPTS
0.91 0.02
57.74 0.002
0.89 0.01
5.26 + 0.64
0.92 0.01
35.74 1.10
0.90 + 0.01
32.89 3.29
Table Coverage and prediction areas for the Checkerboard, Moons, Circles, and Spiral datasets computed over 25 random data splits. Almost all methods achieve the desired coverage level of 0.9. We bold the methods for each datasets that have a mean coverage of 0.90+0.01 and has the lowest area for the dataset. Thus, JAPAN is bolded in all datasets.
Area Plots Different Methads Across Toy Datasets
Checkerboard
Moons
5.5
50-
2.-
POP
Gircles
Spiral
22
35.
20
15
14
creas
Figure 5: Computed areas per dataset (Checkerboard. Moons, Circles, Spiral) for all the methods computed 25 random data splits. Each bar's length signifies the standard deviations, with methods sorted by average area for each dataset. JAPAN consistently has the lowest areas and hence more informative.
A key limitation of this approach is that the resulting prediction areas are not adaptive to the input context. Since the flow models global distribution over outputs, all test-time regions are drawn from the same marginal. Nonetheless, we find that the method can produce surprisingly compact and valid prediction regions, especially when compared baselines that rely on fixed geometric constraints.
We evaluate this setup on two tabular datasets: Energy and RF4D. Results are shown in Table 5 On the Energy dataset, JAPAN achieves the most compact prediction areas among all methods that satisfy the target coverage range of 0.90 + 0.01. On RF4D, area efficiency is less consistent: although JAPAN remains competitive, the average area skewed upward due to a single run with unusually large volume. This effect is even more pronounced in baselines like CONTRA and PCP, which also use latent-space conformity but without density thresholding. In contrast to JAPAN, these methods produce highly inflated regions due to poorly calibrated or geometrically mismatched latent representations. Meanwhile, standard baselines such as CQR, MCQR, also performed worse than JAPAN, despite not having the additional flow modelling. Table summarises the coverage and
16

---

Published as a conference paper at ICLR 2026
Area Plots For Different Methods Across Multivariate Datasets Energy
120
35
30
25
20
15
10
PCF
CONTRA
MIUE
NoSe
RCPP
immiaccppua
JAMET
Capuscopts
CONTRA
Pgpp
MCOR
NiF
Spir
JaNrr
RCP
ConuliCTE
RF4D
SCM
15300
10200
5300
70000
-5009
60000
50000
-10209
40000
CONTRA
PCR
JAMEY
MCCR
Variiic Ca
Copstacin
CreNN
conme
PCPP
MCON
CopolacPiyy
CFRIN
Car
Dist-jiit
RR.E
Dist
Figure 6: Prediction areas per dataset (Energy. RF2D, RF4D, SCM) for all the methods computed over 25 random data splits. Each bar's 's length signifies the standard deviations, with methods sorted by average area for each dataset. JAPAN consistently has the lowest areas and hence is more informative.
Area Plots For Different Methods Across Time Series Datasets
COVID-19
Particle
2.5
800
500
409
1.0
203
ManiacCypula
canrr
PGPP
CopularP7
Datt spali
HisE7
OPAW
NCRR
CONTRE
MOET
Mcoi
Santaraamaaa
Bist Reltt
cail
Drone
Pedestrian
2.5
3.5
CONTRe
RCE
RCP
MCDW
JUKTT
ConsacoPS
CDN
c'een
MRR
coNTR
POP
McON
tasstacaeste
CoPUACPPS
JANET
ALEE
col
crecn
Figure 7: Prediction areas per dataset (COVID, Particle-1, Drone, Pedestrian) for all the methods computed over 25 random data splits. Each bar's length signifies the standard deviations, with methods sorted by average area for each dataset. JAPAN consistently has the lowest areas and hence more informative.
area metrics computed over 5 random splits. Methods with the most compact regions among those achieving near-nominal coverage are highlighted.
17

---

Published as a conference paper at ICLR 2026
Method
Cov. (Energy)
Area (Energy)
Cov. (RF4D)
Area (RF4D) x 10
CONTRA
0.90 + 0.02
2794.48 + 2420.62
0,90 + 0.01
933333.20 + 1229039.0
PCP
0.89 = 0.01
1612.79 2192.88
0.89 + 0.00
486646.40 = 389587.7
NLE
0.87 = 0.03
129.10 38.56
0.90 + 0.01
0.37 = 0.39
RCP
0.87 = 0.02
124.46 29.45
0.00
0.77 = 0.87
MCQR
0.88 + 0.03
1201.33 = 802.06
0.00
9149.48 + 22.60
Dist-Split
0.55 0.21
4
0.02
0.67 = 1.32
CQR
0.91 = 0.02
2093.45 + 2367.83
0.90
1096.55 + 4508.45
JAPAN
0.90 0.02
770.56474.5
0.00
25.86 + 31.53
CopulaCPTS
0.90
252.95 # 67.94
0.01
11,41 + 4.80
VanilaCopula
0.86 1566
175.12 + 48.13
0.89
0.74 + 0.09
CFRNN
0.88 = 0.02
207.45 33.69
0.92
1.34 + 0.27
JANET
0.87 = 0.04
207.44 64.69
0.88 0.01
6.37 + 0.75
Table 5: For Unconditional modelling: Coverage and prediction volume for the Energy and RF4D datasets at e = 0.10. RF4D volumes are rescaled by 10 for clarity. Methods with the lowest volume among those achieving coverage in [0.89, 0.91] are bolded.
Algorithm 2 JAPAN with Unconditional Normalising Flows
Require: Training set Duain {2,, W) calibration set Deal {12g, 8/)] significance level C 10
Ensure: Prediction set I. satisfying marginal coverage 2 1
1: Fit normalising flow model to the marginal distribution of predictions ÿ; =xxx) using (ú:)
2: for each calibration pair (zj, E Deal do
3: Compute conformity score: aj = log (Vú)
4: end for
5: Compute threshold as the k-th smallest value in {az} }リ=1 where k = (1 - 1/m) ml
6: At test time: For any new input x, the prediction set using Proposition 3 |uchh that:
re = {yE log (v) 2 T
A.3 CONDITIONAL MODELLING
We evaluate the conditional modelling extension of JAPAN on two tabular datasets: Energy and RF4D. In this setup, we first train a base predictive model μ to output point estimates á uxx), and then train a conditional normalising flow to model pvy 9). At inference, prediction regions are constructed by thresholding the log-likelihood of y under the conditional flow (see Algorithm B).
Table summarisee the coverage and areas for all methods. JAPAN consistently achieves competitive or sowest-volume prediction sets while maintaining valid coverage. On the Energy dataset, JAPAN achieves the most compact regions among all methods that fall within the target coverage range of 0.90 + 0.01, outperforming all baselines efficiency. In contrast, methods such as CQR, JANET, and CopulaCPTS yield significantly larger regions despite comparable coverage.
On the more challenging RF4D dataset, the benefits of conditional modelling are even more evident. JAPAN achieves volume of just 0.69 x 10 nearly order of magnitude smaller than most baselines, including residual-based and copula-based approaches. The flow model trained on p(y 9) able to adaptively reshape the predictive density correct for errors in the base predictor, yiclding highly compact and calibrated prediction areas. Notably, flow-free basclines such as MCQR, CFRNN, and CopulaCPTS result in excessive volumes or instability (e.g., MCQR reports region volumes the order of 17)
These results highlight the effectiveness of using JAPAN after base model. Even though the flow trained separately from the base model, its ability to correct for systematic biases ý enables significantly efficient conformal sets than any competing method.
18

---

Published as a conference paper at ICLR 2026
Method
Cov. (Energy)
Area (Energy)
Cov. (RF4D)
Area (RF4D) x 10
CONTRA
0.89 + 0.03
53.87 + 14.34
0.88 + 0.01
2.28 + 0.67
PCP
0,89 = 0.02
176.66 + 163.27
0.89 # 0.01
8,94 + 1.05
NLE
0.87 + 0.04
82.75 9.49
0
91.46 + 5.35
RCP
0.87 0.02
124.46 29,45
0.90
77.15 # 8.72
MCQR
0.87
212.58 145.09
0.00
2.71e+17 3.22e++17
Dist-Split
0
151.26 56.65
0.00
20.69 = 3.14
CQR
0.90 0.02
232.93 + 130.44
0.00
20.75 2.67
JAPAN
0.90
50.10 15.36
0.89
0.69 0.04
CopulaCPTS
0.89
2 69.06
0.00
109.84 28.84
VanilaCopula
0.04
175.12 + 48.13
0
74.11 9,46
CFRNN
0.88 + 0.02
207.45 # 33.69
0.92 # 0.01
134.10 # 27.16
JANET
0.90
185.62 38.79
0.89 = 0.00
97.41 28.28
JAPAN Posterior
0.88 + 0.02
60.85 + 5.23
0.89 + 0.01
67.87 16.40
Table 6: For Conditional modelling: Coverage and prediction volume for the Energy and RF4D datasets at 0.10. RF4D volumes are rescaled by 10s for clarity. Methods with the lowest volume among those achieving coverage in [0.89, 0.91] are bolded. JAPAN has the lowest volume with the likelihood-brae -asedd threshold, whereas JAPAN with the posterior threshold fares better than con-flow-based baselines.
Algorithm JAPAN with Conditional Normalising Flow on Predicted Inputs
Require: Training set Dirain {(r, y=) calibration set Deal = {(0,, 95) }= significance level
C (0, 1)
Ensure: Prediction set T. (x) satisfying marginal coverage 2 1
1: Train a base model μ on Diain to predict ú μ(x;))
2: Train conditional normalising flow on pairs (Ūs, y:) to estimate p(y
3: for each calibration point (xzz Uj) E Daal do
4: Compute új (u(x;)
5: Compute conformity score: aj = log p(y; (5)
6: end for
7: Compute threshold as the k-th smallest value in {az]ī= where k = [e (1- 1/m) (mmm)-7 ml
8: At test time: For any new input x; define the prediction using Proposition 3 such that:
r (x) {y E y log u(x)) 2 re)
A.4 POSTERIOR MODELLING
We evaluate the posterior modelling extension of JAPAN the Energy and RF4D datasets, using normalising flow trained to model the inverse conditional density p(0 y). where ğ = p(x) is the output of fixed base predictor. The posterior log-density used a conformity score construct the calibrated prediction set as shown in Algorithm 4
As shown in Table this approach achieves valid coverage on both datasets and produces reasonably compact prediction regions. On Energy, JAPAN Posterior yields volume of 60.85, which is competitive with likelihood-basee JAPAN and smaller than most non-flow baselines. On RF4D, however, the posterior-based method underperforms significantly against its likelihood counterpart, with a volume of 67.87 compared to 0.69 (rescaled), but remains more effective than most of the non-flow baselines, except for CQR and has high coverage.
Intuitively, posterior-based conformity can be interpreted as testing how plausible an input (or a prediction ú) is given candidate output y. This mirrors likelihood- scoring in spirit: true (x, y) pairs are expected to have high posterior values, and low www-potterior values are assigned to implausible combinations. However, posterior calibration may not be practical in saciiiee-paticularrlyy when the input dimensionality is high, and sampling becomes expensive (c.g., for area estimation), or when the posterior distribution p(æ; y) is more complex or less structured than the forward conditional
19

---

Published as a conference paper at ICLR 2026
ply =)) In contrast, likelihood- based conformity scores allow efficient estimation via latent- space sampling. (see Proposition 3
Nevertheless, re-emphasise that posterior-based JAPAN is promising alternative in settings where conditional modelling is ill-posed (c.g., discrete y), or when the posterior distribution is available via generative model.
In Figure 8 we visualise an example of a 2D Crescent-shaped density. Both methods achieve the desired coverage, but the posterior model yields slightly smaller area (4.80 5.57), suggesting potential gains in efficiency in some cases. However, the posterior region appears less symmetrie, exhibiting stronger coverage along the upper edge of the density. It is still not clear if this asymmetrical nature can offer any advantages, or there is any added theoretical value in using the posterior-based conformity score.
True Density
JAPAN Likelihood
JAPAN Posterior
Coverage: Aree: 5.57
Coverage:
4,00
(a) Empirical
(b) Likelihood
(c) Posterior
Figure 8: Visualisation of JAPAN Likelihood vs JAPAN Posterior Crescent density dataset = 0.1. Both achieve the desired coverage, but the Posterior density has a smaller region for the given example, while the Likelihood density looks more symmetrie.
Algorithm 4 JAPAN with Posterior Modelling (p(x y) or p(u(z) v))
Require: Training set Duran {0 9:))% calibration set Deal {(z,, 0;) significance level
1)
Ensure: Prediction set (x( based on posterior log-likelihood conformity
1: Option A (full input posterior): Train a flow model p(z y) using {(xs, 34)
2: Option B (prediction posterior): Train a predictive model μ, then fit p(u(z) y) using {(((zz), vs)
3: for each calibration pair (zi)));) Deal do
4: For Option A: set aj = p(x, ys)
5: For Option B: compute ú; = p(r,) and set aj = log p(új v1)
6: end for
7: Compute threshold Te as the k-th smallest value {a 1m where k [c- 1/m) ml
8: At test time: For any input x, define:
(x = (u E y log p( w) 2 T.) or log log((p(x) v) 2 Te
depending on the variant used, and computed using Proposition 3
A.5 LATENT MODELLING
In the original JAPAN framework, the conformity score is computed with the full log-likelihood produced by a normalising flow as in Eq. (T):
log x)
log pz parppz (h(,, + (y)
20

---

Published as a conference paper at ICLR 2026
where the first term accounted for the density of the transformed variable = h(y, a) in the latent distribution (e.g., standard Gaussian), and the second term accounted for the local change in area/volume under the flow transformation.
Here, we intentionally discard the area/volumer correction term and define the conformity score using only the log-density of the transformed variable in the latent space (see Algorithm
log Pz (eppp2(((g (h
This corresponds to thresholding in the latent space without correcting for the local distortion introduced by the flow mapping. Intuitively, this results in spherical region centred at the origin in latent space, which maps to an arbitrarily shaped region in data space. However, the region may inadvertently include low-density regions of the target distribution: low-density regions that were simply squeezed into high-density latent locations may be mistakenly assigned high conformity.
This design has close connection to residual-based methods. The Gaussian base distribution naturally induces a conformity score that penalises distance from the centre, much like f2-based residual scoring. Indeed, if the flow is polume--reserving (i.e., log det Oh| = or if the area/volume change term is ignored, this formulation becomes equivalent to CONTRA (Fang et al. 2024). which bounds conformal balls in the latent space as their distance-based threshold corresponds working with the density of transformed variable in the latent space. Hence, CONTRA can be viewed as special case of JAPAN that uses only the transformed latent density, without accounting for how that density is warped back into data space. While this setup can produce valid regions, it often leads to large areas similar to CONTRA when modes are well-separated.
Formally, recall that h denote the normalising-flow mapping conditioned on J, so that
= huy,, x),
pz,
where pz the base distribution. We choose pz to be an isotropic Gaussian,
pz(2)= (2T))//2 exp (- 1-211 || 11=12)
so log-density is
1 log 1og = 111111 11 log(2n).
CONTRA's conformity score (latent-ball version) can be written
aCONTRA (y Φ) = -1| [[/ x2
Note that we change the original non-conformity score to conformity score for the comparison by changing the sign.
JAPAN's latent-density score is:
1 =1) - log((
d
aJ)pAN(y z) log (oggz(hh(y,a)) = - x)21
Then for any and x, we have
QJAPAN 2) = aCCOTRA( (x) constant.
Crucially, for any two candidate outputs V1, V2 (and fixed x))
aCONTRA (V1 x) < aCCOTRA(())
(JAPAN xx) 2 QJAPAN(32 z).
That is, the Euclidean norm in latent space and the latent Gaussian log-density induce exactly the same ranking over candidate y's and equivalently the same prediction regions/areas.
In Figure 9 we show comparisons of different variants on JAPAN with visualisation of areas computed at O.I for the Crescent-shaper density. Importantly, the coverage and the computed areas for CONTRA and JAPAN-Latent are identical, showing empirically that CONTRA is a special case within JAPAN
21

---

Published as a conference paper at ICLR 2026
True bersity
JAPAN Orip
JAPAN lau-Adaptive
JUAN Latert
CONTRA
Figure 9: Visualising prediction regions across JAPAN variants. From left to right: the true density, JAPAN Original, JAPAN Tau-Adaptive, JAPAN Latent, and CONTRA. All methods achieve comparable marginal coverage (near 0.89), but the shape and compactness of the prediction regions vary. JAPAN Orig applies fixed threshold across contexts, while Tau-Adaptive adjusts the threshold locally to account for varying volume distortion. JAPAN- and CONTRA both threshold in latent space, resulting in same areas and coverage showing JAPAN in latent equívalent to CONTRA. Ada-ivae yields improved coverage and competitive area, closely matching the data geometry.
Algorithm 5 JAPAN with Latent Modelling
Require: Training set Dirain = {(r,, yi) H=) calibration set Danl = {(1,, w)));" significance level E 1)
Ensure: Prediction set I" (z) using latent-space density as conformity score
1: Train a conditional normalising flow h1,, .r), where 2 = huy,, xr) and = ~ pz
2: for each calibration pair (æy, Uj) do
3: Compute latent variable: 2j h FT))
4: Compute conformity score: a₁ log pz (2z)
5: end for
6: Compute threshold as the k-th smallest value in {a where k = LE- (1 - 1/m) m]
7: At test time: For input x, define:
(e(x) (y E E logpz(b(u,z)) 2 Te)
A.6 ADAPTIVE (x)
While JAPAN produces context-adaptive regions and supports disjoint prediction sets, does not guarantee conditional coverage. This limitation arises because the conformal threshold т, is computed globally, independent of the conditioning variable æ. As result, although the shape of the region adapts with æ, its size may be suboptimal or overly conservative in certain regions of the input space.
To address this, we propose tau-adaptive extension of JAPAN, in which the threshold itself is a function of x, i.c., T (z). Our approach builds on the fact that the log-likelihood in a normalising flow decomposes into two terms, from Eq. (T):
log ()) log (h(,, + + xx)
latent density
area change
We use this decomposition to split the conformity score into two components. The unconditional component is based on the latent density and is calibrated globally using the standard conformal approach. This yields a latent-space threshold, which corresponds to the same threshold used in our latent modelling extension Algorithm 5 CONTRA Fang al. (2024). Formally, we compute the global latent-space threshold F, from the calibration scores log pz (h(y), (T))), and additionally store the corresponding latent vector z, for the e-quantile (i.c., the k-th smallest score). At test time, we use this fixed latent value to construct context-specins threshold in the data space:
T (x) F. + (h (0h--1(zp)) (z,, x))
This formulation ensures that the time threshold accounts for local flow transformations by evaluating how the celibration-latenn 2, maps back to data space under the test-time context æ. Further details on the algorithm are in Algorithm
22

---

Published as a conference paper at ICLR 2026
Intuitively, this adjustment of the threshold using the area-change component reflects how the flow expands or contracts volume around æ. When the Jacobian suggests strong expansion (i.e., high uncertainty), the threshold is lowered, allowing larger region; conversely, in contracting regions (i.e., high confidence), the threshold is raised. It can also be viewed as aligning the latent-space decision boundary with the true conditional-density contour in the data space, and thus restoring degree of conditional coverage. We may even get nearly optimal prediction regions under certain assumptions (e.g., when ply 2) transforms linearly in x and thus maps the true thresholds the same latent vector, and the flow estimates conditional density sufficiently well across z). It improves upon other context- addattiv methods, which only modulate shape by shifting and not scaling, and implicitly assume area-preserving transformations. An alternative to using the latent vector obtained through global calibration can be find the local threshold using the neighbours the context, which also brings JAPAN closer to conditional coverage or more specifically, local coverage. Algorithm 7 shows the full algorithm for this proposal.
we show the prediction regions obtained by the tau-adaptive extension on the Crescent dataset. While the resulting region is slightly less compact than those produced by CONTRA or other JAPAN extensions (with area 6.99 compared to 6.46 and 6.46, 5.99, respectively). it achieves slightly improved coverage (0.91), reflecting the influence of threshold modulation through local volume expansion. This suggests that adapting the threshold based on flow geometry can provide useful mechanism for refining predictive regions, even if it may not always lead to the smallest volume. Unconditional coverages can easily form tighter regions by excluding high-variance predictions, but this comes at the cost of non-uniform calibration conditioned on context. Although the tau tau-ataptive variant introduces modest increase in area, may offer practical pathway toward better conditional coverage.
In Figure 9
Illustration on a Toy Example: To further highlight the role of the tau-adaptive extension, we design simple heteroscedastic regression problem. Let
x~ 40, 1],
y az xe,
(0,11
so that the conditional distribution is
y | x ~ N(ax, x²).
This setting induces a conditional variance that increases with æ, making it a natural testbed for evalu- ating conditional coverage. Since the latent mapping is fixed across 2, the tau-adaptive mechanism can be applied analytically, and the adaptive threshold Te Tz() adjusts to the local variance.
As shown in Table 7 JAPAN with tau-adaptive thresholding achieves the desired conditional coverage (0.90) while also adapting the prediction region width appropriately: smaller when rr = 0.5 and larger when 1, consistent with the variance structure of the model.
Table 7: The conditional coverage for toy example, where our assumptions (latent mapping and hence 7 is the same for different values of I)) are satisfied.
Method
Width
Coverage
JAPAN (r = 1)
3.99 + 0.01
0.90
JAPAN (x = 0.5) 1.64
0.90
Since our primary goal is to highlight the flexibility of flow- conformity scoring, defer deeper investigation into conditional coverage guarantees to future work.
A.7 ABLATION: EFFECT OF UNDERFITTING IN FLOW MODELS
We begin by recalling the role of conformity scores in our framework. In traditional inductive conformal prediction, the predictor μ and the scoring rule A are defined separately, e.g. using aesidual-based scores such as a; = - 1y/ fuzz)). In contrast, JAPAN has multiple settings, for e.g the first setting as in Section 3.3 of JAPAN unifies these components: the normalising flow is not only used to obtain a predictive model, but also directly provides the conformity score via its estimated ligeliiel,o,d, log p(y: I..). This removes the need for separate residual-based metric
23

---

Published as a conference paper at ICLR 2026
Algorithm 6 T-Adaptiive JAPAN
Require: Training set Drrain {Ti Ui) |- calibration set Deal {0z Mj) w;));= significance level
CE 1)
Ensure: Prediction set (2) with geometry-aware adaptive thresholding
1: Train conditional normalising flow: h(y, with base density pz
2: for each calibration pair (zz, Vj) do
3: Compute latent representation 2j Iz)
4: Set conformity score: aj = log (zz)
5: end for
6: Compute global latent threshold F. as the k-th smallest value {az}ā= where t: = Le
1/m) m. and store the corresponding latent vector Ze.
7: At test time: For input z, define the data-space threshold:
Te = fe + + (h (Ze, x),
8: Predictive region:
r.(x) ty E y logpz(h(y), + 2z), 2. (z)2))
Algorithm 7 Te (dda)--dapiiv JAPAN with Neighbour-based local adaptivity
Require: Training set Dirain = {z1 y:) 1 P calibration set Deai {(z,, U3) 1 th significance level 6 1)
Ensure: Prediction set (x) with locally adaptive thresholding
1: Train conditional normalising flow: h(y,x), with base density pz
2: for each calibration pair (z. U)) Deal do
3: Compute latent variable 2j ,,))
4: Store: latent conformity Og (zg correction term 4; D I.y))
5: end for
6: Compute global latent threshold F, the k-th smallest in {a,} where k = [e- (cc-11114 (111mmmmmm) -m]
7: At test time:
For input
Compute ((x) mean of o over the K- K-nearest calibration inputs xy
Set threshold T (2:) = + ((x)
Define:
Te(x) = fy e y log (h + 2(y, 2 [z())}
24

---

Published as a conference paper at ICLR 2026
and ties prediction and conformity into single probabilistic model, whereas the second setting, as in Algorithm 3 used a separate base predictor.
A natural question is how the performance of JAPAN degrades when the underlying flow is underfitted, particularly in high-dimensional settings. While coverage guarantees remain valid regardless of the flow quality, the efficiency of the prediction regions, i.e., their area, may degrade. The extent of this degradation depends on how poorly the flow model approximates the true distribution. Typically, normalising flows trained via maximum likelihood tend assign higher density frequently observed samples, which tends to preserve ranking structures reasonably well even under underfitting.
We first revisit the setting where both predictor μ and conformity function A are defined jointly by the normalising flow. The central question here is whether the prediction regions degrade if the normalising flow is not well-trained. The answer parallels what happens in standard conformal prediction when the predictor is underfitted: in both cases, a poorly trained model results in unreliable conformity scores and hence inefficient or wide prediction sets. Poor flow training is a limitation, but it only widens the prediction sets to account for that uncertainty, which is inherent to all conformal methods.
In fact, all of our main baselines also use normalising flows (for the mean predictor), but rely on residuals or geometric heuristics for scoring. If the flow model is inaccurate in those cases, the reliability of the prediction regions also deteriorates. In this ablation, where we reduce the number of layers to half, we see in Table that JAPAN still achieved more efficient sets compared to baselines, even though the flow was underfitting. We posit this is because the model still learns reasonably well to put high density on the frequently occurring outcomes.
Table 8: Underfitted modelling: Coverage and prediction region area on the Energy and RF4D datasets. Lower area indicates better efficiency at the same coverage.
Method
Coverage (Energy)
Area (Energy)
Coverage (RF4D)
Area (RF4D)
JAPAN
0.89 = 0.01
25.61 + 2.44
0.90 = 0.01
407.68 + 44.19
CONTRA
0.89 + 0.01
28.13 + 6.18
0.90 = 0.01
690.99 + 108.35
CopulaCPTS
0.88 + 0.01
40.63 + 12.33
0.91 = 0.01
2991.77 + 5123.00
Importantly, in this setting, the conformity function A and the predictive component ls were jointly estimated by the normalising flow. This coupling makes it difficult to isolate underfitting effects on A or j individually. In contrast, most baseline methods relied on fixed predictors and separately defined conformity scores, simplifying such analysis.
To further investigate this, we draw attention to the second setting, i.e. the experimental setting in Table which decouples the predictive mean from the conformity score. The results showed that applying the normalising flow improves uncertainty estimation of the base model. We follow the same setting here, and purposely underfit the normalising flow. Note that the normalising flows are only used here in JAPAN, and thus the prediction of other models (baselines) remains constant that JAPAN still yielded much better results, even though it was
from Table 6 In Table 9 underfitting with only half the layers used for training.
Table 9: Underfitted modelling after decoupling scoring and predictor model: Coverage and prediction region area the Energy and RF4D datasets. We chose two non-flow- flow-base baselines that worked the best in this setting. Lower area indicates better efficiency at the same coverage.
Method
Coverage (Energy)
Area (Energy)
Coverage (RF4D)
Area (RF4D) x 1000
JAPAN
0.90 + 0.01
76.19 + 4.92
0.89
0.01
9.91 + 1.12
Dist-Split
0.89 + 0.03
151.26 + 56.65
0.83
0.01
20.69 + 3.14 77.15 + 8.72
RCP
0.87 + 0.02
124.46 + 29.45
0.90 + 0.01
In summary, while underfitting reduces efficiency across all conformal methods, JAPAN degrades less. Because the log-likes likeelooooood estimated by the flow simultancously governs both prediction and conformity, JAPAN remains competitive even when the flow is not fully expressive.
25

---

Published as a conference paper at ICLR 2026
Table 10: Samples required for comparable accuracy on the test set with sampling schemes sampling. d = 4 and d 8.
Method
Samples (d = 4)
Time (minutes)
Samples (d = 8)
Time (minutes)
Proposition 3 (ours)
1,000
0.4
3,000
2.5
Naive MC
100,000
1,000,000
32
B PROOFS
Proposition 1
Optimality Under Rank Preservation]
Let ply | z) denote the true conditional density, and let fely x) denote model-estimated conditional density. Consider the prediction region defined at coverage level 1 c:
(x) = {v p(y a:) 2 r"). where PP(y E r (1)) = -,6
and the conformal prediction region defined via the model as:
[x(x))= (v fo x) 2 T}, where Pyy E r rx(x)) 1
Assume:
There exists strictly monotone increasing function g RR + R+ such that fo z) g(p(y z)), Vy.
The conditional densities p(y x) are homogencous across H, i.c. there exists a base density & and a translative transformation T3 such that p x) = (TT( for all æ.
Then the model-defined conformal region achieves the same area the true optimal region: (z)) (r))((" for all
The homogenity condition assume that the shape of the conditional densities p(y x)) remains the same across different values of up translation based transformations. More generally, this means that the conditional distributions share the same spread or covariance structure but may differ in location. This assumption is standard in settings where the conditional variability is constant but the conditional mean changes. Examples include Gaussian location families with fixed covariance, conditional distributions that differ only by shift parameter.
Remark:
Proof. We know that the prediction region at coverage level 1 - under p(y x) is:
(x) (y p(y æ)) T"}, with [ p(y x) dy 1
(x)
By assumption, foly w) = g(p(y x)) for strictly increasing function g. Hence, its inverse g- is also strictly increasing, and:
ply x)) g (fo(() x))
Then the model-based conformal region is:
T. (z) = {v fo 2 (v plu x)) 2
Since both regions achieve the same coverage level under the same true distribution p(y 2), they must correspond to the same threshold:
(T)
T=g(T")
implying:
[z)
Finally, under the homogeneity assumption, the structure of the density and its level sets invari- ant (up to transformation) across different values of z. Therefore, area computations under the transformation Te preserve the equality:
Arca(T = (z))
26

---

Published as a conference paper at ICLR 2026
Structural Equivalence Under Monotonic Transformation: Even without assuming homogene- ity of the conditional densities pyy 2). if the model density fo(y z) is a strictly monotonic transformation of the true density (i.e., fo x)) = g(p(y æ)) for some fixed, monotone g), then the prediction regions defined by thresholding either density are identical for each æ, and hence their expected areas match:
T (x (2)
[Area((T (x))) 1R, (r)))()
This does not imply that the model-based region is optimal, but rather that it is structurally equivalent the one derived from the true density under the same thresholding scheme.
Implication for JAPAN: Since the logarithm is strictly monotonic transformation, this proposition also justifies using log-likelihood (i.e., log p(y x)) as conformity score in JAPAN, without loss of structural validity in the regions.
Proposition 2 Approximate Optimality under Approximate Rank Preservation] Let p(y x)) be the true conditional density and fe x) an estimated conditional density. Suppose there exists a strictly monotone increasing function IR+ + R+ and a positive constant ó such that for all y and æ:
ffo x) - g x)) S
Assume further that the conditional densities p(y a)) homogeneous across x, i.e., there exists base density & and a translative transformation T. such that:
ply x) = (Ta(v)) for all x.
Define the optimal and estimated prediction regions respectively as: r (x) Ey: p(y æ) 2"),, re(x) = Ey folu 2
both satisfying coverage:
Ply ((yE1)())) Ply (x))
Then the area difference is bounded and approaches zero as d + 0: Area(T - Area(T) < C(6), with C(²) - 0 as 5 -
Proof. By definition of coverage, we have:
J ply z) dy = =1 1 and [ ply x)) dy =
r.(m))
Since foly 2) approximates g(p(y) x)) uniformly within 6, we have:
[e(x) {v foly x) 2 r} g((() x)) 2 T).
Thus, the thresholds 7" and (T) must be close. Formally, continuity and strict monotonicity of and g imply:
T - (T)
with e0) + 0 as 8 + O.
Consider the symmetric difference of the two regions:
(*))A (x) {y ply x)) [min(7" ,9" (T)), max(7", ,9 (7)))]-
This set corresponds to a thin boundary region around level sets defined by ply æ)), with width at most (8))
The area difference is bounded by the area of this boundary set: (Area(T(())) - Area(T" < (rra(T(())AT,aa)))
By the homogeneity assumption, the shape and structure of level sets across different æ are invariant (up to transformation), so the area difference remains controlled uniformly. Hence:
Ax) < C(6), with C(ő) - as + 0.
27

---

Published as a conference paper at ICLR 2026
Proposition 3 [Area (or Volume) Computation] Let pz(2) be samples from base distribution, and let y = (2, x) denote the inverse flow transformation conditioned on z, and exp(@(2, 22 z)) denote the area/volume correction under the inverse transformation. Then the area (or volume) of the conformal prediction region (x) can be approximated as:
Arca(Te(s)) = N 1=1
1
Σ
N
(p(h x)) x)) 2 TT]
exp(o(2,x))
Pz(zi)
where 2, pz for i ....., N, and T2 is the conformal threshold.
Proof. Let F()) == 1 {p(h- æ) x) 2 Te eexxexx(((22xx) Then the area of the conformal region can be written as an integral over the data space:
{ JRD.
Area(T.(z))
1 LE(y x) 2 Te) dy.
Changing variables via the inverse flow y = (2, we obtain:
Area(T.(2)) RRD
dz.
Expressing this as an expectation under 2 ~ pz(2), we get:
F(2)
1 (iph (2,x)) (.r)) 2r] exp(ó(2,x))
Area(T.(s)) Egopp
Esopz
[2((z)
Pz(z)
This expectation can be estimated with N Monte Carlo samples {zp p2) yielding:
x)) 2 2 exp()(z.x))
Area(T.
Σ
1
Pz(2;)
This estimate accounts for area/volume transformation under h-1 and evaluates only those points whose mapped outputs lie within the conformal region. When using discrete Normalising Flows, we can write down exp(o(2, x)) as |det Jn- 1 (Z, x)). where Jn 1 is the Jacobian under the inverse transformation. Then, the computed area takes the form:
Area(Te
1
Σ 1=1
N
{p(h (0(h'(22,00)) (x)2 27.
Idet (Zi,
=
N
P2(=;)
B.1 ALTERNATIVE VIEW OF VOL.UM//ARAA ESTIMATOR IMPORTANCE SAMPLING
Whilst sampling from the latent space is effectively equivalent to sampling from the data distribution, and contributes the efficiency of the estimation process. We can also view the volume/area estimator as an importance sampling estimator, which helps explain both its correctness and practical efficiency.
Using the Shange-of-variabee formula, the volume of the prediction region
Tx(x) = tu log foly (x) 2 }}
can be written as
Area(T. -fr
del Jn (n,+(2 dz.
[(T)(x)))
Here, h the transformation from the latent space to the data space, and Jn-1 is the Jacobian of the transformation.
By introducing the latent Gaussian Pz(z) and multiplying and dividing by pz this becomes the expectation
28

---

Published as a conference paper at ICLR 2026
det Jn (2, 1 10h e e.(co))
Area(T( = Ez-pz
pz(z)
which is exactly the form of an importance sampling estimator with proposal distribution pz. The weights
w(2, 2) det JA-1 (2, x)
(z(z)
the corresponding importance weights, and the indicator selects points in the prediction region. Our Monte Carlo estimator is therefore an Importance Sampling estimator:
M Σ
1
Area(T.(x)) M w022 x) 1 (2,, E (x)),),
24 ~ Pz.
This Importance Sampling viewpoint also explains the efficiency of the estimator. Because normalis- ing flows Jypically, map the data distribution to Gaussian latent pz. and the data variable y can be expressed as function of 2, the Importance Sampling viewpoint of using the Gaussian as the proposal distribution. This avoids the inefficiency of naive grid-sampling or uniform sampling in the output space. In other words, the proposal distribution pz is well-aligned with the transformation-induuee distribution of the data via flows, which also reduces the variance of the estimator, as the weights are naturally obtained from the volume change component that can be analytically calculated via flows.
Proposition 4 (Bounded Misranking Leads to Bounded Excess Area). Lel ply x) denote the true conditional density, and let fo x)) denote the estimated conditional density from model. Define the conformal prediction region and the optimal region at coverage level E as:
(x(x)) by foly xx) r},
[7)x) ty ply x)) 2T}
with thresholds T and T* chosen so that:
P(y E (x)))) = Ply E (x)) =1
We use Arcal-) to denote the Lebesgue measure of a region Du. Define the misranking probability mass:
u(x)= 11 plyn p(y2 x:) dyn dyz, where Ma (11183)) p(un 1 > p(u2 ()), folur (x) < Soly2 x)).
Then there exists a function C(-) with C(u) - 0 as 0 such that:
(Area(T((z)) - Area(["(z)) s C(u(z)).
Proof. Let consider the symmetric difference between the two regions:
(r(x) (a: = (T.(x) (TT()))))) \ Tx()))
Any point y in this symmetric difference satisfies:
(p(y) x) 7" (fo(( - <.
which means that y lies above threshold and below the other precisely due to ranking discrepancy between fo and p near the threshold boundary.
Now consider the set of all such y for which this misalignment occurs. These generated by misranked pairs- - i.e., pairs (y1, y2) such that p(y1 x)) > p(yz x)) but fo(y1 x) < fo x).
The total measure of such misranked pairs is quantified by p(x:).
When p(x) is small, the threshold r will be close to g(T*) (assuming fo approximates a monotonic transformation of p), and the region boundaries will be close. This implies that the symmetric
29

---

Published as a conference paper at ICLR 2026
difference between T. (x) and (x) is thin geometrically, it lies in strip of points near the threshold.
Therefore, the area of the symmetric difference is bounded by function C(u(z)), which goes to zero as ju(x)= + 0:
Area(T((2)AT)(o)) C(u(a)).
Since:
Area (r))(I {rre([[())][[[[]))
conclude the result.
C EXPERIMENTAL AND ARCHITECTURAL DETAILS
C.I TARFLOW FOR TIME SERIES MODELLING
We adapt the TARFLOW architecture (Zhai et al. 2024 originally proposed for image patch modelling, to the time series setting by treating the forecast window as a sequence of tokens and conditioning the flow on the past context. Each time series forecast target y E R is modelled autoregressivell using a stack of Transformer bbasdd flow blocks. These blocks, as in TARFLOW, contain causal Transformer followed by diagonal affine transformations applied permuted inputs. Unlike the image setting, where the input is a sequence of visual patches, we treat the temporal axis the natural ordering and apply causal masking accordingly.
To incorporate conditioning on the historical context æ E R2'xd, we introduce an additive conditioning path. Specifically, we project the input æ through a context encoder and map it sequence lgd of the same length and shape as the target y. This projected context is then added elementwise the latent representations inside the Transformer flow block:
1+ = f2(22) + context encoder(I:),
where f' is the Transformer autoregressive block, / is the layer index, and context encoder(x) denotes the conditioning signal projected to match the forecast window.
Briefly, each flow block f' follows a block autoregressive architecture. The core operation consists of three components:
1. Causal Transformer: standard Transformer with causal masking (along the time axis) processes 2' to produce context-based representations. The attention mechanism ensures each timestep z! can only attend to past and current positions.
2. Diagonal Affine Transform: The output of the Transformer is split into scale and shift vectors, (((2))) and 3 (2(gg) which parameterize an affine transformation:
++ (1
exp(->'(E<2))
This results in triangular Jacobian with efficient deterenermnnant computation. 3. Permutation: After each flow block, fixed permutation reversing the order of the variables is applied the sequence dimensions to improve expressiveness across layers.
This additive conditioning mechanism preserves causality, introduces minimal architectural overhead, and allows the model to adaptively shift the transformation of y based on the input history. This structure preserves TARFLOW's benefits: exact likelihoods via tractable Jacobians, causal modelling of temporal dependencies, and a parallelisable attention-based architecture.
C.2 TABULAR EXPERIMENT DETAILS
For the tabular benchmarks, we closely follow the setup used in CONTRA (Fang et al., 2024). Specifically, we employ a discrete normalising flow architecture (Dinh et al. 2017), which enables efficient sampling and density evaluation. Each coupling layer uses two separate neural networks to parameterise the scale and shift transformations. These networks consist of two hidden layers with 512 hidden units and ReLU activations. We use coupling layers and a batch size of 512 across all datasets. Training performed for 200 epochs using the Adam optimiser with learning rate of
30

---

Published as a conference paper at ICLR 2026
10m3, along with an exponential learning rate scheduler that decays by a factor of 0.999 after each epoch.
Each dataset is split into 60% training, 20% calibration, and 20% test sets. For evaluating region areas in normalising flow-based methods (JAPAN, CONTRA, PCP), we generate 3,000 Monte Carlo samples per test point. We report coverage and area statisties based on 25 random seeds for the main experiments, and 5 random seeds for the various JAPAN extensions.
For baselines that require additional modelling, involves training a model on residuals, we use 60% of the calibration set for fitting a Lasso-based regressor (from scik ticlkit-leate rn) and the remaining 40% for calibration. Hyperparameters selected using cross- validation. For CopulaCPTS, we similarly use 60% of the calibration set to train the copula model, and the remaining 40% for conformal calibration. For all JAPAN extensions evaluated on the Energy and RF4D datasets, we use simple linear regressor from scik k i -t-1a lea rn as the fixed base model p(z) generate point predictions used in the conditional and posterior variants.
the following procedures: for JANET, which
C.3 TIME SERIES EXPERIMENT DETAILS
For the time series benchmarks, we use autoregressive flow architecture adapted from TARFlow (see Section C.T.. Each model is trained for 200 epochs using a batch size of 128 and a learning rate of 10-3, with an exponential learning rate scheduler (decay rate of 0.999) applied at each epoch. We first validated the hyperparameters on a simulated dataset and then reused them across all real- world experiments.
The attention-baan based flow model includes 10 autoregressive flow layers. Within each attention block, use 2 attention heads, with 32 52-dimensional head size and 64-dimensional time-step embeddings (via expansion). Each residual block the flow is expanded by a factor of in hidden size. The input conditioning signal is processed through context encoder comprising 2 stacked LSTM layers, each with a hidden dimension of 64.
For methods with secondary models such JANET, allocate 60% of the calibration set train the residual regressor and reserve 40% for calibration. A Lasso-based regression model from scikit-leazn used, with hyperparameter tuning performed via voos-vvlliation. Similarly, for CopulaCPTS, we use 60% of the calibration set to fit the copula model and the remaining 40% for calibration. We evaluate all methods across 25 random seeds for the main experiments and report mean coverage and prediction area across runs.
C.4 TOY DATASETS EXPERIMENT DETAILS
We evaluate method across five synthetic 2D distributions: Spiral, Checkerboard, Moons, Circles, and Crescent. Each dataset consists of 10,000 samples, split into 6,000 for training, 2,000 for calibration, and 2,000 for testing. To visualise the empirical density, we also generate Ti samples for each dataset.
1000
For the density model, with residual blocks and hidden layer dimensionality of 32. All models are trained using the Adam optimiser with batch size of 512 for 200 epochs. We use learning rate of 10 -3 with exponential learning rate scheduler (decay rate of 0.999) applied at each epoch.
normalising flow architecture composed of 4 coupling layers, each
Compute details: We conducted all experiments using NVIDIA A100 GPUs. Including failed runs, hyperparameter tuning, and iterative experimental development, the total compute time across all models is estimated to be approximately 30 GPU-days.
D DATASET DETAILS
D.1 TABULAR DATASETS
We evaluate our models on several standard tabular datasets frequently used in probabilistic regression and conformal prediction literature:
31

---

Published as a conference paper at ICLR 2026
Dataset
dim(x)
dim(y)
Type
n1 (Train)
n2 (Calibration)
ng (Test)
Energy
8
2
Tabular
460
154
154
2D RF
20
2
Tabular
5403
1801
1801
4D RF
20
Tabular
5403
1801
1801
SCM20D
61
2
Tabular
3000
1000
1000
COVID
100 X 1
10
Time Series
200
100
80
Particle-1
55 x 2
4 x 2
Time Series
2000
500
500
Drone
57 2
Time Series
600
200
200
Pedestrian
16 X 2
4 X 2
Time Series
891
200
200
Table 11: Summary of dataset structures used in our experiments. Datasets are grouped into tabular and time series types. Dimensions and splits are given for input (2)), output (y), and each of the train, calibration, and test sets. For time series datasets, b reflects the time steps the dimensionality of the time series.
Energy (Tsanas & Xifara 2012): This dataset consists of measurements from various buildings to predict heating cooling loads. includes eight features such as relative compactness, surface
wall area, and roof area. The task is to predict the energy efficiency of buildings based on their design parameters.
RF2D and RF4D benchmarks derived from the Mississippi River network. The dataset consists over year of hourly river flow observations from eight measurement sites. In the 2D RF setup, the model is trained to predict river flows 48 hours ahead at 2 downstream sites using past flows. In 4D RF the forecasting extended to 4 downstream sites.
iictomitros-xooufie et al. (016): These are multivariate regression
et al. 2016): This dataset comes from the 2010 Trading Agent Competition in Supply Chain Management, It contains 5,000 records, with the goal of predicting the average price of a supply chain product 20 days into the future based on set of observed features.
SCM20D Espyromitros-Xionee
D.2 TIME SERIES DATASETS
We evaluate our method diverse set of synthetic and real-world time series forecasting datasets:
COVID (Stankeviciute et al. 2021): ): We follow the experimental setup of English et al. (2024), predicting daily case counts in UK regions. The input consists of 100 days of historical data, and the model forecasts the next 10 days. The dataset includes 200 training series, 100 for calibration, and 80 for testing.
Particle-1 (Sun & Yu 2024): These datasets simulate interacting particle systems using the setup of (Kipf et al. 2018). For both, the model predicts future positions 11+1 where
1 -= 55 and 7 = 4. We added Gaussian noise with standard deviation 0.01 to the dynamics. Each dataset includes 5,000 samples, with a 60/20/20 train/calibration/est split. Drone (Sakai et al. 2018): Simulated drone trajectories are generated using the Python- Robotics simulator. At each time step, the input is a 10-dimensional vector (h = 57) and the prediction target is 3-dimensional vector (y: E )
Pedestrian (van den Berg et al., 2008): Forecasts human trajectories based on data from the ORCA simulator. Each sequence consists of 2D position measurements over T 20 time steps for 1,291 pedestrians. This dataset captures multimodal future behaviour and dense agent interactions.
E CALIBRATION PLOTS AND AREAS UNDER DIFFERENT THRESHOLDS
To assess the stability of different methods, we evaluate their calibration performance across a range of significance levels [0.05, 0.95]. Specifically, we plot calibration curves measure empirical coverage and visualise the corresponding prediction areas at each level, providing insight into both validity and efficiency under varying uncertainty thresholds.
32

---

Published as a conference paper at ICLR 2026
Dataset)
Prodiction Areas Coverage Dataset)
Empriiall Covergge
58
Asstroge
Nominal
Nominal
(a) Moons dataset. Coverage and Area for various significance
Cal lbbation (Circle Dataset)
Prediction Nominal Coverage (Circle Dataset)
2.0
30
Emplical Coverage
(ALE
Awcrage Arpa
20
64-
50
Nominal Coverage
Nominal
(b) Circle datasct. Coverage and Area for various significance levels
Figure 10: Calibration and efficiency curves across toy datasets. The left panels shows the calibration curve. plotting empirical versus nominal coverage for significance levels e [0.05, 0.95]. A perfectly calibrated method lies along the dashed diagonal. The bottom panel shows the corre- sponding average prediction region areas. JAPAN (shown with bold line) achieves consistently valid coverage across all levels while maintaining lower or comparable region volumes relative to other methods, highlighting its balance between calibration and efficiency.
Figures 10 to 15 show these curves across multiple datasets. In each plot, the left panel presents the calibration curve, where perfect calibration corresponds to the dashed diagonal. The bottom panel reports the average prediction area across the test set at each significance level. JAPAN is highlighted with bold grey line for easy comparison.
JAPAN consistently maintains empirical coverage close the desired level across all thresholds. demonstrating strong calibration stability. In contrast, several baselines-includingg and Copula-based methods exhibit over-coverage at higher significance levels, indicating instability and excessive conservativeness. Notably, JAPAN also achieves the smallest or among the smallest prediction areas across almost significance levels.
including CQR, Dist-Split,
USE OF LARGE LANGUAGE MODELS
We used large language models (LLMs) solely for writing assistance, limited to polishing grammar and improving clarity of presentation. No LLMs were used for research ideation, theoretical results, experiments, or analysis.
33

---

Published as a conference paper at ICLR 2026
Calibration Curve (Checkerboard Dataset)
Prediction Areas vs. Nominal Coverage /Checkerboard Dataset)
NCOR
13
50-
034
Eme Coverage
56
20
20-
52
0.0
0.6
0.8
1.0
0.2
05 Nominal Coverage 11 d
0.8
Nominal Coverage (1
(a) Checkerboard dataset. Coverage and Area for various significance levels
Callbration Curve (Spiral Dataset)
Prediction Areas vs. Nominal Coverage (Spiral Dataset)
1.0
35
DON
FCP
NLz
0.8
Oge
Empiniall Coverage
25
CopulacTTS
Viri a0roub
06
20
04
35
10
0.2
0.0
6.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
Nominal Coverage (1-8)
Nominal Coverage (1
(b) Spiral dataset. Coverage and various significance levels
Figure 11: Calibration and efficiency curves across toy datasets. The left panels shows the calibration curve, plotting empirical versus nominal coverage for significance levels e E [0.05, 0.95]. A perfectly calibrated method lies along the dashed diagonal. The bottom panel shows the corre- sponding average prediction region areas. JAPAN (shown with bold line) achieves consistently valid coverage across all levels while maintaining lower comparable region volumes relative to other methods, highlighting its balance hetween calibration and efficiency.
34

---

Published as a conference paper at ICLR 2026
Calibration (5CM20D Dataset)
Prediction Areas Nominal Coverage (5CM20D Dataset)
1.-
250020
200020
Emerical Coverage
CPEA
(AL
350030
Aeverage
6.-
200000
30020
50-
Nominal Coverage
Nominal Coverage
(a) SCM20D dataset. Coverage and Area for various significance levels
Callbration Curve (Energy Dataset)
Prediction Areas Nominal Coverage (Energy Dataset)
120
HOOK
Creeu
Emplrcal Coverage
JNN
Average
64
0.2
3.0
0.4
0.8
Nominal Coverage (1
Nominal Coverage 12 [1
(b) Energy dataset. Coverage and Area for various significance levels
Figure 12: Calibration and efficiency curves across tabular datasets. The left panels shows the calibration curve, plotting empirical versus nominal coverage for significance levels [0.05, 0.95]. perfectly calibrated method lies along the dashed diagonal. The bottom panel shows the corre- sponding average prediction region areas. JAPAN (shown with bold line) achieves consistently valid coverage across all levels while maintaining lower or comparable region volumes relative to other methods, highlighting its balance between calibration and efficiency.
35

---

Published as a conference paper at ICLR 2026
Prediction Coverage
1.-
ConTA
R27
Cap SCCPTS
Empitial Coverage
JAETT
24-
Mecrage Ar3
54-
10-
5.-
Nominal Coverage
Nominal
(a) RF2D dataset. Coverage and Area various significance levels
Calibration Curve (RF4D Dataset)
Prediction Areas Nominal (RF4D Dataset)
12000
CONTM
Mas
MCOR
Emelical Coverage
Average
64-
0.
0.2
0.4
0.8
1.0
02
0.4
0.6
Nominal Coverage (1-
Nominal Coverage [1-
(b) RF4D dataset. Coverage and Area for various significance levels
Figure 13: Calibration and efficiency curves across tabular datasets. The left panels shows the calibration curve, plotting empirical versus nominal coverage for significance levels [0.05, 0.95]. perfectly calibrated method lies along the dashed diagonal. The bottom panel shows the corre- sponding average prediction region areas. JAPAN (shown with bold line) achieves consistently valid coverage across all levels while maintaining lower or comparable region volumes relative to other methods, highlighting its balance between calibration and efficiency.
36

---

Published as a conference paper at ICLR 2026
Cal Calibration Dataseti
Prediction Nominal Coverage Dataset)
1.-
Emerical Coverage
Avarage Aesia (41000)
6.-
50-
Nominal Coverage 11
Nominal Enverage (1-
(a) Drone dataset. Coverage and Area for various significance levels
Calibbattion Curve (Pedestrian Dataset)
Prediction Areas Nominal Coverage (Pedestrian Dataseti
NOOK
Emplrcal Coverage
Average Araa
20
JNN
25
64-
6.
2.4
62
04
0.8
Nominal Coverage 11
Nominal Coverage 12:41
(b) Pedestrian dataset. Coverage and Area for various significance levels
Figure 14: Calibration and efficiency curves across time series datasets. The left panels shows the calibration curve, plotting empirical versus nominal coverage for significance levels [0.05, 0.95]. perfectly calibrated method lies along the dashed diagonal. The bottom panel shows the corre- sponding average prediction region areas. JAPAN (shown with bold line) achieves consistently valid coverage across all levels while maintaining lower or comparable region volumes relative to other methods, highlighting its balance between calibration and efficiency.
37

---

Published as a conference paper at ICLR 2026
Callbration
Prediction Nominal Coverage (COVD Dataseti
1.-
1752
1502
0
Emerical Coverage
1252
(ALI
Avergge Area
3002
6.-
50-
Nominal Coveraag
Nominal Coverage
(a) Covid dataset. Coverage and Area for various significance levels
Calibration Curve (Particle-1 Dataset)
Prediction Areas vs. Nominal Coverage (Particle-1 Dataset)
HOOK
Emplrcal Coverage
Avoraige Area (x1000)
onA
JAN
64-
63
0.4
0.0
Nominal Coverage 11
Nominal Coverage (1-
(b) Particle-1 dataset. Coverage and Area for various significance levels
Figure 15: Calibration and efficiency curves across time series datasets. The left panels shows the calibration curve, plotting empirical versus nominal coverage for significance levels [0.05, 0.95]. perfectly calibrated method lies along the dashed diagonal. The bottom panel shows the corre- sponding average prediction region areas. JAPAN (shown with bold line) achieves consistently valid coverage across all levels while maintaining lower or comparable region volumes relative to other methods, highlighting its balance between calibration and efficiency.
38


Reviews:
[{"review_id": "09XFaufsg8", "reviewer": "Reviewer_MeP1", "summary": "This paper proposes JAPAN, a novel conformal prediction framework that defines flexible prediction regions using normalizing flows. Unlike CONTRA that relies on latent spaces to form conformity scores, JAPAN directly uses the loglikelihood as scores.\n\nThis enables the prediction sets to capture arbitrary, disjoint, and non-convex, non-connected shapes, while satisfying finitesample,\ndistribution-free marginal coverage guarantees.\n\nThe approach is evaluated and performed well on synthetic datasets, multivariate regression. The method has also been extended to handle interesting time series tasks.", "questions": "My questions are embedded in my comments for Strengths and Weaknesses above. Thanks.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Presentaiton-wise, the paper is well-structured and highly informative.\n\nEssentially, JAPAN can provide flexible prediction sets because it is based on estimating conditional densities, and normalizing flow is a good tool to achieve this.\n\nThe application to Time series is very interesting. (I am not an expert of CP used in this area of application. I do see there are multiple papers cited. Could the authors briefly talk about how the exchangeability-backed CP theory works in times series applications? Any tips on checking the validity of marginal coverage in theory and in practice? Given a new timeseries, how do I know if JAPAN produced CP are trustworthy?)\n\nThe paper provides some theory of coverage guarantees, rank preservation, and a volume estimator that supports efficient sampling.\n\nThe method is evaluated on over mulitple datasets, including time series data, and compared with many alternative methods.\nJAPAN consistently achieves lower prediction volume while maintaining coverage level.", "weaknesses": "The authors talked bout how to estimate the area of the proposed CP using Monte Carlo, a sensible solution. In contrast, there was a lack of mention of how to identify points y in the proposed CP, \\Gamma_\\epsilon(x). That is, even if one has the analytical form of the (log) probability density function, looking for all the points with density higher than a given threshold is itself a computing challenge. (A related problem is contour finding). I am myself not an expert, but I wonder what is the computing complexity of this step, and if this could quickly become prohibitively expensive when the dimensionality of outcome increases? Is it possible to take advantage of the smooth property of \\hat{p} to get more efficient solutions? How this step is achieved and visualized in the current examples and how this can be efficiently achieved in higher-dimensional examples is something worth more discussion. \n\nThe paper does not specify whether JAPAN, CONTRA, and PCP use the same training setup. For\na fully controlled comparison, it would be helpful to clarify whether the same flow architecture and\ntraining procedure were used across methods.\n\n Like most conformal methods, JAPAN does not guarantee conditional coverage in theory. But since JAPAN estimates conditional p(y|x) at its core,  perhaps some empirical studies can be added to show if it does reasonably well in conditional coverage?\n\n The method requires training normalizing flows and sampling extensively to approximate the region.\nRuntime and efficiency were not clearly discussed in the main text.\n\n\nPresentation-wise, the paper lacks consistency and rigoricity in symbol usage and mathematical languages. Most are minor issues, but worth a good round of proofreading. I am only giving three examples below, among other similar problems:\n\nIn the statement of Proposition 3, I think “z” was abused to denote both a sample of size N (z_1,\\cdots, z_n) and just one instance/ in h^{-1}(z, x)? Also, both z and z_i appeared in the Area estimation formula. And mathematically, it’s not tight to write z \\sim p_X(z), with the same symbol z on both sides…\n\nThe introduction, followed by a review of NF and inductive conformal prediction is generally well-done. The detailed description provided for ICP is actually for a most popular case of ICP, the split conformal. So that part could be a little misleading.\n\nSec 3.1 \\hat{p} was not officially defined before it appears inside the conformal prediction set.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3RMNP1kSS5", "reviewer": "Reviewer_LvS2", "summary": "The paper introduces JAPAN, a flow-based conformal prediction framework that uses density estimates (e.g., log p(y|x) from normalizing flows) as conformity scores to build compact, potentially disjoint prediction areas with finite-sample marginal coverage. It uses density-based conformity scores rather than residuals to avoid rigid geometries, enabling multimodal and context-adaptive regions.\nThe method constructs prediction areas by thresholding the estimated density and compute their volume efficiently via sampling in the flow’s latent space. It also provide theoretical motivation for efficiency under rank-preserving mappings and bounded misranking. The paper demonstrates versatility through variants: unconditional density, conditional on base predictions, posterior density, latent-space scoring, and a “tau-adaptive” variant to adjust thresholds with context. It conducts extensive experiments empirically evaluating across multivariate regression and time-series forecasting, reporting consistent calibration and smaller areas than strong baselines.", "questions": "Could you include a formal argument showing JAPAN in latent is equivalent to CONTRA, \nor CONTRA is a special case within JAPAN, and under what conditions?\n\nWhat is the variance of the area estimator (Proposition 3) ? This estimator appears unstable due to the appearance\nof the density function in the denominator.\n\nCould you please clarify the second condition in the proof of Proposition 1 (page 24), which states that the conditional densities \np(y∣x) are homogeneous across x? Additionally, could you provide some examples where this condition is satisfied?\n\nJAPAN appears to only provide marginal guarantee. What about conditional converage?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The strengths of the paper include: \n1. JAPAN’s “geometry-free” density-thresholded sets can be disjoint and adapt to x, and the framework accommodates several practical scoring variants (conditional, posterior, latent), which broadens applicability.\n2. Proposition 3 gives a tractable estimator via latent-space sampling and change-of-variables of the area/volume.\n3. The experimental results consistently show comparable coverage and significantly smaller areas than baselines like CQR, PCP, CONTRA, copula-based methods, etc.", "weaknesses": "The main weakness is that the theoretical assumptions in the paper are strong and somewhat informal: Propositions 1 and 2 rely on strict monotonic transformations and density homogeneity assumptions. Proposition 4 uses a misranking mass argument without explicit finite-sample rates or dimension dependence; these read more as plausibility results than tight guarantees. Also, the claimed equivalence to CONTRA is supported empirically but not formalized.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xPDpGBa6uX", "reviewer": "Reviewer_EZpP", "summary": "The author proposes JAPAN, a multidimensional conformal prediction method that is based on Normalising flows to efficiently construct prediction regions. The author argues that current multidimensional conformal prediction methods are usually imposing restrictive geometric constraints like rectangular or spherical which often fail to capture more complex, multimodal distributions, thus leading to overly conservative sets that lose efficiency. By training a flow model, it allows for complex shape of regions which focus on the high-density areas, thus greatly improving efficiency. Besides, the author provides theoretical guarantees showing that when the density estimate is correct, the resulting prediction area is optimal. Then the author validates the performance of JAPAN by comparing it with various benchmark conformal methods.", "questions": "1. I notice there is other conformal prediction method for multidimensional time series ''Flow-based Conformal Prediction for Multi-dimensional Time Series'', how would you compare their method with yours?\n\n2. It seems that the theory part is lacking the coverage guarantee. Would the author explain more in detail about this?\n\n3. While establishing the optimality of the prediction region, the author assumes that the flow model can approximate the true density. Would the author explain more about when the assumption might approximately hold?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Strength:\n1. Provide theoretical guarantee showing that the region is optimal when the models correctly estimate the distribution\n2. Comprehensive benchmarking in the experiment showing the efficiency and accuracy of the method", "weaknesses": "Weakness:\n1. Lack coverage guarantees which are important in conformal prediction\n2. Not comparing with other related flow-based conformal prediction method", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "CAU5iMfQrV", "reviewer": "Reviewer_AaZM", "summary": "The paper introduces a conformal prediction framework that constructs prediction sets in the output space using the log-likelihood from normalizing flows as the conformity score. This design avoids fixed geometric assumptions and enables multi-modal, disconnected prediction regions. The authors provide a rank-preserving efficiency bound and propose a latent space estimator for computing region area. Experiments on several regression and time-series datasets show smaller regions at comparable coverage.", "questions": "See above", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. This paper replaces residual based conformity scores with output space likelihoods from a normalizing flow is a simple yet powerful idea that removes restrictive geometric assumptions and enables nonconvex, multi-modal prediction regions.\n2. The work subsumes prior CP variants (residual, latent, posterior) as special cases, and the paper demonstrates several practical extensions including conditional andτ-adaptive versions.\n3. The experiments span both regression and time series tasks, using eight diverse datasets, and consistently show smaller prediction areas at target coverage.", "weaknesses": "1. The paper provides intuitive rank-preserving arguments but stops short of offering stronger formal results, there is no finite-sample or conditional coverage bound. Besides the author should incorporate more theoretical evidence.\n2. The efficiency claim hinges on the flow preserving likelihood ranking. There is no evidence that this holds across architectures, capacities, or under mild misspecification. A few diagnostic plots or ablations could make the claim much more convincing.\n3. Although the experiments cover many datasets, several baselines (notably RCP, MCQR, and CopulaCPTS) are numerically unstable or poorly tuned. As a result, it’s unclear whether the observed area reductions stem from the proposed idea or from implementation details.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)