Valid and Expressive Copulas for Irregular Multivariate Time Series

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Valid and Expressive Copulas for Irregular Multivariate Time Series

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

arXiv:2605.23632v1 [cs.LG] 22 May 2026

Valid and Expressive Copulas for Irregular Multivariate Time Series Christian Klötergens Institute of Computer Science University of Hildesheim Hildesheim, Germany [email protected]

Tom Hanika Institute of Computer Science University of Hildesheim Hildesheim, Germany

Lars Schmidt-Thieme Institute of Computer Science University of Hildesheim Hildesheim, Germany

Vijaya Krishna Yalavarthi Institute of Computer Science University of Hildesheim Hildesheim, Germany

Abstract We introduce CoPFITi, a copula model for probabilistic forecasting of irregular multivariate time series (IMTS). Our model combines the expressivity of normalizing flows for univariate marginals with the consistency and flexibility of a Gaussian Mixture Copula for the joint dependency structure. Our experiments show that copula-based approaches, which decouple the marginals from the joint, yield better marginal models than architectures that directly fit the full joint. With CoPFITi we propose the first IMTS copula that is marginalization consistent by construction, and establish a new state of the art in joint IMTS density modeling.1

Introduction

Sparse and irregularly sampled multivariate time series (IMTS) arise in many real-world domains, including healthcare, climate science, and sensor networks, where variables are observed asynchronously and at non-uniform timestamps. In these settings, predicting future values is inherently uncertain due to noise, missing observations, and partial observability of the underlying system. As a result, for many decision making tasks, accurate forecasting requires not just point predictions but full probabilistic forecasts. Crucially, these forecasts must capture the joint distribution over all queried variables and time points. Dependencies across channels and time carry important information, and modeling them independently as marginal distributions leads to incoherent and possibly misleading predictions. The irregular structure of IMTS introduces an additional challenge: predictions must be marginalization consistent [Yalavarthi et al., 2025b]. Since each query may involve a different subset of variables and timestamps, the model must ensure that predictions over any subset agree with those obtained from larger joint predictions. This requirement is difficult to satisfy unless the model admits tractable and well-defined marginal distributions. Copula models provide a natural way to address this challenge by separating univariate marginals from the multivariate dependency structure. This decomposition can guarantee marginalization consistency when it is properly constructed. It also enables specialization of both components. Marginals can be learned independently of the dependency model. This often leads to more accurate univariate predictions than joint training. At the same time, the dependency model can focus entirely 1 The code can be found here: https://anonymous.4open.science/r/CoPFITi-81E1

Preprint.

(a) Ground Truth

(b) Marginals

(d) GM-C

y1 4

−2

Density

−4

y1 4

−2

Density

−4

−4 −2

Figure 1: Demonstrating the advantage of a Gaussian Mixture Copula (GM-C) over a Mixture of Gaussian Copulas (Mix-GC). (b) shows the univariate marginal density functions of the 2-D distribution shown in (a). (c) and (d) depict attempts to model the distribution with Mix-GC and GM-C using the respective marginal distributions. For both models and datasets we used mixtures of 5 components.

on capturing interactions between variables, without needing to model marginal behavior. Currently, TACTiS [Drouin et al., 2022] and TACTiS-2 [Ashok et al., 2023] are the only copula-based approaches suitable for IMTS. However, they fail to strictly preserve the univariate marginals, which violates a defining property of copulas [Nelsen, 2006]. We propose CoPFITi (Copulas for Probabilistic Forecasting of Irregular Time Series), a copula model tailored to IMTS. Our approach constructs the dependency structure in a latent space using a Gaussian mixture model, while fully decoupling it from independently trained univariate marginals. This design ensures marginalization consistency by construction while remaining flexible and expressive. Our contributions are as follows: 1. We introduce CoPFITi, to the best of our knowledge, the first copula-based framework for IMTS. CoPFITi constructs the dependency structure in a latent space using a Gaussian mixture model, while fully decoupling it from the marginal distributions. 2. We propose MargFlow, a simple yet strong model for univariate marginal distributions using Deep Sigmoidal Flows. In our evaluation on established benchmark tasks on four datasets, we show that MargFlow achieves the best marginal likelihood. 3. CoPFITi attains joint likelihoods that are significantly better than non-copula baselines and are on par with or often better than TACTiS-2. In contrast to TACTiS-2, CoPFITi is consistent.

Background

2.1

Copulas

Throughout, uppercase letters denote cumulative distribution functions (CDFs) and the corresponding N lowercase letters denote their densities. A Copula is a CDF with domain [0, 1] and uniform marginals. Sklar’s theorem [Sklar, 1959, Nelsen, 2006] states that any joint distribution F with marginals Fn admits a copula C such that: F (y) = C F1 (y1 ), . . . , FN (yN ) . (1) Intuitively, copulas allow separating the dependency structure from the univariate marginal distributions. If the marginals are continuous, C is unique. If F is continuous with density f and copula 2

density c, then: N Y f (y) = c F1 (y1 ), . . . , FN (yN ) fn (yn ).

(2)

n=1

Copulas from an auxiliary density. Let g be a continuously differentiable joint density with marginals gn . Then cg denotes the density of the copula induced by g: −1 g G−1 N 1 (u1 ), . . . , GN (uN ) cg (u) = , u ∈ [0, 1] . (3) QN −1 n=1 gn Gn (un ) The corresponding copula depends only on the dependence structure of g, since its marginals are removed by the probability integral transform. A well-established instance of this construction is the Gaussian copula, where g is a multivariate Gaussian with standard normal marginals, i.e. zero mean and unit variances, so that its covariance matrix coincides with a correlation matrix Σ. Standardizing the marginals is convenient: the dependence structure is then fully captured by Σ, and the marginal transformations Gn and G−1 n reduce to the standard normal CDF and its inverse. However, the Gaussian copula can only represent the dependence structure of a multivariate Gaussian. Why a mixture of Gaussian copulas is not enough. A natural attempt to increase flexibility PK is to take a mixture of Gaussian copulas, i.e. a convex combination c(u) = j=1 πj cΣj (u) of Gaussian copulas with distinct correlation matrices Σj . Since each component cΣj is a valid copula, the resulting mixture remains a valid copula. This construction captures dependence structures that arise as a superposition of several linear-correlation regimes. For instance, the “X”-shaped density in Figure 1 (top row) can be well-approximated by a mixture of two Gaussians with oppositesign correlations. This construction is, however, fundamentally restricted: each component cΣj is induced by a zero-mean latent Gaussian, so all components share a common center in the latent space RN , and the Gaussian copula is by construction invariant to marginal location. Consequently, a mixture of Gaussian copulas can only interpolate between different linear dependence patterns sharing a common center, and cannot represent dependence structures induced by multimodal latent densities. Multi-modal structures such as the ring-shaped density in Figure 1, whose probability mass concentrates on a manifold away from the origin, therefore lie outside the expressive class of any mixture of Gaussian copulas, regardless of the number of components [Khaled and Kohn, 2023]. Gaussian Mixture Copulas. A Gaussian Mixture Copula (GM-C) [Bilgrau et al., 2016, Rajan and Bhattacharya, 2016, Tewari, 2023] is not a mixture of Gaussian copulas. GM-C first defines a Gaussian mixture model in RN , gθ (z) =

K X

πj N (z; µj , Σj ),

(4)

j=1

and maps un ∈ [0, 1] to the latent space via the marginal inverse CDFs: zn = G−1 θ,n (un ) where Gθ,n denotes the marginal CDF of the Gaussian mixture gθ . The copula is then obtained by applying (3) using the mixture density gθ . In contrast, a mixture of Gaussian copulas applies (3) to each Gaussian component gj separately, yielding copulas cgj , which are then combined via a convex sum. The PK difference is therefore the order of operations: j=1 πj cgj ̸= cPK . Mixing after the copula j=1 πj gj transformation forces all components to share the same marginal transformations G−1 n (one per component, applied independently), whereas mixing before the transformation lets the marginal transformations G−1 θ,n depend on the full mixture, so that the induced copula inherits the multimodal structure of gθ . This enables GM-Cs to represent complex, non-elliptical dependence patterns and, in principle, approximate a broad class of dependence structures arbitrarily well. 2.2

Probabilistic IMTS Forecasting ⊤

We aim to model the multivariate probability density of N many future values y = (y1 , . . . , yN ) ∈ RN of an irregular multivariate time series (IMTS). Due to the sporadic and irregular sampling across 3

both time and channels, this density is conditioned on two inputs: (i) a query Q specifying where, in future, predictions are required, and (ii) the observed history X . The forecasting distribution is thus: p(y | Q, X ) : RN → R≥0 .

(5)

N The query Q consists of N target locations, Q = (tn , cn ) n=1 , where (tn , cn ) ∈ R × {1, . . . , C}, with tn denoting a continuous timestamp and cn the channel index. The component yn corresponds to the target value at the n-th query location. The observed history is given as a collection of M M triplets, X = (tm , cm , ym ) m=1 , where (tm , cm , ym ) ∈ R × {1, . . . , C} × R. Here, tm is the observation time, cm the channel index, and ym the observed value. Note that for any two IMTS, query lengths and/or history lengths can vary. Furthermore, the representation does not impose any intrinsic ordering on the observations or query points. As a result, the sequences X and (Q, y) are naturally interpreted as sets rather than ordered lists, implying that any model operating on them should be permutation-invariant. A suitable model must (i) handle variable-sized inputs X and Q, (ii) capture the full joint distribution over y, and (iii) satisfy consistency conditions implied by the Kolmogorov extension theorem [Øksendal, 2003]. In particular, it should obey marginalization consistency: for any subset of indices S ⊆ {1, . . . , N } with complement −S, Z p(y | Q, X ) dyS = p(y−S | Q−S , X ). (6) This condition requires that predictions over any subset of query points remain consistent, regardless of whether additional query points are included or marginalized out [Yalavarthi et al., 2025b].

MargFlow: A Simple and Strong Model for Univariate Marginals

Constructing a copula-based model requires a dedicated model for the univariate marginals. Existing models [De Brouwer et al., 2019, Biloš et al., 2021, Schirmer et al., 2022] for marginals rely on Neural-ODEs and are restricted to Gaussian distributions. While Neural-ODEs have been shown to be ineffective in forecasting [Klötergens et al., 2024], limiting marginals to Gaussian is very restrictive. To fulfill the necessity of a competitive marginal model, we combine a state-of-the-art encoder with expressive Deep Sigmoidal Flows (DSF) [Huang et al., 2018] and name the resulting model MargFlow. Encoder. To parameterize the DSFs of MargFlow we need an encoding en for every element of the query (tn , cn ) ∈ Q: en = Enc((tn , cn ), X ) (7) Since the final joint model can only be marginalization consistent if every component is, MargFlow itself must satisfy this property. This in turn requires an encoder that fully separates the elements of a query Q, so that the predicted distribution at any query point depends only on the observation history X and that single query point, not on the rest of Q. The encoder of CircuITS [Klötergens, 2026] is built to do exactly this, as CircuITS itself is designed to be marginalization consistent. As CircuITS is the strongest currently available marginalization-consistent model on our benchmarks, we adopt its encoder unchanged for MargFlow. Normalizing Flows for marginal CDFs. We use the query encodings (en ) to model the respective marginal CDFs Fn via a DSF [Huang et al., 2018]. A DSF is a strictly monotone scalar transformation L Tθn : R → (0, 1) parameterized by θn = {a(ℓ) , b(ℓ) , w(ℓ) }ℓ=1 , defined as a composition of L sigmoidal blocks. A single block of width M acts on a scalar input x as ! M X −1 s(x; a, b, w) = σ wm σ(am x + bm ) (8) m=1

P σ is the logistic sigmoid, am > 0, bm ∈ R, and wm ≥ 0 with m wm = 1. Positivity of am and the convex combination over sigmoids ensure that the inner sum is a strictly increasing function valued in (0, 1), so applying σ −1 yields a strictly increasing map R → R. Stacking L such blocks gives a strictly increasing map R → R, and a final logistic sigmoid squashes the output to (0, 1), producing 4

ŷn MARGINAL LIKELIHOOD

MargFlow Encm

n fn (yn )

DSF un = Fn (yn )

CONDITIONAL JOINT LIKELIHOOD

X ×

CoPFITi

GMM

ûn = FGMM,n (ẑn )

ẑ Encc

p(y | Q, X )

Σ1 , . . . , Σ K µ1 , . . . , µ K

COPULA DENSITY

−1 zn = FGMM,n (un )

f (z) QN GMM n fGMM,n (zn )

Figure 2: CoPFITi architecture. → is used to indicate the computation of the joint likelihood. In contrast, we use 99K to indicate the sampling procedure. → represents the conditioning via the encoder, which is performed during likelihood estimation and sampling. a valid CDF. The conditional CDF at query point n is then Fn (yn ) = Tθn (yn ), where the block parameters θn are produced from the query encoding en via a small MLP θn = MLPDSF (en ) . The positivity of a(ℓ) is enforced by a softplus activation and the simplex constraint on w(ℓ) is enforced by a softmax. Because the map is strictly monotonic, Tθn is a one-dimensional normalizing flow that pushes the target yn forward to un = Fn (yn ) ∈ (0, 1), distributed uniformly under the model. The likelihood of a Normalizing Flow is obtained by the change of variable formula, ∂Tθn fn (yn ) = pU Tθn (yn ) · (9) (yn ) , ∂yn {z } | =1

where the first term vanishes because the base distribution pU is uniform on (0, 1), so the likelihood reduces to the Jacobian of Tθn . Conditioning θn only on en is what makes MargFlow query-separable, and thus marginalization consistent.

CoPFITi

We now introduce CoPFITi, our guaranteed-valid copula model for IMTS. Recall from Section 2.1 that a Gaussian Mixture Copula is induced by a single base joint density fGMM (z) = PK N j=1 πj N (z; µj , Σj ) on R , whose own one-dimensional marginal CDFs FGMM,n define the copula via Equation (3). Designing CoPFITi therefore reduces to one task: predicting the parameters K {πj , µj , Σj }j=1 of this latent Gaussian mixture, conditional on the observation history X and the query Q, in such a way that the resulting copula is marginalization consistent. Although we evaluate CoPFITi together with MargFlow from Section 3, the construction is agnostic to the specific choice of marginals and would work with any marginalization-consistent univariate model that exposes a tractable CDF and PDF. In Figure 2 we illustrate how CoPFITi and MargFlow work together to model the joint likelihood of IMTS forecasting queries and how to sample from the learned distributions. Encoder. CoPFITi reuses the encoder of Section 3 unchanged, both to keep the architecture small and, crucially, to inherit its query-separability property. However, CoPFITi and MargFlow apply two fully separated instances of this encoder with potentially different hyperparameters. As an intermediate step this encoder produces a global summary of the IMTS H̃ ∈ RC×D , where C represents the number of channels. We aggregate H̃ into a global summary h̄ ∈ RD obtained by pooling via linear attention with a single attention-query. 5

Mixture weights. The mixture weights π = (π1 , . . . , πK ) are inferred from this global summary by a 2-layer MLP followed by a softmax function: π = softmax MLPπ (h̄) . (10) Since π only depends on h̄, which only depends on the observations X, it is not informed about which subset of queries is requested. This is intended and helps us to guarantee that CoPFITi is marginalization consistent. Mean vectors. For each mixture component j ∈ {1, . . . , K}, the latent Gaussian has its own mean vector µj ∈ RN , with one entry per query point. We predict the n-th entry of µj from the corresponding query embedding en via a small 2-layer MLP whose output dimension is K, µ1 , µ2 , . . . , µK n = MLPµ (en ) ∈ RK . (11) This factorization is what gives CoPFITi its expressivity beyond mixtures of Gaussian copulas: because µj ̸= µj ′ in general, the components of the latent GMM can have spatially separated modes, and the induced copula can therefore represent multi-modal dependencies. Standard deviations. The diagonal entries of Σj correspond to the per-query, per-component K standard deviations σj,n . We predict the full collection {σj,n }j=1 for query point n from en using a separate 2-layer MLP with a softplus output activation to enforce positivity, σ1,n , . . . , σK,n = softplus MLPσ (en ) ∈ RK (12) >0 . Importantly, since each [µj ]n and [σj ]n depend only on en , removing query points other than n from Q leaves [µj ]n and [σj ]n unchanged, which is necessary for marginalization consistency. Covariance matrices. The off-diagonal entries of Σj encode the dependence structure between query points within component j, and are the most delicate part of the construction. We construct ⊤ each Σj ∈ RN ×N in two steps. Stacking the query embeddings into E = [e1 , . . . , eN ] ∈ RN ×D , N ×H a 2-layer MLP applied row-wise produces K per-component feature matrices Uj ∈ R . We form their Gram matrices, add an identity regularizer, and normalize to correlation matrices, Gj = Uj Uj⊤ + I,

−1/2

Rj = DGj

−1/2

Gj DGj ,

DGj = diag(Gj )

(13)

Adding I to the Gram matrix guarantees Gj to be positive definite. Finally, we scale Rj on both sides by the predicted standard deviations to obtain the final covariance matrix, Σj = Dj Rj Dj ,

Dj = diag(σj,1 , . . . , σj,N ).

(14)

2 By construction, Σj is symmetric positive definite and has the predicted variances σj,n on its diagonal. Crucially, the entry Σj,nm depends only on en and em , so dropping any other query point from Q simply removes the corresponding row and column of Σj , exactly as required when marginalizing a Gaussian.

Marginalization consistency. Combining the three observations above, πj does not depend on Q, while µj and Σj depend on Q only through the corresponding entries. Marginalizing out a subset of query points thus reduces to the standard marginalization rule for a multivariate Gaussian, and propagates through the mixture and Equation (3), so CoPFITi is marginalization consistent by construction. A formal proof is given in Section A. Training. Following previous work [Yalavarthi et al., 2025a,b], we optimize CoPFITi by minimizing the normalized joint negative log-likelihood (njNLL), defined in (15). For copula-based approaches the predicted joint log-density decomposes as in (16) (see (2)), with un = Fn (yn | Q, X ). Because −1 MargFlow provides un directly, we set zn = FGMM,n (un ) and evaluate the copula term as in (17), PK where fGMM (z) = j=1 πj N (z; µj , Σj ) is the latent Gaussian mixture density and fGMM,n is its 6

n-th one-dimensional marginal. LnjNLL (θ) =

1 |B|

−

(Q,X ,y)∈B

1 log p̂(y | Q, X ), |y|

log p̂(y | Q, X ) = log c(u1 , . . . , uN ) +

N X

log fn (yn | Q, X ),

(15)

(16)

n=1

log c(u1 , . . . , uN ) = log fGMM (z) −

N X

log fGMM,n (zn ).

(17)

n=1

We use a two-stage optimization scheme [Ashok et al., 2023]: First, we train MargFlow on its own by deactivating the copula and treating all query points as independent. Then, we freeze its weights and train the copula. In the second stage, the marginal term is constant, so minimizing njNLL is equivalent to maximizing the copula log-density. Evaluating the log copula density requires −1 zn = FGMM,n (un ), but the GMM marginal CDF cannot be inverted analytically. We obtain zn by a safeguarded Newton-bisection scheme in the forward pass (Section B), and compute gradients −1 of FGMM,n with respect to the GMM parameters analytically [Tewari, 2023]. Lemma 1 states the resulting closed forms, and Section C derives them step by step; in the backward pass we compose these expressions with the softmax and softplus parameterizations. −1 Lemma 1 (Gradients of the GMM inverse CDF). Let zn = FGMM,n (un ), and let Φ and φ denote the standard-normal CDF and PDF. Then the gradients of the GMM parameters are: zn −µj,n zn −µj,n zn −µj,n Φ π φ π φ (zn − µj,n ) j j σj,n σj,n σj,n ∂zn ∂zn ∂zn , , =− = = 2 3 f ∂πj fGMM,n (zn ) ∂µj,n σj,n fGMM,n (zn ) ∂σj,n 2 σj,n GMM,n (zn )

Sampling. Given X and Q, we first run the encoder to parameterize the GM-C. We then draw a sample ẑ ∈ RN from this mixture by first sampling a component index j ∼ Categorical(π) and then ẑ ∼ N (µj , Σj ). Next, we transform ẑ to the unit cube using the GMM’s own marginal CDFs FGMM,n , that is ûn = FGMM,n (ẑn ), where each FGMM,n is itself a one-dimensional Gaussian mixture CDF and is available in closed form. Finally, we invert the flow of the marginal model on each coordinate to obtain ŷn = Fn−1 (ûn ). Limitations. Each Σj has a rank-H off-diagonal Gram structure, so CoPFITi cannot represent near-full-rank dependence when H ≪ N . In practice this ceiling is rarely binding (real-world dependencies are rarely full-rank) and improves scaling in N (see Section E and Section F). Second, to guarantee marginalization consistency, the mixture weights π depend solely on the history X , not on the query Q.

Related Work

Deep Learning for probabilistic IMTS forecasting. The majority of the IMTS forecasting literature concerns simple point forecasting [Yalavarthi et al., 2024, Li et al., 2025, Luo et al., 2025, Zhang et al., 2024]. To address the need for uncertainty quantification, architectures such as GRUODE-Bayes [De Brouwer et al., 2019], Continuous Recurrent Units (CRU) [Schirmer et al., 2022], and Neural Flows [Biloš et al., 2021] modeled the forecasting targets as independent Gaussians. Recent Normalizing Flow [Papamakarios et al., 2021]-based architectures aim to capture the full joint distribution of an IMTS. ProFITi [Yalavarthi et al., 2025a] models complex joint dependencies but lacks marginalization consistency, frequently yielding contradictory forecasts across variable subsets. Conversely, MOSES [Yalavarthi et al., 2025b] guarantees consistency by applying separable flows over a latent multivariate Gaussian mixture model. CircuITS [Klötergens, 2026] addresses this expressivity trade-off by explicitly modeling inter-channel dependencies via probabilistic circuits [Choi et al., 2020] and intra-channel dynamics using Gaussian Copulas. While CircuITS is marginalization-consistent and theoretically a universal approximator, circuits are notoriously inefficient at parameterizing dense continuous dependencies (such as linear correlations). 7

Copulas for Time Series. Copulas are a foundational tool for modeling time series dependencies in finance and classical statistics [Patton, 2012, Größer and Okhrin, 2022]. However, these applications typically assume a single, global distribution across the dataset, which does not hold in our domain. Machine Learning-based copulas for time series have historically been centered around regular time series and were restricted to Gaussian Copulas [Wilson and Ghahramani, 2010, Salinas et al., 2019, Wen and Torkkola, 2019]. TACTiS [Drouin et al., 2022] and TACTiS-2 [Ashok et al., 2023] attempt to address this gap by presenting an attentional copula that can be applied to irregular time series. However, these models do not preserve the marginals by architectural design. Instead, they rely entirely on convergence theory to learn this property.

Experiments

We empirically evaluate MargFlow and CoPFITi against state-of-the-art IMTS forecasting baselines, examining marginal likelihood, joint likelihood, the contribution of the GM-C, the comparison against TACTiS-2, and the validity gap that follows from CoPFITi’s marginalization-consistent construction. Datasets. Our evaluation uses the four benchmark IMTS datasets that appear throughout the paper: USHCN [Menne et al., 2016], PhysioNet-2012 [Silva et al., 2012], MIMIC-III [Johnson et al., 2016], and MIMIC-IV [Johnson et al., 2023]. Crucially, the well-established preprocessing, binning, and split protocol is adopted directly from prior work [Biloš et al., 2021, Yalavarthi et al., 2025a,b, Klötergens, 2026] to provide a fair comparison. This protocol trains and evaluates each model five times on different train, validation, and test splits. Importantly, every model is evaluated on the same five splits, so the runs are paired across methods. We refer to Section D for more details. Baselines. Our comparison covers both univariate and joint baselines. On the univariate side, the baselines are GRU-ODE [De Brouwer et al., 2019], NeuralFlows [Biloš et al., 2021], and CRU [Schirmer et al., 2022]. On the joint side, we use ProFITi [Yalavarthi et al., 2025a], Gaussian Process Regression [Dürichen et al., 2015], MOSES [Yalavarthi et al., 2025b], and CircuITS [Klötergens, 2026]. As a copula baseline we use TACTiS-2 [Ashok et al., 2023]. Since it was not originally designed for sparse IMTS, we adapt it to our setting; we describe this adaptation in detail in Section D.8. Details on hyperparameters are given in Section D.3. Because we adopt the identical protocol used in prior work, baseline numbers for GRU-ODE, NeuralFlows, CRU, GPR, ProFITi, MOSES, and CircuITS are taken directly from Yalavarthi et al. [2025a,b], Klötergens [2026].

Model

USHCN

Physionet

MIMIC-III

MIMIC-IV

Joint

ProFITi GPR MOSES CircuITS Joint-Abl

-3.324 ± 0.206 1.235 ± 0.096 -3.355 ± 0.156 -3.717 ± 0.201 -3.663 ± 0.210

-0.016 ± 0.085 1.161 ± 0.065 -0.271 ± 0.028 -0.287 ± 0.037 -0.176 ± 0.470

0.408 ± 0.030 1.341 ± 0.009 0.163 ± 0.026 0.095 ± 0.107 0.108 ± 0.238

0.500 ± 0.322 1.161 ± 0.010 -0.634 ± 0.017 -0.731 ± 0.063 -0.677 ± 0.044

Univariate

Table 1: Comparing mNLL (Marginal Negative Log-Likelihood) across datasets. Lower values indicate better performance. We present the mean and standard deviation of 5 runs. The best model is marked in bold.

GRU-ODE NeuralFlows CRU MargFlow

0.776 ± 0.172 0.775 ± 0.180 0.762 ± 0.180 -3.948 ± 0.294

0.504 ± 0.061 0.492 ± 0.029 0.931 ± 0.019 -0.479 ± 0.021

0.839 ± 0.030 0.866 ± 0.097 1.209 ± 0.044 -0.158 ± 0.205

0.876 ± 0.589 1.796 ± 0.050 OOM -0.833 ± 0.077

MargFlow predicts better marginals. Table 1 reports the marginal Negative Log-Likelihood (mNLL, defined in (38), Section D.5) of MargFlow alongside the baselines. On all four datasets, MargFlow has the lowest mNLL. To assess the effect of the isolated marginal training, we also report the Joint-Ablation (Joint-Abl), which uses the same architecture as CoPFITi but trains all weights jointly under njNLL. We want to highlight that the mNLL of Joint-Abl is significantly worse than that of MargFlow trained in isolation and on par with the best joint baselines MOSES and CircuITS. We present learning curves of Joint-Abl and the MargFlow + CoPFITi setup in Section G. 8

Table 2: Comparing njNLL (Normalized Joint Negative Log-Likelihood) across datasets. Lower values indicate better performance. We present the mean and standard deviation of 5 runs. The best model is marked in bold. USHCN

Physionet

MIMIC-III

MIMIC-IV

Univariate

GRU-ODE NeuralFlows CRU MargFlow

0.766 ± 0.159 0.775 ± 0.152 0.761 ± 0.191 -3.969 ± 0.305

0.501 ± 0.001 0.496 ± 0.000 1.057 ± 0.007 -0.442 ± 0.021

0.961 ± 0.064 0.998 ± 0.111 1.234 ± 0.076 -0.505 ± 0.208

0.823 ± 0.318 0.689 ± 0.087 OOM -1.819 ± 0.065

Joint

ProFITi GPR MOSES CircuITS Joint-Abl

-3.226 ± 0.225 2.011 ± 1.376 -3.357 ± 0.176 -3.789 ± 0.218 -3.907 ± 0.269

-0.647 ± 0.078 1.367 ± 0.074 -0.491 ± 0.041 -0.550 ± 0.013 -0.489 ± 0.399

-0.377 ± 0.032 3.146 ± 0.359 -0.305 ± 0.027 -0.574 ± 0.080 -0.508 ± 0.309

-1.777 ± 0.066 2.789 ± 0.057 -1.668 ± 0.097 -2.113 ± 0.044 -2.019 ± 0.049

Copula

Model

TACTiS-2 CoPFITi (Mix-GC) CoPFITi

-4.254 ± 0.294 -4.078 ± 0.308 -4.135 ± 0.304

-0.732 ± 0.115 -0.728 ± 0.018 -0.745 ± 0.019

-0.780 ± 0.282 -0.779 ± 0.216 -0.835 ± 0.218

-2.053 ± 0.170 -2.080 ± 0.064 -2.201 ± 0.058

Copulas with good marginals outperform jointly trained models. Equipped with MargFlow’s marginals, the copula models consistently outperform all jointly trained baselines (including JointAbl) in terms of joint likelihood (Table 2). The sole exception occurs on MIMIC-IV, where CircuITS marginally edges out two of the three copula variants; nevertheless, CoPFITi retains a strictly superior marginal likelihood even in this case. CoPFITi can utilize the expressivity of GM-C. We investigate the impact of using a GM-C over a Mix-GC in two ways. First, we introduce the ablation CoPFITi-Mix-GC, in which the GM-C is replaced by a Mix-GC. We refer to Section D.6 for more details. Second, we compare CoPFITi against TACTiS-2. Comparing two copulas based only on the standard deviation over the five runs is insufficient, because the evaluation protocol uses different splits for each run. In addition, the performance of MargFlow varies between runs. Hence, the paired Corrected Resampled ttest [Nadeau and Bengio, 1999] shown in Table 3 provides more insight into significance. CoPFITi has a p-value < 0.05 versus CoPFITi-Mix-GC on all datasets but USHCN, where it is close to 0.05. Compared to TACTiS-2, CoPFITi achieves a lower mean njNLL on PhysioNet-2012, MIMIC-III, and MIMIC-IV, though the t-test cannot rule out random noise on any of these three datasets. Table 3: Bounds of p-values of corrected resampled Student’s t-test [Nadeau and Bengio, 1999] across datasets. Model

USHCN

Physionet

MIMIC-III

MIMIC-IV

CoPFITi / CoPFITi-Mix-GC CoPFITi / TACTiS-2

0.0538 0.0006

0.0017 0.8563

0.0006 0.3759

0.0000 0.2255

Quantifying that CoPFITi is more valid than TACTiS-2. To quantify the validity gap, we draw 1000 samples from MargFlow and from the full copula setup, and compare the per-dimension univariate Wasserstein-1 distance between the resulting empirical copulas (see Section D.7 for the definition). Lower values indicate a closer match in between distributions. The Control row is not a separate copula model: it is obtained by sampling twice from the same MargFlow instance using different random seeds. Therefore, it represents the irreducible sampling error. As shown in Table 4, CoPFITi is either within standard deviations of the control or very close to it. The WD of TACTiS-2, however, is orders of magnitude higher, indicating that it did not learn to perfectly preserve the marginal distributions. The experiment above probes only the univariate marginalization consistency, which is the defining property of a copula. Due to the page limit, we refer to Section H for a more in-depth analysis of TACTiS-2’s marginalization inconsistency. 9

Table 4: Wasserstein distances ×10−3 for 1000 samples from MargFlow (alone) and the MargFlow + copula set up averaged over 5 splits. The Control row quantifies the sampling error calculated based on sampling twice from MargFlow with different random seeds. Model

USHCN

Physionet

MIMIC-III

MIMIC-IV

TACTiS-2 CoPFITi (GM-C)

27.60 ± 14.47 5.083 ± 2.444

10.78 ± 2.978 0.535 ± 0.131

15.72 ± 8.570 0.524 ± 0.079

2.939 ± 1.133 0.381 ± 0.081

Control

5.485 ± 2.541

0.633 ± 0.142

0.408 ± 0.043

0.274 ± 0.022

Conclusion

We introduced CoPFITi, the first marginalization-consistent copula model for irregular multivariate time series. By inferring a Gaussian Mixture Copula through a query-separable encoder that predicts mixture weights, means, and covariance matrices in a way that commutes with marginalization, CoPFITi matches the expressivity of attentional copulas while providing structural validity guarantees. MargFlow, a specialized DSF-based marginal model, produces the most accurate marginals on every benchmark, and feeding these marginals into copulas yields accurate estimation of the joint likelihood. Looking ahead, CoPFITi naturally extends to other conditional joint density estimation problems with variable target sets, most notably missing-value imputation in tabular data.

References Arjun Ashok, Étienne Marcotte, Valentina Zantedeschi, Nicolas Chapados, and Alexandre Drouin. TACTiS-2: Better, Faster, Simpler Attentional Copulas for Multivariate Time Series. In The Twelfth International Conference on Learning Representations, October 2023. Anders Ellern Bilgrau, Poul Svante Eriksen, Jakob Gulddahl Rasmussen, Hans Erik Johnsen, Karen Dybkaer, and Martin Boegsted. GMCM: Unsupervised Clustering and Meta-Analysis Using Gaussian Mixture Copula Models. Journal of Statistical Software, 70:1–23, April 2016. ISSN 1548-7660. doi: 10.18637/jss.v070.i02. Marin Biloš, Johanna Sommer, Syama Sundar Rangapuram, Tim Januschowski, and Stephan Günnemann. Neural Flows: Efficient Alternative to Neural ODEs. In Advances in Neural Information Processing Systems, volume 34, pages 21325–21337. Curran Associates, Inc., 2021. Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu. Recurrent Neural Networks for Multivariate Time Series with Missing Values. Scientific Reports, 8(1):6085, April 2018. ISSN 2045-2322. doi: 10.1038/s41598-018-24271-9. Y Choi, Antonio Vergari, and Guy Van den Broeck. Probabilistic circuits: A unifying framework for tractable probabilistic models. UCLA. URL: http://starai. cs. ucla. edu/papers/ProbCirc20. pdf, 6, 2020. Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. GRU-ODE-Bayes: Continuous Modeling of Sporadically-Observed Time Series. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. Alexandre Drouin, Étienne Marcotte, and Nicolas Chapados. TACTiS: Transformer-Attentional Copulas for Time Series. In Proceedings of the 39th International Conference on Machine Learning, pages 5447–5493. PMLR, June 2022. Robert Dürichen, Marco A. F. Pimentel, Lei Clifton, Achim Schweikard, and David A. Clifton. Multitask Gaussian Processes for Multivariate Physiological Time-Series Analysis. IEEE Transactions on Biomedical Engineering, 62(1):314–322, January 2015. ISSN 1558-2531. doi: 10.1109/TBME.2014.2351376. Tilmann Gneiting and Adrian E Raftery. Strictly Proper Scoring Rules, Prediction, and Estimation. Journal of the American Statistical Association, 102(477):359–378, March 2007. ISSN 0162-1459. doi: 10.1198/ 016214506000001437. Joshua Größer and Ostap Okhrin. Copulae: An overview and recent developments. WIREs Computational Statistics, 14(3):e1557, 2022. ISSN 1939-0068. doi: 10.1002/wics.1557. Chin-Wei Huang, David Krueger, Alexandre Lacoste, and Aaron Courville. Neural Autoregressive Flows. In Proceedings of the 35th International Conference on Machine Learning, pages 2078–2087. PMLR, July 2018. Alistair E. W. Johnson, Tom J. Pollard, Lu Shen, Li-wei H. Lehman, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G. Mark. MIMIC-III, a freely accessible critical care database. Scientific Data, 3(1):160035, May 2016. ISSN 2052-4463. doi: 10.1038/sdata.2016.35. Alistair E. W. Johnson, Lucas Bulgarelli, Lu Shen, Alvin Gayles, Ayad Shammout, Steven Horng, Tom J. Pollard, Sicheng Hao, Benjamin Moody, Brian Gow, Li-wei H. Lehman, Leo A. Celi, and Roger G. Mark. MIMIC-IV, a freely accessible electronic health record dataset. Scientific Data, 10(1):1, January 2023. ISSN 2052-4463. doi: 10.1038/s41597-022-01899-x. Mohamad A. Khaled and Robert Kohn. On approximating copulas by finite mixtures, February 2023. Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization, January 2017. Christian Klötergens, Vijaya Krishna Yalavarthi, Randolf Scholz, Maximilian Stubbemann, Stefan Born, and Lars Schmidt-Thieme. Physiome-ODE: A Benchmark for Irregularly Sampled Multivariate Time-Series Forecasting Based on Biological ODEs. In The Thirteenth International Conference on Learning Representations, October 2024. Christian Klötergens. Probabilistic Circuits for Irregular Multivariate Time Series Forecasting, January 2026. URL http://arxiv.org/abs/2604.27814. Boyuan Li, Yicheng Luo, Zhen Liu, Junhao Zheng, Jianming Lv, and Qianli Ma. HyperIMTS: Hypergraph Neural Network for Irregular Multivariate Time Series Forecasting. In Forty-Second International Conference on Machine Learning, June 2025. Ilya Loshchilov and Frank Hutter. Decoupled Weight Decay Regularization, January 2019.

Yicheng Luo, Bowen Zhang, Zhen Liu, and Qianli Ma. Hi-Patch: Hierarchical Patch GNN for Irregular Multivariate Time Series. In Forty-Second International Conference on Machine Learning, June 2025. M. J. Menne, Jr Williams, and R. S. Vose. Long-Term Daily and Monthly Climate Records from Stations Across the Contiguous United States (U.S. Historical Climatology Network). Technical Report osti:1394920; cdiac:NDP-019; doi:10.3334/CDIAC/CLI.NDP019, Environmental System Science Data Infrastructure for a Virtual Ecosystem (ESS-DIVE) (United States); CDIAC, January 2016. Claude Nadeau and Yoshua Bengio. Inference for the Generalization Error. In Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999. Roger B. Nelsen. An Introduction to Copulas. Springer Series in Statistics. Springer, New York, NY, 2006. ISBN 978-0-387-28659-4. doi: 10.1007/0-387-28678-0. Bernt Øksendal. Stochastic Differential Equations. Universitext. Springer, Berlin, Heidelberg, 2003. ISBN 978-3-540-04758-2 978-3-642-14394-6. doi: 10.1007/978-3-642-14394-6. George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing Flows for Probabilistic Modeling and Inference. Journal of Machine Learning Research, 22(57): 1–64, 2021. ISSN 1533-7928. Andrew J. Patton. A review of copula models for economic time series. Journal of Multivariate Analysis, 110: 4–18, September 2012. ISSN 0047-259X. doi: 10.1016/j.jmva.2012.02.021. William H. Press. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press, September 2007. ISBN 978-0-521-88068-8. Vaibhav Rajan and Sakyajit Bhattacharya. Dependency clustering of mixed data with Gaussian mixture copulas. In Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, IJCAI’16, pages 1967–1973, New York, New York, USA, July 2016. AAAI Press. ISBN 978-1-57735-770-4. David Salinas, Michael Bohlke-Schneider, Laurent Callot, Roberto Medico, and Jan Gasthaus. High-dimensional multivariate forecasting with low-rank Gaussian Copula Processes. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. Mona Schirmer, Mazin Eltayeb, Stefan Lessmann, and Maja Rudolph. Modeling Irregular Time Series with Continuous Recurrent Units. In Proceedings of the 39th International Conference on Machine Learning, pages 19388–19405. PMLR, June 2022. Ikaro Silva, George Moody, Daniel J Scott, Leo A Celi, and Roger G Mark. Predicting in-hospital mortality of ICU patients: The PhysioNet/Computing in cardiology challenge 2012. In 2012 Computing in Cardiology, pages 245–248, September 2012. M. Sklar. Fonctions de répartition à N dimensions et leurs marges. Annales de l’ISUP, VIII(3):229, 1959. Ashutosh Tewari. On the Estimation of Gaussian Mixture Copula Models. In Proceedings of the 40th International Conference on Machine Learning, pages 34090–34104. PMLR, July 2023. Ruofeng Wen and Kari Torkkola. Deep Generative Quantile-Copula Models for Probabilistic Forecasting, July 2019. Andrew G Wilson and Zoubin Ghahramani. Copula Processes. In Advances in Neural Information Processing Systems, volume 23. Curran Associates, Inc., 2010. Vijaya Krishna Yalavarthi, Kiran Madhusudhanan, Randolf Scholz, Nourhan Ahmed, Johannes Burchert, Shayan Jawed, Stefan Born, and Lars Schmidt-Thieme. GraFITi: Graphs for Forecasting Irregularly Sampled Time Series. Proceedings of the AAAI Conference on Artificial Intelligence, 38(15):16255–16263, March 2024. ISSN 2374-3468. doi: 10.1609/aaai.v38i15.29560. Vijaya Krishna Yalavarthi, Randolf Scholz, Stefan Born, and Lars Schmidt-Thieme. Probabilistic Forecasting of Irregularly Sampled Time Series with Missing Values via Conditional Normalizing Flows. Proceedings of the AAAI Conference on Artificial Intelligence, 39(20):21877–21885, April 2025a. ISSN 2374-3468. doi: 10.1609/aaai.v39i20.35494. Vijaya Krishna Yalavarthi, Randolf Scholz, Christian Klötergens, Kiran Madhusudhanan, Stefan Born, and Lars Schmidt-Thieme. Reliable Probabilistic Forecasting of Irregular Time Series through MarginalizationConsistent Flows. In The Fourteenth International Conference on Learning Representations, October 2025b. Weijia Zhang, Chenlong Yin, Hao Liu, Xiaofang Zhou, and Hui Xiong. Irregular Multivariate Time Series Forecasting: A Transformable Patching Graph Neural Networks Approach. In Forty-First International Conference on Machine Learning, June 2024.

Marginalization Consistency of CoPFITi

In this appendix we give a formal proof that CoPFITi is marginalization consistent by construction. The argument reuses the framework of separable flows and mixtures of separable flows introduced by Yalavarthi et al. [2025b]. We first restate the two lemmas of that paper that we will invoke, adapted to the notation of the present work, and then apply them to CoPFITi. N

We adopt the notation of Section 4: X denotes the observation history, Q = {(tn , cn )}n=1 the query set, Qn = (tn , cn ) the n-th query point, en the per-query embedding produced by the MargFlow encoder, Fn the MargFlow marginal CDF at query n, and FGMM,n the n-th one-dimensional marginal of the latent Gaussian mixture density fGMM . We say that a conditional density p̂(y | Q, X ) on R|Q| is marginalization consistent if for every index k ∈ {1, . . . , |Q|}, Z p̂(y | Q, X ) dyk = p̂(y−k | Q−k , X ), (18) R

where Q−k and y−k denote Q and y with the k-th entry removed; the property extends to arbitrary subsets by induction. A.1

Cited Lemmas from Yalavarthi et al. (2025b)

Lemma 2 (Separable flows preserve marginalization consistency; Yalavarthi et al., 2025b, Lemma 3.1). Let f ( · | Q, X ) : R|Q| → R|Q| be a conditional flow that is separable, i.e. of the form f (z | Q, X ) = ϕ(z1 | Q1 , X ), . . . , ϕ(z|Q| | Q|Q| , X ) (19) for some univariate function ϕ that is invertible in its first argument. If the base density p̂Z (z | Q, X ) is marginalization consistent, then the pushforward density ∂f −1 (y | Q, X ) p̂(y | Q, X ) = p̂Z f −1 (y | Q, X ) Q, X · det ∂y

(20)

is also marginalization consistent. Lemma 3 (Mixtures with query-independent weights preserve marginalization consistency; Yalavarthi et al., 2025b, Lemma 3.2). Let p̂1 , . . . , p̂D be conditional densities that are each marginalization consistent, and let w : Seq(X ) → ∆D be a weight function that depends only on the observation history X . Then the mixture p̂(y | Q, X ) =

D X

wd (X ) p̂d (y | Q, X )

(21)

d=1

is marginalization consistent. For completeness, we recall the proof idea of Lemma 2: separability makes the Jacobian of f −1 diagonal, so its determinant factors across coordinates, the integral over yk collapses to an integral over zk via the change-of-variables theorem, and marginalization consistency of p̂Z then yields the claim. Lemma 3 follows by interchanging the integral over yk with the finite sum over mixture components, which is permitted because wd does not depend on Q. Full proofs are given in Appendices A.1 and A.2 of Yalavarthi et al. [2025b]. A.2

Main Result

Theorem 4. CoPFITi is marginalization consistent: for every observation history X , every query set Q, and every n ∈ {1, . . . , N }, the joint density predicted by CoPFITi satisfies (18). Proof. The proof proceeds in three steps. 13

Step 1: CoPFITi is a query-separable transformation of a latent base density. The generative process of CoPFITi maps a latent vector z, drawn from the latent Gaussian mixture base density fGMM (z | Q, X ), to the target domain via the per-coordinate transformation yn = Fn−1 FGMM,n (zn ) . (22) Both Fn and FGMM,n are parameterized entirely by the per-query embedding en , which by construction depends only on X and on the single query point Qn = (tn , cn ). Defining ϕ(zn | Qn , X ) := Fn−1 (FGMM,n (zn )), the joint transformation has the separable form required by Lemma 2. By that lemma, marginalization consistency of CoPFITi reduces to marginalization consistency of the latent base density fGMM . Step 2: Each Gaussian component of the latent mixture is marginalization consistent. The latent base density is the Gaussian mixture fGMM (z | Q, X ) =

K X

πj N (z; µj , Σj ).

(23)

j=1

As established in Section 4, the entry [µj ]n depends only on en , and the entry Σj,nm depends only on en and em . Removing a subset of query points from Q therefore deletes the corresponding entries of µj and the corresponding rows and columns of Σj , leaving the remaining entries unchanged. This coincides with the standard marginalization rule for a multivariate Gaussian, so each component N (z; µj , Σj ) is marginalization consistent in the sense of (18). Step 3: The mixture weights are query-independent, so the latent GMM is marginalization consistent. The mixture weights π = (π1 , . . . , πK ) are inferred from a global summary vector h̄ that is a function of the observation history X alone and is uninformed of the query Q. Together with Step 2, the two hypotheses of Lemma 3 are satisfied, so fGMM is marginalization consistent. Combining Steps 1 to 3, the latent base density of CoPFITi is marginalization consistent and is pushed forward to the target domain by a query-separable transformation. Lemma 2 then yields that the full joint density predicted by CoPFITi is marginalization consistent.

Numerical inversion of the GMM marginal CDF

−1 This appendix describes the numerical procedure we use to evaluate zn = FGMM,n (un ) in the forward pass of the copula training stage (Section 4). The marginal CDF of a one-dimensional Gaussian mixture, K X z−µ FGMM,n (z) = πj Φ σj,nj,n , (24) j=1

is strictly increasing and smooth in z, but its inverse has no closed form. We therefore solve FGMM,n (z) = u for z numerically. The key requirements are (i) batched evaluation across all query points and mini-batch elements simultaneously on GPU, (ii) robustness in the extreme tails where the PDF is near zero, and (iii) sufficient accuracy so that the analytic gradients of Lemma 1, which are evaluated at the returned z, are reliable. A pure Newton iteration meets (i) and (iii) when initialized well, but can overshoot or stall in the tails; a pure bisection is robust but converges only linearly. We combine the two into a safeguarded Newton-bisection solver, which is a standard choice for monotone scalar root finding [Press, 2007]. Initial bracket. Because each component is Gaussian, the probability mass of FGMM,n outside the interval min µj,n − 10 max σj,n , max µj,n + 10 max σj,n (25) j

−24

is below Φ(−10) ≈ 7.6×10 , well under floating-point precision. We use (25) as the initial bracket [LOW, HIGH] for every u ∈ (0, 1) and initialize the iterate at its midpoint x ← (LOW + HIGH)/2. Throughout the iteration we maintain the invariant FGMM,n (LOW) ≤ u ≤ FGMM,n (HIGH), so the true root is always contained in the bracket. 14

Iteration. At each step we evaluate the CDF and PDF at the current iterate x, both available in closed form from the GMM parameters, and perform the following updates, all batched along the query and mini-batch axes: 1. Bracket update. For each element, if FGMM,n (x) < u then the root lies above x, so we set LOW ← x; otherwise we set HIGH ← x. This preserves the bracketing invariant. 2. Newton step. Using F ′ (z) = fGMM,n (z), the Newton update is xnew = x −

FGMM,n (x) − u . fGMM,n (x)

The denominator is clamped from below by 10−8 to avoid division by zero in regions where the mixture density is numerically negligible. 3. Safeguard. We accept the Newton step only if it falls inside the current bracket, xnew ∈ [LOW, HIGH]. Otherwise (the Newton step overshot, or the local gradient was too flat for Newton to be informative), we fall back to a bisection step x ← (LOW + HIGH)/2 for that element. Combining the two rules in a single torch.where keeps the update fully vectorized: every element either takes the Newton step or a bisection step, depending on whether its Newton proposal was in-bracket. Convergence. Inside the bracket, Newton’s method on a smooth, strictly increasing F converges quadratically, while the bisection fallback contracts the bracket by a factor of two whenever Newton would have left it. The number of iterations is fixed (we use the same budget for every element of the batch, which keeps the kernel launches static and the whole loop compilable), and in practice a small number of steps is enough to drive |FGMM,n (x) − u| below the precision required by the −1 downstream gradients. Crucially, because the gradients of FGMM,n are obtained analytically via implicit differentiation (Section C), we do not need to backpropagate through the iterations: the loop is wrapped in a no_grad context, and only the converged z is fed into the closed-form expressions of Lemma 1. Why this matters for training. Treating the solver as a black box that returns z and pairing it with the analytic gradients of Lemma 1 has three consequences. First, the memory footprint of the backward pass is independent of the number of Newton iterations, since none of the intermediate iterates are retained on the autograd tape. Second, the safeguarded update means the solver does not diverge in the tails, where the GMM density is small and a naive Newton step would otherwise overshoot far outside the support of the marginals. Third, because the bracket (25) is constructed directly from the predicted GMM parameters, no dataset-dependent hyperparameters are introduced.

Derivation of the analytic iCDF gradients

This appendix is intended to be self-contained. Our goal is to compute, by hand, the gradients of the −1 marginal inverse CDF zn = Fθ,n (un ) with respect to the mixture parameters θ. We do this because the inverse CDF has no closed form: at training time, zn is obtained numerically (via Newton’s method). Naively back-propagating through that solver would require unrolling all of its iterations, which is expensive and numerically fragile. The classical alternative, implicit differentiation, sidesteps the solver entirely: it expresses the gradient of zn directly in terms of quantities we already have at the converged solution. The result is a handful of simple algebraic expressions, given in Proposition 1 of Section 4. The derivation below shows where each of those expressions comes from. What does the inverse CDF do? At a high level, Fθ−1 takes a number u ∈ (0, 1), interpreted as a probability, and returns the value z at which a Gaussian-mixture random variable accumulates exactly that much probability mass below it. That is, z is defined implicitly by Fθ (z) = u,

Fθ (z) = Pr[Z ≤ z | θ].

There is no algebraic formula that solves this equation for a general mixture, so we use Newton’s method to find z numerically. But once z is found, the equation Fθ (z) = u is an exact identity that we can differentiate symbolically. That is the whole trick. 15

Which partials do we need? The pseudo-observation un is treated as an input to the inverse CDF and is upstream of the mixture parameters; in our backward pass we do not need a gradient with 2 respect to it. The mixture parameters θ = (πj , µj,n , σj,n ) are predicted by the marginal network, 2 so the partial derivatives we actually need are ∂zn /∂πj , ∂zn /∂µj,n , ∂zn /∂σj,n . This restriction makes the derivation a bit cleaner, since u is treated as a constant whenever we differentiate the defining identity. Throughout the derivation n is fixed and we drop the subscript n wherever doing so does not cause 2 confusion: write z = zn , u = un , µj = µj,n , σj = σj,n , σj2 = σj,n . The marginal CDF and PDF are then K X z − µj Fθ (z) = πj Φ tj (z) , tj (z) := , (26) σj j=1 fθ (z) =

K X dFθ πj (z) = φ tj (z) , dz σ j=1 j

(27)

where Φ and φ are the standard-normal CDF and PDF. The variable tj (z) is just the standardized distance of z from the centre of component j, measured in units of that component’s standard deviation, this notation will keep the algebra below compact. P A note on the simplex constraint. The mixture weights satisfy πj > 0 and j πj = 1, but we do not enforce the simplex constraint inside this derivation: we treat each πj as a free positive scalar in (26), take partial derivatives with the other πj ′ held fixed, and let the softmax (or whichever projection produces π) handle the constraint outside. This is exactly what the code does: PyTorch’s autograd composes the softmax Jacobian with ∂z/∂πj during back-propagation, so working with unconstrained πj here gives the correct end-to-end gradient once it is composed with the constraint-respecting layer. C.1

Setup: the implicit equation and its differential

The forward pass solves Fθ (z) = u (28) for z = z(θ, u). Equation (28) is the defining identity: at the converged Newton iterate (up to the chosen tolerance, which we take to be machine precision in float32) it holds as a true equality. This is the only place in the derivation where we use that Newton’s method returned a true root: we needed (28) to be an identity in order to differentiate it. The trick of implicit differentiation is to take the derivative of both sides of an identity and read off the gradient of the implicit quantity (here z) without ever solving the equation symbolically. Concretely, we differentiate (28) with respect to a generic mixture parameter ξ ∈ {πj , µj , σj2 }. Since u does not depend on ξ, the right-hand side is zero: ∂ Fθ z(θ, u) = 0, ∂ξ ∂Fθ ∂z ∂Fθ (z) + = 0. ∂z ∂ξ ∂ξ z fixed | {z }

(29)

fθ (z)

The second line uses the chain rule. The variable ξ enters Fθ (z) in two ways: through the explicit dependence of Fθ on its parameters θ at a fixed argument z, and through the implicit dependence of the argument z itself on ξ (because changing θ changes the root of (28)). The first term collects the implicit channel; the second term collects the explicit channel. Solving (29) for ∂z/∂ξ gives the master formula we will specialize three times: ∂z 1 ∂Fθ = − . ∂ξ fθ (z) ∂ξ z fixed

(30)

Two remarks make this formula useful in practice. First, the denominator fθ (z) > 0 everywhere, because each Gaussian component has full support on the real line, so the formula is well-defined and 16

never blows up in finite precision provided the σj are bounded away from zero (which is enforced by the softplus parametrization in our network). Second, every term on the right-hand side is something we already know how to compute at the converged z, the value of the mixture PDF and the explicit partial of the mixture CDF, both available in closed form. So this single line replaces the entire backward pass through the Newton solver. C.2

Specializing for each parameter

We now plug each of the three mixture parameters into (30). The pattern is always the same: compute the explicit partial of Fθ at fixed z and divide by −fθ (z). (a) ξ = πP j (gradient with respect to a mixture weight). Looking at (26), only the j-th term in the sum j ′ πj ′ Φ(tj ′ (z)) contains πj , and it is linear in πj . The explicit partial of Fθ at fixed z therefore reads off immediately: ∂Fθ = Φ tj (z) . (31) ∂πj z fixed Substituting (31) into the master formula, Φ tj (z) ∂z . (32) = − ∂πj fθ (z) The minus sign is the natural one. Φ(tj (z)) ∈ (0, 1) is the contribution of component j to the total probability mass at z, so increasing πj by a small amount pushes more probability to the left of z. To keep Fθ (z) = u unchanged we have to compensate by sliding z leftwards, hence the negative sign. The magnitude of the response is set by how peaked the overall density is near z (the 1/fθ (z) factor): in flat regions a small change in mass requires a large displacement of z, in peaked regions a tiny one. (b) ξ = µj (gradient with respect to a mean). From (26), only the j-th term depends on µj , and ∂tj (z)/∂µj = −1/σj (with z held fixed, because tj depends on µj only through the numerator z − µj ). The chain rule, together with Φ′ = φ, gives ∂tj (z) ∂Fθ πj = πj φ tj (z) = − φ tj (z) . (33) ∂µj z fixed ∂µj σj Plugging (33) into (30) cancels the minus sign in the master formula and yields πj φ tj (z) /σj ∂z . (34) = ∂µj fθ (z) Intuitively, shifting the j-th component to the right (increasing µj ) drags its probability mass to the right, so the cumulative mass below any fixed z decreases. To restore Fθ (z) = u we have to move z rightwards as well, matching the positive sign of the gradient. The size of the response is governed by how much component j contributes at z, namely πj φ(tj (z))/σj , divided by the total mixture density at z. Components that are far from z exert almost no influence, because φ(tj (z)) is exponentially small. (c) ξ = σj2 (gradient with respect to a variance). The slightly more delicate piece is ∂tj (z)/∂σj2 . Differentiating tj (z) = (z − µj ) σj−1 at fixed z and using the chain rule for σj−1 as a function of σj2 , dσj−1 z − µj ∂tj (z) tj (z) = (z − µj ) = (z − µj ) − 12 σj−3 = − = − , (35) 2 2 ∂σj dσj 2 σj3 2 σj2 where in the last step we used (z − µj )/σj = tj (z) to fold the prefactor back into tj and keep the formula compact. Hence ∂tj (z) ∂Fθ = πj φ tj (z) ∂σj2 z fixed ∂σj2 πj φ tj (z) tj (z) = − . (36) 2 σj2 Substituting (36) into the master formula, the minus sign once more cancels and we get πj φ tj (z) tj (z) πj φ tj (z) (z − µj )/σj ∂z = = . (37) ∂σj2 2 σj2 fθ (z) 2 σj2 fθ (z) 17

Summary. We started from a single defining identity, Fθ (z) = u, and differentiated it three times to read off how z responds to each mixture parameter. Each gradient ended up as a small, closed-form expression involving only the mixture PDF fθ (z), the standard-normal CDF Φ, the standard-normal PDF φ, and the standardized residual tj (z), all of which are computed once during the forward pass. The cost of the backward pass through the inverse CDF is therefore the cost of one forward evaluation, independent of how many Newton iterations the forward pass took. This is the key practical payoff of implicit differentiation in our setting.

Experimental Details

This appendix details the data-splitting protocol, training configuration, and hyperparameter selection used in our experiments. We use the four benchmark datasets (USHCN, PhysioNet’12, MIMIC-III, MIMIC-IV) following the preprocessing protocols described by Yalavarthi et al. [2025b], including the binning intervals, observation horizons, and channel selections reported there. We refer the reader to that work for the full dataset characteristics; here we focus on the protocol that governs how we train, evaluate, our model and its baselines. D.1

Dataset Statistics

For convenience we summarize the four datasets in Table 5, including the per-instance query-count statistics (minimum, average, maximum). USHCN [Menne et al., 2016] is a climatological dataset of daily measurements from 1,100 U.S. weather stations over four years (1996–2000), tracking five variables (snow precipitation, rain precipitation, snow depth, minimum temperature, maximum temperature). Following De Brouwer et al. [2019], we transform USHCN into an IMTS by randomly dropping 95% of the recorded observations. PhysioNet’12 [Silva et al., 2012] contains ICU records of 12,000 patients monitored over a 48-hour period across 37 vital signs, binned into hourly intervals [Che et al., 2018]. MIMIC-III [Johnson et al., 2016] provides ICU stays from Beth Israel Deaconess Medical Center, with 96 variables aggregated into 30-minute bins [De Brouwer et al., 2019]. MIMICIV [Johnson et al., 2023] records ICU stays from a tertiary academic medical center in Boston with 102 variables binned at the 1-minute resolution [Biloš et al., 2021]; this fine-grained binning yields the largest number of observations and queries among the four datasets. Table 5: Summary of dataset characteristics and preprocessing. The last three columns report the minimum, average, and maximum number of forecasting query points N per instance. Query points N Dataset USHCN PhysioNet’12 MIMIC-III MIMIC-IV D.2

Samples

Channels

Duration

Binning

min

N̄

max

1,100 (Stations) 12,000 (Patients) 21,000 (Patients) 18,000 (Patients)

5 37 96 102

4 Years 48 Hours 48 Hours 48 Hours

1 Week 1 Hour 30 Mins 1 Min

3 1 1 1

3.3 19.8 10.3 7.8

6 53 85 79

Data Splits

We adopt the splitting protocol of Yalavarthi et al. [2025b] without modification, both to ensure a fair head-to-head comparison with previously reported numbers and to remove any degree of freedom that could be tuned in our favor. Each dataset is partitioned into training, validation, and test sets in a 70:20:10 ratio. We repeat this partitioning under five distinct random seeds. We use the same seeds as experiments done in previous work [Yalavarthi et al., 2025b, Klötergens, 2026]. The test seed is tied to the fold index (e.g., Split 1 uses Seed 1), so the same five test sets are used by every model in our experiments. All reported numbers are means and standard deviations across these five folds. Interpreting the standard deviation. The five folds differ in which forecasting windows they expose to the model, and some folds contain inherently harder forecasting horizons than others (e.g., a higher proportion of unobserved channels at prediction time, or query points that fall further from the conditioning window). Consequently, the per-fold scores of every model on a given dataset shift up or down together as the fold becomes harder or easier, which inflates the across-fold standard deviation 18

relative to the much smaller within-fold gap between models. The relevant comparison is therefore the gap between two models on the same fold, aggregated across folds, rather than the absolute magnitude of either model’s standard deviation. As a rule of thumb, when a model outperforms a baseline by approximately one reported standard deviation in this protocol, the improvement is consistent across folds and statistically meaningful. D.3

Hyperparameter Selection

For each model and dataset, we draw 10 hyperparameter configurations from the search space below using random search, train each configuration on fold 0, and select the configuration with the best validation njNLL. The selected configuration is then retrained from scratch on all five folds; the test numbers we report come from this retrained configuration. We use the same procedure and the same search-space sizes for all baselines, taking their search spaces from the corresponding original publications. This keeps the per-model tuning budget identical across the comparison. Tables 6 to 8 list the search spaces for MargFlow, CoPFITi, and our re-implemented TACTiS-2 baseline (see Section D.8). CoPFITi and TACTiS-2 share the same encoder, so the encoder-side hyperparameters (copula_n_heads, copula_hidden_dim) are drawn from identical grids; the remaining hyperparameters are specific to each copula module. Table 6: Hyperparameter search space for MargFlow, the marginal model used as the first stage of CoPFITi and TACTiS-2. Hyperparameter

Description

Values

marg_weight_decay marg_n_heads marg_hidden_dim flow_mlp_layers flow_mlp_dim flow_layers flow_hid_dim

Weight decay of the marginal encoder Attention heads of the marginal encoder Hidden dim of the marginal encoder MLP layers of the conditioner network MLP hidden dim of the conditioner network Number of normalizing-flow layers Hidden dim of each normalizing-flow layer

{10−4 , 10−3 } {1, 2, 4} {32, 64, 128} {1, 2, 3} {32, 64, 128} {1, 2, 3} {10, 20}

Table 7: Hyperparameter search space for CoPFITi. Hyperparameter

Description

Values

copula_components copula_n_heads copula_hidden_dim corr_net_hidden_dim

Number of mixture components Attention heads of the encoder Hidden dim of the encoder Hidden dim of the MLPs producing µ, σ, and the Gram-matrix vectors

{1, 3, 5, 7, 10} {1, 2, 4} {32, 64, 128} {16, 32, 64, 128}

Table 8: Hyperparameter search space for our re-implemented TACTiS-2 baseline. The copula-side ranges match those reported by Ashok et al. [2023]; the encoder-side ranges (copula_hidden_dim, copula_n_heads) match CoPFITi. Hyperparameter

Description

Values

copula_attn_heads copula_attn_layers copula_attn_dim copula_mlp_layers copula_mlp_dim copula_resolution copula_dropout copula_hidden_dim copula_n_heads

Attention heads of the attentional copula Attention layers of the attentional copula Attention dim of the attentional copula MLP layers of the attentional copula MLP hidden dim of the attentional copula CDF discretization resolution Dropout rate of the attentional copula Hidden dim of the encoder Attention heads of the encoder

{2, 4, 8} {1, 2, 3, 4} {8, 16, 32} {1, 2, 3} {32, 64, 128} {8, 16, 32} {0.0, 0.1, 0.2} {32, 64, 128} {1, 2, 4}

D.4

Training Setup

All models are trained with the AdamW optimizer [Kingma and Ba, 2017, Loshchilov and Hutter, 2019] using an initial learning rate of 10−3 , weight decay of 10−3 , and a batch size of 64. We schedule the learning rate with reduce-on-plateau (factor 0.5, patience of 5 validation epochs without improvement) and train for at most 2000 epochs with early stopping (patience of 30 epochs on the validation njNLL). All experiments use float32 precision. For the MIMIC-IV dataset we reduce the batch size to 32 when memory pressure requires it; we apply this reduction uniformly across all models on that dataset. We train on a single NVIDIA A100 GPU per fold. D.5

Evaluation Metrics

We evaluate every model on the held-out test split of each fold using the normalized joint negative log-likelihood (njNLL) and the marginal negative log-likelihood (mNLL). The njNLL is defined in (15); the mNLL evaluates the predictive density at each query coordinate independently and averages over query points, mNLL(y | Q, X ) =

N 1 X qry − log p yn (tqry n , cn ), X . N n=1

(38)

We report the mean and standard deviation of both metrics across the five folds. njNLL is the primary metric (it is also the training objective) and mNLL serves as a diagnostic for marginalization consistency: a model that is internally consistent should not exhibit a large gap between its joint and marginal scores. D.6

The CoPFITi-Mix-GC Ablation

Recall from Section 2.1 that a GM-C is induced by a single multi-modal base density fGMM (z) = PK N j=1 πj N (z; µj , Σj ) on R via (3) using the GMM’s own marginal CDFs. A Mix-GC instead is a convex combination of K independent Gaussian copulas, c(u1 , . . . , uN ) =

K X

πj cRj (u1 , . . . , uN ),

(39)

j=1

where cRj is the Gaussian copula induced by a correlation matrix Rj . Each Gaussian copula component is centered at the origin in latent space and uses standard-normal marginals, so the per-component means µj and standard deviations σj,n that appear in the GM-C construction are no longer present, only the off-diagonal correlation structure remains. Architectural changes. The CoPFITi-Mix-GC ablation reuses every component of CoPFITi described in Section 4 unchanged, except for the parts that produce the per-component means and standard deviations. Concretely: • The encoder, the global summary h̄, and the per-query embeddings en are identical to CoPFITi. • The mixture-weight head π = softmax(MLPπ (h̄)) is identical to CoPFITi. • The mean head MLPµ ((11)) and the standard-deviation head MLPσ ((12)) are removed: each Gaussian copula component is centered at zero with unit variance by definition, so µj,n and σj,n are not free parameters. • The correlation network that produces the per-component feature matrices Uj ∈ RN ×H is unchanged. We form the Gram matrices Gj = Uj Uj⊤ + I and normalize them to correlation −1/2

−1/2

matrices Rj = DGj Gj DGj exactly as in CoPFITi, but stop there: Rj is now the full per-component covariance, not just its correlation structure. The hyperparameter search space therefore matches that of CoPFITi (Table 7) with the corr_net_hidden_dim entry reinterpreted as the hidden dimension of the only remaining MLP. 20

D.7

Univariate Wasserstein-1 Distance

The Wasserstein distances reported in Table 4 are per-dimension univariate Wasserstein-1 distances between K = 1000 samples drawn from the model and the corresponding ground-truth targets, averaged over the forecasting dimensions of each instance and then over the test set, and finally averaged over the five splits. We compare two model setups: MargFlow alone, and the full MargFlow+copula setup, either CoPFITi or TACTiS-2. The Control row quantifies the irreducible sampling error: it is computed by drawing two independent sets of K = 1000 samples from MargFlow under different random seeds and applying the same estimator to those two sample sets. For two probability measures P and Q on R with cumulative distribution functions FP and FQ , the Wasserstein-1 distance admits the closed form Z W1 (P, Q) = FP (t) − FQ (t) dt. (40) K

R K

Given two equal-sized samples {ai }i=1 and {bi }i=1 drawn from P and Q, the empirical Wasserstein1 distance reduces to the average absolute gap between order statistics, K X c1 {ai }, {bi } = 1 W a(i) − b(i) , K i=1

(41)

where a(1) ≤ · · · ≤ a(K) and b(1) ≤ · · · ≤ b(K) denote the sorted samples. For each forecasting model dimension d of an instance, we apply this estimator directly to the K = 1000 model samples {yd,i } GT and the corresponding ground-truth values {yd,i }, and average the resulting per-dimension distances across dimensions, instances, and splits. D.8

Adapting TACTiS-2 to sparse IMTS

The natural copula baseline for CoPFITi is the attentional copula introduced in TACTiS [Drouin et al., 2022] and reused in TACTiS-2 [Ashok et al., 2023]. The two methods share exactly the same attentional copula module; what changes between them is the training scheme. In the original TACTiS, a single shared encoder feeds both the marginals and the copula, and all components are trained jointly under the joint negative log-likelihood. In TACTiS-2, the marginals and the copula are produced by two separate encoders and trained in two stages: the marginals are trained first, then frozen, and the copula is trained on top of the frozen marginals. Ashok et al. [2023] report that this two-stage protocol consistently and substantially outperforms the joint TACTiS training across all of their tasks. The two-stage protocol is also exactly the protocol used by CoPFITi, so comparing against the TACTiS-2 form of the attentional copula gives the most direct copula-vs-copula comparison. One ingredient of the original TACTiS-2 does not transfer to our setting unchanged: its encoder. The TACTiS-2 encoder was designed for moderately irregular but still relatively dense multivariate time series, in which every channel is observed at most timestamps. The IMTS benchmarks we use are far sparser, with channels that may be observed only a handful of times over the whole observation window, and the original TACTiS-2 encoder underperforms badly in this regime. To isolate the contribution of the copula module itself, we therefore re-implement the attentional copula on top of the same encoder that CoPFITi uses (described in Section 3). Concretely, our TACTiS-2 baseline keeps the original attentional copula and the original two-stage training protocol, and only swaps the encoder for one that is suited to sparse IMTS. This is also the configuration in which it performs best on our benchmarks, so we use it everywhere we refer to TACTiS-2 in the main paper.

Sensitivity Analysis

We study the sensitivity of CoPFITi to the two hyperparameters that are specific to the GM-C copula module: the number of mixture components K and the per-component correlation-network hidden dimension H (i.e. the width of the MLPs that emit µj , σj , and the rows of the Gram-matrix factors Uj ∈ RN ×H , see Table 7). For each sweep, we vary the analyzed hyperparameter while fixing all remaining hyperparameters to the values selected by the random search described in Section D.3. Each configuration is trained from scratch on all five splits, and we report the mean test njNLL across the five splits. Number of mixture components K. Figure 3 sweeps K ∈ {1, 3, 5, 7, 10} on all four datasets. The case K = 1 collapses the GM-C to a single non-zero-mean Gaussian copula and is clearly inferior on every dataset. Moving from K = 1 to K = 3 produces a large drop in njNLL on all four benchmarks, and the curve essentially flattens beyond K ≥ 3. This pattern is consistent with the expressivity argument in Section 2.1: a single Gaussian copula cannot represent multimodal dependence. USHCN

Physionet

Test NJNLL

−3.9

MIMIC-III

−0.60

−0.6

−0.65

−4.0

−0.75

−2.0 −2.1

−0.8

−0.80

−4.2

−1.9

−0.7

−0.70

−4.1

MIMIC-IV −1.8

−2.2

−0.9 1

−2.3

Figure 3: Sensitivity of CoPFITi to the number of mixture components K. Test njNLL (lower is better) averaged over the five splits, with all remaining hyperparameters fixed to the values selected by the random search of Section D.3. Correlation-network hidden dimension H. Figure 4 sweeps H ∈ {16, 32, 64, 128}. Across all four datasets the test njNLL is essentially flat with respect to H. The only visible deviation is a small uptick on USHCN at H = 128, which we attribute to mild overfitting on the smallest dataset. It is worth highlighting that H = 16 is sufficient on MIMIC-III and MIMIC-IV even though many test instances have N > H in those datasets (the maximum query count N is 85 for MIMIC-III and 79 for MIMIC-IV; see Table 5). In this regime, each per-component covariance Σj = Uj Uj⊤ + I with Uj ∈ RN ×H has off-diagonal block of rank at most H, which is the limitation we flag in the main text: CoPFITi cannot represent full-rank or near-perfect dependence when H ≪ N . The flat sensitivity curves indicate that this rank ceiling is not an issue in practice, and may even be beneficial, since in reality not all variables depend on each other. USHCN

Physionet

Test NJNLL

−4.00 −4.05 −4.10 −4.15 −4.20 −4.25

MIMIC-III

−0.65

−0.75

−0.70

−0.80

−0.75

−0.85

−0.80

−0.90

MIMIC-IV

−2.15 −2.20 −2.25 −2.30

16 32

128

16 32

128

16 32

128

16 32

128

Figure 4: Sensitivity of CoPFITi to the correlation-network hidden dimension H. Test njNLL (lower is better) averaged over the five splits, with all remaining hyperparameters fixed to the values selected by the random search of Section D.3.

Computational Complexity

The dominant cost in evaluating the CoPFITi log-density at a single query of size N is the percomponent evaluation of the multivariate Gaussian density N (z; µj , Σj ). For each component j ∈ {1, . . . , K}, this requires solving a linear system in Σj and computing log det Σj , both of which are typically obtained from a Cholesky factorization of Σj at a cost of O(N 3 ). A naive implementation therefore scales as O(K N 3 ) per query. This worst-case cost is, however, never realized for CoPFITi. By construction (see Section 4), each covariance is built as −1/2

Σj = Dj Rj Dj ,

Rj ∝ DGj

−1/2

(Uj Uj⊤ + I) DGj ,

(42)

where Uj ∈ RN ×H and Dj , DGj are diagonal. The inner matrix Uj Uj⊤ + I is a rank-H update of the identity, so the Sherman-Morrison-Woodbury identity and the matrix determinant lemma yield −1

(Uj Uj⊤ + I)

−1

= I − Uj (I + Uj⊤ Uj )

Uj⊤ ,

(43)

log det(Uj Uj⊤ + I) = log det(I + Uj⊤ Uj ),

(44)

both of which reduce to operations on the H × H matrix I + Uj⊤ Uj . Combined with the diagonal rescalings, evaluating log N (z; µj , Σj ) for one component costs O(N H 2 + H 3 ), and the full logdensity scales as O K (N H 2 + H 3 ) rather than O(KN 3 ). If we operate in a regime, where H ≪ N , the cubic dependence on N is turned into a linear one. The remaining components of the forward pass are negligible by comparison: the encoder is O(M D2 + N D2 ) in the history and query lengths, the parameter MLPs are O((N + 1)D2 ), the Newton inversion of the GMM marginal CDFs is O(N K) per Newton step, and the marginal flow contributes O(N LM ) for L DSF blocks of width M . Empirical training time. Table 9 reports the average wall-clock training time per epoch for the three copula models considered in our experiments, measured on a single NVIDIA A100 GPU and averaged across the five splits of each dataset. Unlike the Mix-GC ablation, evaluating the −1 GM-C requires the numerical marginal iCDF FGMM,n in the forward pass via Newton iteration (Lemma 1), which adds a small overhead and is one of the reasons why CoPFITi is slightly slower than CoPFITi-Mix-GC. The gap is, however, modest, and CoPFITi remains faster than TACTiS-2 across all four datasets. Table 9: Average wall-clock training time per epoch (in seconds) on a single NVIDIA A100 GPU, averaged over 5 splits. Model

USHCN

Physionet

MIMIC-III

MIMIC-IV

TACTiS-2 CoPFITi-Mix-GC CoPFITi (GM-C)

0.28 0.16 0.22

6.1 1.9 2.2

4.2 3.0 4.0

11.3 2.4 3.5

Learning Curves

To complement the marginal-quality analysis in Section 6, we report the validation learning curves of the Joint-Ablation (Joint-Abl) and the decoupled MargFlow + CoPFITi setup for MIMIC-III in Figure 5. Joint-Abl shares the architecture of CoPFITi but trains all weights jointly under njNLL, whereas the decoupled setup first fits MargFlow in isolation and then trains the copula on top of the frozen marginals. All models were trained until the validation njNLL did stopped improving for 30 epochs. Learning curves: MargFlow+CopFITi vs. Joint-Abl (MIMIC-III) (a) Marginal NLL

(b) Normalized Joint NLL 1.0

1.0

0.5

0.6

NJNLL

MNLL

0.8

0.4 0.2 0.0

0.0 −0.5 −1.0

−0.2 0

250

500

750

1000

1250

1500

1750

Epoch MargFlow train MargFlow val

250

500

750

1000

1250

1500

1750

Epoch CopFITi train CopFITi val

Joint-Abl train Joint-Abl val

MargFlow → CopFITi

Figure 5: Validation learning curves of the Joint-Ablation (Joint-Abl) and the decoupled MargFlow + CoPFITi setup on MIMIC-III. Training the marginals in isolation reaches a lower marginal NLL than training all weights jointly under njNLL.

On TACTiS-2’s Marginalization Consistency

The Wasserstein gap of Table 4 probes only the univariate face of marginalization consistency: whether each per-dimension marginal of the model’s samples matches its MargFlow marginal. We single this slice out in the main paper because uniform univariate marginals are the defining property of a copula, and any deviation already violates the Sklar decomposition the model is built on. Marginalization consistency in its full form is strictly stronger: it asks that the model produce the same marginal on any subset of dimensions, regardless of which other dimensions are queried alongside. CoPFITi satisfies this stronger condition by construction (Section A), whereas TACTiS-2’s attentional copula guarantees neither the univariate nor the multivariate version, since every query subset induces a different autoregressive factorization and attention conditioning. We isolate the multivariate failure on a controlled 3-d toy distribution below: TACTiS-2’s bivariate marginal queried directly disagrees visibly with the same bivariate marginal obtained by marginalizing its 3-d joint, whereas CoPFITi’s two routes coincide up to sampling noise. We complement the validity-gap measurements of Table 4 with a small but diagnostic toy experiment that isolates the source of TACTiS-2’s marginalization gap from any confounder in the IMTS setting. Setup. We construct a 3-dimensional toy distribution with two symmetric Gaussian clusters indexed by a quasi-discrete cluster variable y2 . Cluster A is centered at (y1 , y2 , y3 ) ≈ (−2.5, −2, +2.5), and cluster B at (y1 , y2 , y3 ) ≈ (+2.5, +2, −2.5). Knowing y2 localizes y1 , so the true 2-d marginal p(y1 , y2 ) is unimodal per cluster, while the 2-d marginal p(y1 , y3 ) is a strongly bimodal anti-diagonal mixture. We train both models on the full 3-d distribution. Both models share the same pre-trained 1-d marginals p(yi ), so the only component that differs between them is the copula. Any disagreement in the predicted joints is therefore attributable to the copula alone. For each model, we obtain the 2-d marginals p(y1 , y2 ) by sampling in two different ways: 1. Direct: query only the dimensions (y1 , y2 ), with the third dimension absent from the query. 24

(a) Ground truth (3-D) (b) CopFITi; dir.

10.0

0.0

−5.0

−4

−10

−7.5

0 −2

−5

−2.5

10 −10

−10.0

5 −5

0 5

10 −10

−5

WD = 0.042

WD = 0.926

(d) GT (2-D)

(e) CopFITi; marg.

Density

2.5

5.0

7.5

(f ) TACTiS-2; marg.

4 4

2 0

0 2

−2

−4

−4 −4

−2

−4

−2

−4

−2

Figure 6: The Attentional Copula fails to be marginalization consistent. (a) shows the ground truth 3-d distribution. Panels (b) and (e) show CoPFITi’s direct and marginalized views of p(y1 , y2 ), which match by construction. Panels (c) and (f) show TACTiS-2’s direct and marginalized views, which disagree visibly; the arrows report the corresponding WD. 2. Marginalized: query the full 3-d distribution (y1 , y2 , y3 ). A marginalization consistent copula must produce the same distribution under both routes. We quantify the gap with the Wasserstein distance (WD) between the two sample sets, indicated by the arrows between rows 0 and 1 of Figure 6. Why TACTiS-2 fails. The attentional copula factorizes the joint as an autoregressive product of conditional CDFs over the queried dimensions, with attention conditioned on the specific index set S being modeled, Y pθ yi y<i , S . (45) pθ (yS ) = i∈S

Both the autoregressive ordering and the attention weights depend on S, so each query set instantiates a distinct factorization of the same network. During training on the full 3-d distribution the model only ever observes S = {1, 2, 3}, so the only factors that receive a gradient signal are those induced by that query, namely pθ (y3 ), pθ (y2 | y3 ), and pθ (y1 | y2 , y3 ). Evaluating the bivariate marginal under the direct route asks the same network for an entirely different set of factors, pθ (y2 ) and pθ (y1 | y2 ), neither of which is ever instantiated during training and neither of which is constrained by any loss. The attentional copula is therefore free to map this out-of-distribution conditioning to an arbitrary distribution, and no architectural mechanism ties its output back to the marginal of the trained 3-d joint. Empirically, the marginalized view inherits the bimodal structure that the copula had to learn for y3 under S = {1, 2, 3}, smearing the otherwise per-cluster unimodal p(y1 , y2 ), while the direct view is governed by these untrained factors; panels (c) and (f) consequently disagree and the WD is large.

Sampling Results

We complement the likelihood-based evaluation in the main paper with three sample-based metrics: Mean Squared Error (MSE) for point forecasts (Table 10), Energy Score (ES) for multivariate distributional accuracy (Table 11), and Continuous Ranked Probability Score (CRPS) for marginal distributional accuracy (Table 12). For each test instance we draw S = 1000 samples ŷ(s) ∈ RN from the predictive distribution conditioned on (Q, X ), and compare them against the ground-truth target y ∈ RN . All three metrics are reported as means across the test set and across the five folds. PS (s) MSE. We use the sample mean ŷ¯n = S1 s=1 ŷn as the point forecast and report the percoordinate squared error averaged over query points, N

MSE(ŷ, y) =

1 X ¯ 2 (ŷn − yn ) . N n=1

(46)

Energy Score. The Energy Score [Gneiting and Raftery, 2007] is a strictly proper scoring rule for multivariate distributions. For a predictive distribution F̂ and target y ∈ RN it is defined as ES(F̂ , y) = EY∼F̂ ∥Y − y∥2 − 21 EY,Y′ ∼F̂ ∥Y − Y′ ∥2 ,

(47)

where Y, Y′ are independent draws from F̂ . Given S samples, we use the standard unbiased estimator S X X ′ 1 c = 1 ES ∥ŷ(s) − y∥2 − ∥ŷ(s) − ŷ(s ) ∥2 . (48) S s=1 2 S (S − 1) ′ s̸=s

CRPS. The Continuous Ranked Probability Score [Gneiting and Raftery, 2007] is the univariate counterpart of the Energy Score. For a one-dimensional predictive CDF F̂n and target yn ∈ R, Z 2 CRPS(F̂n , yn ) = F̂n (z) − ⊮{z ≥ yn } dz = E |Yn − yn | − 21 E |Yn − Yn′ |, (49) R

with Yn , Yn′ ∼ F̂n independent. We report the average across query points, estimated from the (s) per-coordinate samples ŷn as   N S X X X ′ 1 1 1 \ =  |ŷ (s) − yn | − |ŷn(s) − ŷn(s ) | . (50) CRPS N n=1 S s=1 n 2 S (S − 1) ′ s̸=s

Lower values are better for all three metrics. Table 10: Comparing models w.r.t. MSE (Mean Squared Error) across datasets. Lower values indicate better performance. We present the mean and standard deviation of 5 runs. The best model is marked in bold. Model

USHCN

Physionet

MIMIC-III

MIMIC-IV

GRU-ODE NeuralFlows CRU ProFITi GPR MOSES CircuITS

0.410 ± 0.106 0.424 ± 0.110 0.290 ± 0.060 0.308 ± 0.061 0.597 ± 0.110 0.411 ± 0.099 0.306 ± 0.030

0.329 ± 0.004 0.331 ± 0.006 0.475 ± 0.015 0.305 ± 0.007 0.575 ± 0.059 0.307 ± 0.006 0.292 ± 0.000

0.479 ± 0.044 0.479 ± 0.045 0.725 ± 0.037 0.548 ± 0.063 0.862 ± 0.016 0.517 ± 0.057 0.475 ± 0.050

0.365 ± 0.012 0.374 ± 0.017 OOM 0.389 ± 0.015 0.609 ± 0.014 0.342 ± 0.028 0.285 ± 0.001

MargFlow

0.389 ± 0.154

0.310 ± 0.016

0.480 ± 0.047

0.278 ± 0.003

Table 11: Comparing models w.r.t. Energy Score (multivariate) across datasets. Lower values indicate better performance. We present the mean and standard deviation of 5 runs. The best model is marked in bold. Model

USHCN

Physionet

MIMIC-III

MIMIC-IV

NeuralFlows ProFITi MOSES CircuITS TACTiS-2

0.661 ± 0.059 0.452 ± 0.044 0.552 ± 0.044 0.468 ± 0.020 0.477 ± 0.052

1.691 ± 0.001 0.879 ± 0.303 1.599 ± 0.013 1.613 ± 0.000 1.613 ± 0.018

1.381 ± 0.033 1.606 ± 0.168 1.353 ± 0.033 1.489 ± 0.040 1.305 ± 0.031

0.982 ± 0.009 0.808 ± 0.003 0.906 ± 0.029 0.927 ± 0.050 0.825 ± 0.007

MargFlow CoPFITi

0.483 ± 0.054 0.482 ± 0.054

1.606 ± 0.012 1.593 ± 0.010

1.349 ± 0.033 1.340 ± 0.033

0.826 ± 0.005 0.826 ± 0.005

Table 12: Comparing models w.r.t. CRPS (Continuous Ranked Probability Score) on marginals across datasets. Lower values indicate better performance. We present the mean and standard deviation of 5 runs. The best model is marked in bold. Model

USHCN

Physionet

MIMIC-III

MIMIC-IV

NeuralFlows ProFITi MOSES CircuITS

0.306 ± 0.028 0.182 ± 0.007 0.220 ± 0.019 0.182 ± 0.010

0.277 ± 0.003 0.271 ± 0.003 0.260 ± 0.002 0.252 ± 0.001

0.308 ± 0.004 0.319 ± 0.003 0.296 ± 0.005 0.286 ± 0.010

0.281 ± 0.004 0.279 ± 0.012 0.245 ± 0.010 0.221 ± 0.001

MargFlow

0.187 ± 0.022

0.251 ± 0.001

0.289 ± 0.005

0.217 ± 0.001

Broader Impacts

CoPFITi is a methodological contribution to probabilistic forecasting of irregular multivariate time series, evaluated on standard public benchmarks. We do not propose a deployed system, and the released artifacts are model code and training scripts rather than pre-trained generative models, scraped corpora, or systems intended for direct end-user use. Potential positive impacts. Three of the four benchmarks used in our evaluation (PhysioNet-2012, MIMIC-III, and MIMIC-IV) are clinical, and one (USHCN) is environmental. In such settings, well-calibrated probabilistic forecasts with reliable uncertainty estimates are arguably more useful than point predictions: they enable downstream decision-making to incorporate forecast uncertainty rather than treating model outputs as ground truth. CoPFITi’s marginalizationconsistency guarantee additionally ensures that predictions over different subsets of query points cannot contradict one another, which is a basic prerequisite for any model whose outputs may be inspected at multiple resolutions or by multiple downstream consumers. Potential negative impacts. The general risks of probabilistic forecasting apply: forecasts may be over-trusted, miscalibrated under distribution shift, or used in decision pipelines for which the original training distribution is no longer representative. This is particularly relevant in clinical settings, where unwarranted confidence in model output can have direct consequences for patients. We emphasize that CoPFITi is evaluated on retrospective benchmark data and is not a clinical decision-support tool, and that any deployment in a high-stakes domain would require dataset-specific validation, calibration analysis, and appropriate human oversight that go beyond the scope of this work. We are not aware of dual-use concerns or specific misuse scenarios beyond those generic to probabilistic time-series modeling.

NeurIPS Paper Checklist 1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope? Answer: [Yes] Justification: The two claims of the abstract, that CoPFITi is marginalization consistent by construction and that it establishes a new state of the art in joint IMTS density modeling, are formalized and supported in Section 4 (with proof in Section A) and in Section 6, respectively. Guidelines: • The answer [N/A] means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A [No] or [N/A] answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 2. Limitations Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: A dedicated Limitations paragraph at the end of Section 4 discusses the rankH structure of the per-component covariances and the constraint that the mixture weights π depend only on the observation history X , which is required by our marginalizationconsistency guarantee. Guidelines: • The answer [N/A] means that the paper has no limitation while the answer [No] means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate “Limitations” section in their paper. • The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren’t acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 28

3. Theory assumptions and proofs Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Justification: The marginalization-consistency claim is stated in Section 4 and proven in full in Section A. The closed-form gradients of the GMM inverse CDF (Lemma 1) are derived step-by-step in Section C. Guidelines: • The answer [N/A] means that the paper does not include theoretical results. • All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems. • The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 4. Experimental result reproducibility Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: We adopt the publicly available preprocessing, binning, and split protocol of prior work [Biloš et al., 2021, Yalavarthi et al., 2025a,b, Klötergens, 2026] on four standard public benchmarks. The architecture is fully specified in Sections 3 and 4, and hyperparameters and training details are listed in Section D. Guidelines: • The answer [N/A] means that the paper does not include experiments. • If the paper includes experiments, a [No] answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general, releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully. (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset). 29

(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 5. Open access to data and code Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: The four datasets are publicly available and we use the established preprocessing pipeline from prior work, cited in Section 6. An anonymized code release accompanies the submission and contains instructions to reproduce all reported results. Guidelines: • The answer [N/A] means that paper does not include experiments requiring code. • Please see the NeurIPS code and data submission guidelines (https://neurips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so [No] is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //neurips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 6. Experimental setting/details Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer) necessary to understand the results? Answer: [Yes] Justification: Data splits follow the established protocol cited in Section 6. Hyperparameters, optimizer settings, and training schedule for MargFlow and CoPFITi are reported in Section D, and the TACTiS-2 adaptation is described in Section D.8. Guidelines: • The answer [N/A] means that the paper does not include experiments. • The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 7. Experiment statistical significance Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] 30

Justification: Each model is trained and evaluated on five paired train/validation/test splits, and we report the mean and standard deviation across these splits in all result tables. Because the splits are paired across methods, we additionally report a Corrected Resampled t-test [Nadeau and Bengio, 1999] in Table 3, with the interpretation discussed in Section D. Guidelines: • The answer [N/A] means that the paper does not include experiments. • The authors should answer [Yes] if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g., negative error rates). • If error bars are reported in tables or plots, the authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 8. Experiments compute resources Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: All experiments are run on a single NVIDIA A100 GPU. Section F gives a complexity analysis and reports per-epoch wall-clock training times for the three copula models in Table 9. Guidelines: • The answer [N/A] means that the paper does not include experiments. • The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn’t make it into the paper). 9. Code of ethics Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: We have reviewed the NeurIPS Code of Ethics and our work conforms with it in all respects. We use only publicly available, anonymized benchmark datasets and do not collect any new data. Guidelines: • The answer [N/A] means that the authors have not reviewed the NeurIPS Code of Ethics. 31

• If the authors answer [No], they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: Section J discusses both potential positive impacts (better-calibrated probabilistic forecasts in clinical and environmental benchmarks, marginalization consistency as a prerequisite for trustworthy multi-resolution use) and potential negative impacts (over-trust, miscalibration under distribution shift, and the general risks of probabilistic forecasting in high-stakes domains). Guidelines: • The answer [N/A] means that there is no societal impact of the work performed. • If the authors answer [N/A] or [No], they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate Deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 11. Safeguards Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pre-trained language models, image generators, or scraped datasets)? Answer: [N/A] Justification: The paper does not release pre-trained generative models or scraped datasets. The released artifacts are model code and training scripts for a probabilistic forecasting model on existing public benchmarks, which we judge to pose no high risk of misuse. Guidelines: • The answer [N/A] means that the paper poses no such risks. • Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. 32

• We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: The four datasets (USHCN, PhysioNet-2012, MIMIC-III, MIMIC-IV) and all baseline methods are cited at their original sources in Section 6 and used in accordance with their respective terms of use. Guidelines: • The answer [N/A] means that the paper does not use existing assets. • The authors should cite the original paper that produced the code package or dataset. • The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset. • For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset’s creators. 13. New assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: The accompanying anonymized code release contains documentation, configuration files, and scripts to reproduce the reported results. Guidelines: • The answer [N/A] means that the paper does not release new assets. • Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 14. Crowdsourcing and research with human subjects Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? Answer: [N/A] Justification: The paper does not involve crowdsourcing or research with human subjects. We use only publicly available, previously released benchmark datasets. Guidelines: • The answer [N/A] means that the paper does not involve crowdsourcing nor research with human subjects. 33

• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 15. Institutional review board (IRB) approvals or equivalent for research with human subjects Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [N/A] Justification: The paper does not involve research with human subjects. The medical benchmarks used here were collected and de-identified by their original creators under their own ethical approvals. Guidelines: • The answer [N/A] means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 16. Declaration of LLM usage Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigor, or originality of the research, declaration is not required. Answer: [N/A] Justification: LLMs are not part of the core methodology of this paper. Guidelines: • The answer [N/A] means that the core method development in this research does not involve LLMs as any important, original, or non-standard components. • Please refer to our LLM policy in the NeurIPS handbook for what should or should not be described.

Record · ID 222603 · SHA-256 f7bee5a248b61d1a
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.