# SEMF: Supervised Expectation-Maximization Framework for Predicting Intervals

Ilia Azizi
HEC, University of Lausanne
Lausanne, Switzerland
<EMAIL_ADDRESS>
&Marc-Olivier Boldi
HEC, University of Lausanne
Lausanne, Switzerland
<EMAIL_ADDRESS>
&Valérie Chavez-Demoulin
HEC, University of Lausanne
Expertise Center for Climate Extremes (ECCE)
Lausanne, Switzerland
<EMAIL_ADDRESS>
Corresponding author.

###### Abstract

This work introduces the Supervised Expectation-Maximization Framework (SEMF),
a versatile and model-agnostic framework that generates prediction intervals
for datasets with complete or missing data. SEMF extends the Expectation-
Maximization (EM) algorithm, traditionally used in unsupervised learning, to a
supervised context, enabling it to extract latent representations for
uncertainty estimation. The framework demonstrates robustness through
extensive empirical evaluation across 11 tabular datasets, achieving—in some
cases—narrower normalized prediction intervals and higher coverage than
traditional quantile regression methods111Code:
https://github.com/Unco3892/semf-paper. Furthermore, SEMF integrates
seamlessly with existing machine learning algorithms, such as gradient-boosted
trees and neural networks, exemplifying its usefulness for real-world
applications. The experimental results highlight SEMF’s potential to advance
state-of-the-art techniques in uncertainty quantification.

##  1 Introduction

In the evolving field of machine learning (ML), the quest for models capable
of predicting outcomes while quantifying the uncertainty of their predictions
is critical. Estimating prediction uncertainty is particularly vital in high-
stakes domains such as healthcare [1], finance [2], and autonomous systems
[3], where prediction-based decisions can have significant consequences.
Traditional approaches have primarily focused on point estimates, providing
limited insight into the reliability of these predictions. This limitation
underscores the necessity for frameworks that generate precise point
predictions and robust prediction intervals. These intervals offer a range
within which the true outcome is expected to lie with a specified probability,
enabling a deeper understanding of prediction uncertainty. Consequently,
research has increasingly focused on methodologies that extend beyond point
estimation to include uncertainty quantification, thereby enhancing decision-
making in applications that rely on predictive modeling [4].

In this paper, we introduce the Supervised Expectation-Maximization Framework
(SEMF) based on the Expectation-Maximization (EM) algorithm [5]. Traditionally
recognized as a clustering technique, EM is used for supervised learning in
SEMF, allowing for both point estimates and prediction intervals using any ML
model (model-agnostic). SEMF generates representations for latent or missing
modalities, making it particularly relevant for incomplete multi-modal data.
This paper details the methodology behind the framework and proposes a
training algorithm based on Monte Carlo (MC) sampling, also used in
variational inference for Variational Auto-Encoders (VAEs) [6, 7]. Comparisons
with other models across several datasets showcase its contribution to the
existing uncertainty quantification and predictive modeling literature [8].

The remainder of this paper is organized as follows: Section 2 details the
theory and the methodology of SEMF. Section 3 reviews related works in latent
representation learning, uncertainty estimation, and handling of missing data.
Section 4describes the experimental setup, including datasets and evaluation
metrics. Section 5 discusses the results, demonstrating the efficacy of SEMF.
Lastly, Section 6 concludes the paper, and Section 7 outlines the limitations
and potential research directions.

##  2 Method

This section presents the founding principles of SEMF from its parameters,
training, and inference procedure with, at its core, the EM algorithm. This
algorithm, first introduced by [5], is an unsupervised method for handling
latent variables and incomplete data. Invented to maximize the model
likelihood, it builds a sequence of parameters that guarantee an increase in
the log-likelihood [9] by iterating between the Expectation (E) and the
Maximization (M) steps. In the E-step, one computes

$Q(p|p^{\prime})=\mathbb{E}_{Z\sim p^{\prime}(z|x)}\left[\log
p(x,Z)\right]=\int\log p(x,z)p^{\prime}(z|x)dz,$ (1)

where $p^{\prime}$ stands for the current estimates, $\log p(x,z)$ is the log-
likelihood of the complete observation $(x,z)$, and $z$ is a latent variable.
The M-step maximizes this $Q$-function:
$p^{\prime}\leftarrow\arg\max_{p}Q(p|p^{\prime})$ . The sequence is repeated
until convergence.

The EM algorithm has already been used for supervised learning tasks using
specific models [10, 11, 12, 13], and, to the best of our knowledge, mainly
focusing on point prediction through the use of Gaussian Mixture Models
(GMMs). Our work differs by modifying the EM algorithm and MC sampling to
generate prediction intervals with any ML model and, in theory, under any
distribution. This results in SEMF, a framework consisting of an encoder-
decoder architecture and a missing data simulator.

###  2.1 Problem Scenario

Let $x=(x_{1},x_{2},\ldots,x_{K})$ denote $K$ inputs and the output be $y$.
For simplicity, we limit $y$ to be numerical, although it could be categorical
without loss of generality. Component $x_{k}$ is a source: a modality, a
single or group of variables, or an unstructured input such as an image or
text. For clarity, we limit to $K=2$, where $x_{1}$ and $x_{2}$ are single
variables. We assume that only $x_{1}$ may contain missing values, either at
random or partially at random.

Let $p(y|x)$ be the density function of the outcome given the inputs (the fact
that $y$ is continuous can be easily relaxed). A founding assumption, in the
spirit of VAE, is that $p(y|x)$ decomposes into $p(y|x)=\int
p(y|z)p(z_{1}|x_{1})p(z_{2}|x_{2})dz_{1}dz_{2}$, where $z=(z_{1},z_{2})$ are
unobserved latent variables. We assume that $p(y|z,x)=p(y|z)$, that is, $z$
contains all the information of $x$ about $y$, and that
$p(z|x)=p(z_{1}|x_{1})p(z_{2}|x_{2})$, that is, there is one latent variable
per source. These are independent conditionally on their corresponding source.
Finally, if $x_{1}$ is missing, then $p(y|x_{2})=\int
p(y|z)p(z_{1}|x_{1})p(x_{1}|x_{2})p(z_{2}|x_{2})dx_{1}dz_{1}dz_{2}$. The
contribution to the log-likelihood of a complete observation $(y,z,x)$ is
$\log p(y,z|x)=\log p(y|z)+\log p(z|x)$. In the E-step, we compute

$\displaystyle\int\log p(y,z|x)p^{\prime}(z|y,x)dz$ $\displaystyle=$
$\displaystyle\int\log p(y|z)p^{\prime}(z|y,x)dz+\int\log
p(z|x)p^{\prime}(z|y,x)dz.$ (2)

where $p^{\prime}$ is our current estimate. Eq. 2 can be estimated by MC
sampling. Since sampling from $p^{\prime}(z|y,x)$ can be inefficient, we
rather rely on the decomposition
$p^{\prime}(z|y,x)=p^{\prime}(y|z)p^{\prime}(z|x)/p^{\prime}(y|x)$. Thus, we
sample $z_{r}$ from $p^{\prime}(z|x)$, $r=1,\ldots,R$, and, setting
$w_{r}=p^{\prime}(y|z_{r})/\sum_{t}p^{\prime}(y|z_{t})$, approximate the
right-hand side term of Eq. 2

$\displaystyle\int\log p(y,z|x)p^{\prime}(z|y)dz$ $\displaystyle\approx$
$\displaystyle\sum_{r=1}^{R}\left\\{\log p(y|z_{r})+\log
p(z_{r}|x)\right\\}w_{r}.$ (3)

If $x_{1}$ is missing, a similar development leads to

$\displaystyle\int\log p(y,z,x_{1}|x_{2})p^{\prime}(z,x_{1}|y)dz$
$\displaystyle\approx$ $\displaystyle\sum_{r=1}^{R}\left\\{\log
p(y|z_{r})+\log p(z_{r}|x)+\log p(x_{1,r}|x_{2})\right\\}w_{r},\quad$ (4)

where $x_{1,r}$ and $z_{r}$ are respectively sampled from
$p^{\prime}(x_{1}|x_{2})$ and $p^{\prime}(z|x_{1,r},x_{2})$.

###  2.2 Objective Function

Adapting Eq. 3 and Eq. 4 for the observed data
$\\{(y_{i},x_{i})\\}_{i=1}^{N}$, the overall loss function, $\mathcal{L}$, is

$\mathcal{L}(\phi,\theta,\xi)=-\sum_{i=1}^{N}\sum_{r=1}^{R}\left\\{\log
p_{\phi}(z_{i,r}|x_{i,r})+\log p_{\theta}(y_{i}|z_{i,r})+\mathds{1}_{\\{i\in
I_{m}\\}}\log p_{\xi}(x_{1,i,r}|x_{2,i})\right\\}w_{i,r},$ (5)

where $I_{m}$ is the set of those $i$’s such that $x_{1,i}$ is missing. The
models of $p(y|z)$, $p(z|x)$, and $p(x_{1}|x_{2})$ inherit parameters
$\theta$, $\phi$, and $\xi$, respectively. Also, $x_{1,i,r}$ is sampled from
$p_{\xi^{\prime}}(x_{1}|x_{2,i})$, if $x_{1,i}$ is missing, and $z_{1,i,r}$
and $z_{2,i,r}$ are sampled from $p_{\phi^{\prime}}(z_{1}|x_{1,i,r})$ and
$p_{\phi^{\prime}}(z_{2}|x_{2,i})$. Furthermore, for compactness of notation,
$x_{i,r}$ is $(x_{1,i},x_{2,i})$ if $x_{1,i}$ is observed, and
$(x_{1,i,r},x_{2,i})$ if $x_{1,i,r}$ is missing. Finally, the weights are

$w_{i,r}=\frac{p_{\theta^{\prime}}(y_{i}|z_{i,r})}{\sum_{t=1}^{R}p_{\theta^{\prime}}(y_{i}|z_{i,t})}.$
(6)

Eq. 5 shows that ${\cal L}$ is a sum of losses associated with the encoder
model, $p_{\phi}$, for each source, the decoder model, $p_{\theta}$, from the
latent variables to the output, and, if applicable, the model handling missing
data, $p_{\xi}$. At each M-step, ${\cal L}$ is minimized with respect to
$\theta$, $\phi$, and $\xi$. Then, $\theta^{\prime}$, $\phi^{\prime}$, and
$\xi^{\prime}$ are updated, as well as the weights and the sampling. Then the
process is iterated until convergence.

####  2.2.1 Example: $\mathcal{L}(\phi,\theta,\xi)$ under normality

Similar to [7], we develop further ${\cal L}$ under normality assumptions for
the encoder, $p_{\phi}$, and the decoder, $p_{\theta}$. These simple cases are
illustrative, though any other distributions could be adopted, including non-
continuous or non-numerical outcomes.

##### Encoder $p_{\phi}(z|x)$.

Let $m_{k}$ be the length of the latent variable $z_{k}$, $k=1,2$. We assume a
normal model for $Z_{k}$ given $X_{k}=x_{k}$,

$Z_{k}|X_{k}=x_{k}\sim{\cal
N}_{m_{k}}(g_{\phi_{k}}(x_{k}),\sigma_{k}^{2}J_{m_{k}}),$ (7)

where $J_{m_{k}}$ is the $m_{k}\times m_{k}$ identity matrix. In particular,

$\displaystyle\log p_{\phi}(z_{k}|x_{k})$ $\displaystyle=$
$\displaystyle-\frac{m_{k}}{2}\log
2\pi-\frac{1}{2}\log\sigma_{k}^{2}-\frac{1}{2\sigma_{k}^{2}}\sum_{j=1}^{m_{k}}\\{z_{k,j}-g_{\phi_{k},j}(x_{k})\\}^{2},\quad
k=1,2.$ (8)

The mean $g_{\phi_{k}}(x_{k})$ can be any model, such as a neural network,
with output of length $m_{k}$, $k=1,2$. The scale $\sigma_{k}$ can be fixed,
computed via the weighted residuals, or learned through a separate set of
models. It controls the amount of noise introduced in the latent dimension and
is pivotal in determining the prediction interval width for $p(y|z)$. In this
paper, $\sigma_{k}$ is fixed for simplicity.

##### Decoder $p_{\theta}(y|z)$.

We assume a normal model for $Y$ given $Z=z$,

$Y|Z=z\sim\mathcal{N}(f_{\theta}(z),\sigma^{2}).$ (9)

This results in a log-likelihood contribution,

$\log p_{\theta}(y|z)=-\frac{1}{2}\log
2\pi-\frac{1}{2}\log\sigma^{2}-\frac{1}{2\sigma^{2}}\\{y-f_{\theta}(z)\\}^{2}.$
(10)

Again, the mean $f_{\theta}(z)$ can be any model, such as a neural network.

##### Model for missing data $p_{\xi}(x_{1}|x_{2})$.

We use an empirical model for $X_{1}$ given $X_{2}=x_{2}$ where
$p_{\xi}(x_{1}|x_{2})$ put masses only on those non-missing $x_{1}$’s in the
training set. Let $I_{nm}$ be the set of indices such that $x_{1,j}$, $j\in
I_{nm}$, are all the non-missing $x_{1}$ in the training set. Additionally,
for a given $j\in I_{nm}$, $x_{1}[j]$ is the observed $x_{1}$ corresponding to
$j$. For a given $x_{2}$, $p_{\xi}(x_{2})$ is a vector of length $|I_{nm}|$,
the cardinality of $I_{nm}$, with components

$p_{\xi}(j|x_{2})=\frac{\exp\\{h_{\xi,j}(x_{2})\\}}{\sum_{t\in
I_{nm}}\exp\\{h_{\xi,t}(x_{2})\\}},\quad j\in I_{nm},$ (11)

where $h_{\xi}$ is a vector of length $|I_{nm}|$, typically a neural network
with input $x_{2}$ and output on $I_{nm}$. Now, the probability of
$X_{1}=x_{1}$ given $X_{2}=x_{2}$ is

$p_{\xi}(x_{1}|x_{2})=\sum_{j\in
I_{nm}}p_{\xi}(j|x_{2})\cdot\mathds{1}\\{x_{1}=x_{1}^{(nm)}[j]\\}.$ (12)

##### Summary.

Overall, the M-step is

$\phi_{k}^{*}=\arg\min_{\phi_{k}}\sum_{i,r}w_{i,r}\sum_{j=1}^{m_{k}}\\{z_{k,i,r,j}-g_{\phi_{k},j}(x_{k,i,r})\\}^{2},\quad
k=1,2,$ (13)
$\theta^{*}=\arg\min_{\theta}\sum_{i,r}w_{i,r}\\{y_{i}-f_{\theta}(z_{i,r})\\}^{2},$
(14)
$(\sigma^{*})^{2}=\frac{1}{N}\sum_{i,r}w_{i,r}\\{y_{i}-f_{\theta^{*}}(z_{i,r})\\}^{2},$
(15) $\xi^{*}=\arg\max_{\xi}\sum_{t\in D_{x}}w_{t}\log
p_{\xi}(j_{t}|x_{2,t}).$ (16)

When $x_{1,i}$ is missing, the sampling is enriched by simulated $x_{1,i,r}$.
Eq. 16 selects $j^{*}$ from $I_{nm}$ based on $p_{\xi^{\prime}}(j|x_{2})$
according to Eq. 12 and Eq. 11. The above $D_{x}$ is a subset of the data for
$\xi$ is only learned on the missing data part. After discussing the training
steps in Subsection 2.3, the details are given in Algorithm 1 (Appendix A). We
also note that the parameters above can be learned in parallel.

###  2.3 Training

For efficiency purpose, the training set, $\\{1,\ldots,N\\}$, is segmented
into batches $\\{b_{1},\ldots,b_{L}\\}$ on which the index $i$ runs (and thus
the denominator of Eq. 15 must be adapted accordingly). The process iterates
for each batch until the maximum number of steps is reached or an early
stopping criterion is satisfied. The framework requires tuning hyper-
parameters such as the number of MC samples $R$, the number of latent nodes
$m_{k}$, and the standard deviation $\sigma_{k}$ of $Z_{k}$. Monitoring the
point prediction on a hold-out validation is important to combat overfitting
and terminate the training early with a patience hyper-parameter. Moreover,
due to the generative nature of SEMF, the variation resulting from the initial
random seed is measured in Subsection 4.2. Additionally. the model-specific
hyper-parameters ($p_{\phi}$, $p_{\theta}$ and $p_{\xi}$) are also discussed
in the same Subsection.

###  2.4 Inference

The encoder-decoder structure of SEMF entails the simulations of $z_{r}$
during inference, as depicted in Figure 1. In theory, any inference can be
performed for $\hat{y}$, for instance the mean value
$\hat{y}=\frac{1}{R}\sum_{r=1}^{R}f_{\theta}(z_{r})$, where $z_{r}\sim
p_{\phi}(z|x)$ (see Algorithm 2 in Appendix A). For prediction intervals, a
double simulation scheme is used,

$z_{r}\sim p(z|x),\quad\hat{y}_{r,s}\sim p_{\theta}(y|z_{r}),\quad
r,s=1,\ldots,R.$ (17)

Prediction interval at a given level of certainty $\alpha$ follows as

$PI=\text{quantile}\left(\\{\hat{y}_{r,s}\\};\frac{\alpha}{2},1-\frac{\alpha}{2}\right).$
(18)

##### Remark.

We denote the $R$ for inference as $R_{\mathrm{infer}}$.

Figure 1: Inference procedure with the SEMF’s learnable parameters
$\phi_{k},\theta$ and $\xi$. Here, we illustrate the number of inputs $k$ as
$k=1,2$, assuming that $x_{1}$ may contain missing values

##  3 Related Work

###  3.1 Latent Representation Learning

Latent representation learning involves modeling hidden variables from
observed data for various ML tasks, most notably Auto-Encoders (AEs) and VAEs.
AEs are neural networks reconstructing inputs by learning an intermediate
latent representation. The leftmost part of AE is the encoder
$g_{\phi}(\cdot)$, which embeds the input $x$ into a latent variable $z$. It
is passed to the rightmost part of AE, the decoder $f_{\theta}(\cdot)$, for
the reconstruction $\hat{x}=f_{\theta}(z)$ of $x$ through $z$. The parameters
$\phi$ and $\theta$ are trained by minimizing the reconstruction loss. Unlike
AEs, which focus on reconstruction, VAEs use variational methods to fit
distributions of latent variables and of the output [7]. Given a sample $x$
with a latent variable of $z$, VAEs learn the marginal likelihood,
$p_{\theta}(x)$. Direct maximization of this likelihood is difficult [14].
Thus, alternative methods exist, such as maximizing the evidence lower bound
(ELBO), which provides guarantees on the log-likelihood [15].

For supervised and semi-supervised tasks, latent representation learning can
include task-specific predictions [16]. More specifically, models such as AEs
follow the classical encoder-decoder objective while training a predictor
$h_{\psi}(z)$ through an additional layer or model to estimate the output $y$.
This dual objective helps in learning more task-relevant embeddings [17, 18].
Semi-supervised VAEs are similar, with the distinction that they couple the
reconstruction loss of the unlabeled data with a variational approximation of
latent variables. This is effective even with sparse labels [19, 20].
Similarly, the EM algorithm adapts well to minimal supervision [21] and using
labeled and unlabeled data in semi-supervised settings for both single and
multiple modalities [12, 22].

###  3.2 Prediction Intervals

Crucial for estimating uncertainty, prediction intervals in regression are
often derived using methods such as Bayesian approaches [23, 24, 25], ensemble
techniques [26, 27], or quantile regression [28, 29]. Additionally, these
methods can be complemented with conformal prediction, a framework for
calibrating any point predictor to produce prediction intervals [30, 31],
making it highly relevant for enhancing reliability in applications requiring
rigorous uncertainty quantification.

A key component of quantile regression is the pinball loss function, which
effectively balances the residuals to capture the desired quantiles even for
non-parametric models [32], making it ideal for asymmetric distributions where
tail behavior is of critical importance [29]. This loss function is pivotal
not only for single model scenarios but also enhances ensemble methods by
refining their quantile estimate [33]. Conformal prediction further extends
the applicability of these intervals by providing a layer of calibration that
adjusts intervals obtained from any predictive model, ensuring they cover the
true value with a pre-specified probability [34].

###  3.3 Missing Data

Managing missing values is a pivotal aspect when dealing with real-world data.
Naive methods such as discarding instances or mean/median-imputation may be
infeasible or carry the risk of changing the data distribution [35, 36]
introducing bias if the missing values are not missing at random (MAR) or
completely at random (MCAR). The chosen technique for handling missing data
should adhere to the dataset’s characteristics and mechanisms behind the
missing data [37].

More advanced approaches, like the Iterative Imputer from scikit-learn [38],
an implementation of Multiple Imputation by Chained Equations (MICE) [39],
expand on the simple imputations by iteratively modeling each feature with
missing values as a function of other features [40, 41]. Due to its
complexity, MICE is best compared with alternatives such as K-means clustering
and artificial neural networks (ANNs). K-means assigns missing values based on
cluster centroids [42] as opposed to ANNs, which learn complex, non-linear
relationships between variables [43]. ANNs effectively predict or reconstruct
missing values and thus are particularly useful in datasets where
relationships between variables are intricate and not easily captured by
straightforward imputation methods. Accordingly, ANNs have demonstrated
superior performance over MICE and GMM (with EM for missing values) in
scenarios where a large proportion of the data is missing [44].

##  4 Experimental Setup

###  4.1 Datasets

We systematically curate a subset of datasets from the OpenML-CTR23 [45]
benchmark suite to evaluate and carry out our experiments. Initially
comprising 35 datasets, we apply an exclusion criteria to refine this
collection to 11 datasets. The details and overview are in Appendix B. We
remove duplicated rows from all the datasets and carry out the standardization
(scaling) of all predictors, including the outcome, which we transformed to
have zero means and unit variances. The features of these datasets are then
treated as separate inputs to SEMF. In the second stage of our experiments, we
artificially introduce 50% missing values in our datasets for any predictor
except for the first feature of a randomly chosen row, which emulates MAR or
MCAR data. In all our datasets, 70% of the data is used to train all models,
15% as a hold-out validation set to monitor SEMF’s performance, and 15% to
evaluate the models. To combat overfitting, baseline models that benefit from
early stopping are allocated another 15% from the training data. Lastly, it is
essential to note that all data in SEMF are processed batch-wise, without
employing mini-batch training, to ensure consistency and stability in the
training process.

###  4.2 Models

Our baseline consists of both point and quantile regression eXtreme Gradient
Boosting (XGBoost) [46], Extremely Randomized Trees (ET) [47], and neural
networks [48], depicted in Table 1. To ensure consistency in our experimental
setup, we align the families and hyper-parameters of $p_{\phi}$ and
$p_{\theta}$ with our baseline models. For example, in the case of XGBoost in
SEMF, we use $K$ XGBoost’s $g_{\phi_{k}}(x_{k})$, one for each input $x_{k}$,
$k=1,\ldots,K$, and one XGBoost for $f_{\theta}(z)$ with the same hyper-
parameters. We refer to the SEMF’s adoption of these models as MultiXGBs,
MultiETs, and MultiMLPs. When establishing prediction intervals, we set our
uncertainty tolerance at 5% for both the baseline and SEMF (Eq. 18). The
missing data simulator, $p_{\xi}$, is constructed using a shallow neural
network, which employs the $\mathrm{SELU}$ activation function [49]. It is
experimented with two distinct node counts: 50 and 100.

To constrain the breadth of our parameter exploration, the simulator for the
missing model adopts the optimal set of hyper-parameters identified from
analyses involving complete datasets. We target (larger) $\sigma_{k}$ values
that introduce more noise and produce better intervals than point predictions.
The optimal models are then trained and tested with five different seeds, and
the results are averaged. The point prediction, $\hat{y}$, uses the mean
inferred values. We then study the performance of SEMF against mean and median
imputation techniques, five nearest neighbors, and MICE from [38]. The
imputers are used within the point and interval baseline models explained in
the following subsection to form the missing value baseline. Appendix C
contains more details on the hyper-parameters for each SEMF model and dataset.

Table 1: SEMF models, baselines, and hyper-parameters.

SEMF | Point Prediction Baseline | Interval Prediction Baseline
---|---|---
MultiXGBs |  XGBoost [50]
Trees: 100, Maximum depth: 5, Early stopping steps: 10 |  Quantile XGBoost
Same as point prediction baseline, XGBoost
MultiETs |  Extremely Randomized Trees [38]
Trees: 100, Maximum depth: 10 |  Quantile Extremely Randomized Trees [51]
Same as point prediction baseline, Extremely Randomized Trees
MultiMLPs |  Deep Neural Network
Hidden layers: 2, Nodes per layer: 100, Activation functions: ReLU, Epochs: 1000 or 5000, Learning rate: 0.001, Batch training, Early stopping steps: 100 |  Simultaneous Quantile Regression [48]
Same as point prediction baseline, Deep Neural Network

###  4.3 Metrics

The evaluation of point predictions employs Root Mean Squared Error (RMSE),
Mean Absolute Error (MAE), and R-squared ($\text{R}^{2}$). For prediction
intervals, the chosen metrics, following [52, 53], are the Prediction Interval
Coverage Probability (PICP), Mean Prediction Interval Width (MPIW), and
Normalized Mean Prediction Interval Width (NMPIW) as described in Subsection
D.1. This paper also evaluates SEMF on our new metric, termed Coverage-Width
Ratio (CWR),

$\text{CWR}=\frac{\text{PICP}}{\text{NMPIW}},$ (19)

which evaluates the coverage probability ratio to the prediction interval’s
width. CWR provides a refined understanding of the balance between an
interval’s accuracy (coverage) and precision (width). Though a larger value of
this metric is better in higher confidence levels, the marginal increase in
NMPIW is likely higher than that of PICP, resulting in decreasing CWR.

In our case, measuring the performance of SEMF over the baseline models is far
more critical than reviewing absolute metrics in isolation. For any metric
above, except for ($\text{R}^{2}$), this is computed as

$\text{Metric}_{\Delta}(\%)=\left(\frac{\text{Metric}_{\text{SEMF}}-\text{Metric}_{\text{Baseline}}}{\text{Metric}_{\text{Baseline}}}\right)\times
100,$ (20)

on which we base our decisions for selecting the best hyper-parameters as
explained in Subsection D.2.

##  5 Results

We trained and tested 330 models corresponding to the three model
types—MultiXGBs, MultiETs, and MultiMLPs—across 11 datasets with both complete
and 50% missing data, with five seeds for each combination. Table 2 and Table
3 provide the mean and standard deviation for our metrics aggregated over the
five seeds. All the results are in Appendix E. Below, we discuss the
performances with complete data and with incomplete data.

###  5.1 Complete Data

The results of models with complete data are in Table 2. Overall, MultiXGBs
and MultiMLPs performed generally well in producing intervals compared to
their baselines as shown by positive $\Delta$CWR and $\Delta$PICP, notably,
for instance, on cpu_activity, naval_propulsion_plant, and energy_efficiency.
Interestingly, MultiMLPs also attained good relative performance improvements
for point prediction, which can indicate that the chosen $\sigma_{k}$ was too
low for generating performant prediction intervals. MultiETs attained mixed
results as, despite significant improvements in $\Delta$CWR, it often failed
to beat the baseline on $\Delta$PICP.

Table 2: Test results for all models with complete data at 95% quantiles
aggregated over five seeds. For each metric, the mean and standard deviation
of the performance across the seeds are separated by ±. Performance over the
baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
MultiXGBs |  |  |  |  |  |  |  |  |
space_ga | -1%±6% | 6%±2% | -7%±6% | 0.89±0.01 | 1.80±0.12 | 0.20±0.01 | -9%±1% | -10%±2% | 0.60±0.01
cpu_activity | 25%±7% | 9%±1% | 12%±5% | 0.89±0.01 | 0.37±0.01 | 0.07±0.00 | 21%±1% | -1%±3% | 0.98±0.00
naval_propulsion_plant | 156%±15% | 8%±1% | 58%±2% | 0.96±0.00 | 0.39±0.01 | 0.11±0.00 | -14%±4% | -12%±2% | 0.99±0.00
miami_housing | 61%±6% | -7%±2% | 43%±3% | 0.82±0.01 | 0.45±0.01 | 0.06±0.00 | 4%±5% | 3%±3% | 0.91±0.01
kin8nm | 10%±3% | 5%±2% | 4%±4% | 0.90±0.01 | 1.94±0.07 | 0.38±0.01 | -20%±2% | -22%±2% | 0.63±0.01
concrete_compressive_strength | 13%±5% | 26%±6% | -11%±8% | 0.91±0.01 | 1.25±0.03 | 0.25±0.01 | -26%±7% | -40%±13% | 0.85±0.01
cars | 18%±9% | 15%±4% | 2%±9% | 0.89±0.02 | 0.70±0.05 | 0.12±0.01 | -3%±3% | -1%±4% | 0.95±0.00
energy_efficiency | 247%±99% | 22%±7% | 60%±20% | 0.88±0.02 | 0.16±0.04 | 0.05±0.01 | -15%±25% | -16%±19% | 1.00±0.00
california_housing | 31%±5% | 0%±1% | 23%±4% | 0.88±0.01 | 1.17±0.02 | 0.28±0.00 | -1%±1% | -1%±1% | 0.81±0.00
airfoil_self_noise | 61%±45% | 4%±10% | 31%±18% | 0.81±0.07 | 0.89±0.22 | 0.19±0.05 | -64%±36% | -73%±46% | 0.86±0.05
QSAR_fish_toxicity | 10%±4% | 11%±4% | -1%±3% | 0.76±0.01 | 1.46±0.03 | 0.23±0.00 | 3%±4% | 1%±5% | 0.55±0.02
MultiETs |  |  |  |  |  |  |  |  |
space_ga | 9%±4% | -2%±1% | 9%±4% | 0.91±0.01 | 2.04±0.07 | 0.22±0.01 | -16%±2% | -18%±3% | 0.54±0.02
cpu_activity | 27%±3% | -8%±0% | 27%±2% | 0.90±0.01 | 0.47±0.01 | 0.09±0.00 | -14%±2% | -20%±2% | 0.98±0.00
naval_propulsion_plant | 219%±36% | -9%±0% | 71%±3% | 0.91±0.00 | 0.66±0.07 | 0.19±0.02 | -320%±45% | -409%±40% | 0.96±0.01
miami_housing | 69%±4% | -13%±0% | 48%±1% | 0.86±0.00 | 0.65±0.02 | 0.08±0.00 | -10%±2% | -20%±2% | 0.90±0.00
kin8nm | 10%±2% | -9%±0% | 17%±2% | 0.88±0.01 | 2.27±0.05 | 0.44±0.01 | -35%±2% | -38%±2% | 0.48±0.01
concrete_compressive_strength | 9%±10% | -12%±2% | 19%±9% | 0.81±0.02 | 1.18±0.11 | 0.24±0.02 | -67%±7% | -94%±10% | 0.76±0.02
cars | 11%±7% | 16%±2% | -5%±7% | 0.83±0.02 | 0.53±0.04 | 0.09±0.01 | 3%±4% | 1%±2% | 0.95±0.00
energy_efficiency | -4%±6% | 0%±8% | -5%±15% | 0.65±0.04 | 0.08±0.01 | 0.02±0.00 | 3%±4% | 0%±2% | 1.00±0.00
california_housing | 31%±2% | -7%±0% | 29%±1% | 0.90±0.00 | 1.70±0.03 | 0.41±0.01 | -15%±1% | -21%±1% | 0.71±0.01
airfoil_self_noise | 18%±28% | -15%±5% | 25%±18% | 0.83±0.05 | 1.22±0.31 | 0.26±0.07 | -117%±50% | -141%±57% | 0.80±0.09
QSAR_fish_toxicity | 12%±2% | -9%±1% | 19%±2% | 0.82±0.01 | 1.81±0.02 | 0.28±0.00 | -6%±1% | -11%±1% | 0.53±0.01
MultiMLPs |  |  |  |  |  |  |  |  |
space_ga | 8%±1% | 0%±3% | 7%±3% | 0.81±0.02 | 1.25±0.02 | 0.14±0.00 | 0%±1% | -1%±1% | 0.75±0.01
cpu_activity | 8%±3% | -9%±4% | 15%±6% | 0.72±0.02 | 0.24±0.01 | 0.05±0.00 | 7%±2% | 5%±2% | 0.98±0.00
naval_propulsion_plant | 21%±20% | 0%±3% | 16%±15% | 0.92±0.01 | 0.22±0.03 | 0.07±0.01 | -45%±28% | -36%±23% | 1.00±0.00
miami_housing | 9%±4% | 2%±3% | 6%±6% | 0.81±0.02 | 0.41±0.03 | 0.05±0.00 | -3%±2% | 4%±2% | 0.91±0.00
kin8nm | 2%±8% | 8%±6% | -6%±12% | 0.81±0.02 | 0.68±0.03 | 0.13±0.01 | 7%±3% | 7%±2% | 0.93±0.00
concrete_compressive_strength | 98%±15% | -23%±8% | 61%±5% | 0.56±0.04 | 0.28±0.01 | 0.06±0.00 | 12%±6% | 15%±4% | 0.91±0.01
cars | 0%±15% | 6%±13% | -9%±33% | 0.76±0.09 | 0.40±0.10 | 0.07±0.02 | -3%±3% | 0%±3% | 0.95±0.00
energy_efficiency | 82%±20% | 25%±26% | 31%±12% | 0.62±0.05 | 0.06±0.00 | 0.02±0.00 | 32%±4% | 33%±3% | 1.00±0.00
california_housing | -14%±3% | 8%±1% | -26%±4% | 0.89±0.01 | 1.31±0.02 | 0.31±0.01 | 7%±0% | 9%±1% | 0.82±0.00
airfoil_self_noise | 74%±46% | 0%±7% | 39%±17% | 0.75±0.02 | 0.35±0.01 | 0.08±0.00 | 18%±4% | 16%±5% | 0.97±0.00
QSAR_fish_toxicity | -6%±4% | 11%±4% | -19%±5% | 0.76±0.04 | 1.49±0.06 | 0.23±0.01 | 4%±5% | 7%±5% | 0.55±0.04

###  5.2 Missing Data

The results of models with 50% missing data are in Table 3. Note that relative
metrics for SEMF are compared with the best result from any baseline imputer
on that metric, regardless of how the imputer performed on the other metrics.
Overall, the results are less consistent than with complete data. MultiXGBs
maintained strong performance on some datasets, such as
naval_propulsion_plant, but declined on others, like cpu_activity. MultiETs
exhibited narrow predictions, improved absolute PICP in some cases, such as
energy_efficiency, and improved relative point predictions. Conversely,
MultiMLPs, while offering narrow prediction intervals, had higher PICP,
indicating decreased predictive power and increased model uncertainty.

Table 3: Test results for all models with missing data at 95% quantiles
aggregated over five seeds. For each metric, the mean and standard deviation
of the performance across the seeds are separated by ±. Performance over the
baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
MultiXGBs |  |  |  |  |  |  |  |  |
space_ga | -5%±6% | 4%±2% | -12%±7% | 0.86±0.02 | 2.05±0.16 | 0.22±0.02 | -7%±5% | -9%±3% | 0.45±0.07
cpu_activity | -5%±21% | 12%±9% | -37%±52% | 0.89±0.07 | 0.76±0.24 | 0.14±0.05 | -1%±24% | -17%±14% | 0.83±0.08
naval_propulsion_plant | 105%±33% | -16%±7% | 56%±6% | 0.70±0.07 | 0.62±0.06 | 0.18±0.02 | -14%±22% | -29%±39% | 0.77±0.06
miami_housing | 16%±14% | -5%±6% | 10%±18% | 0.82±0.07 | 0.77±0.15 | 0.10±0.02 | -14%±12% | -18%±11% | 0.72±0.11
kin8nm | -5%±2% | 3%±2% | -14%±5% | 0.88±0.01 | 2.46±0.07 | 0.48±0.01 | -3%±3% | -4%±5% | 0.40±0.02
concrete_compressive_strength | 2%±3% | 8%±4% | -10%±4% | 0.77±0.04 | 1.47±0.03 | 0.30±0.01 | -16%±14% | -20%±16% | 0.58±0.09
cars | -15%±15% | 26%±10% | -63%±29% | 0.85±0.04 | 1.58±0.18 | 0.27±0.03 | -30%±25% | -35%±27% | 0.66±0.17
energy_efficiency | 10%±29% | 11%±6% | -24%±48% | 0.92±0.03 | 0.77±0.19 | 0.22±0.05 | -186%±296% | -153%±210% | 0.91±0.09
california_housing | 26%±8% | -2%±3% | 21%±7% | 0.88±0.03 | 1.67±0.12 | 0.40±0.03 | -4%±4% | -5%±5% | 0.69±0.05
airfoil_self_noise | -2%±21% | -3%±4% | -4%±19% | 0.73±0.05 | 1.59±0.33 | 0.34±0.07 | -37%±34% | -43%±40% | 0.41±0.13
QSAR_fish_toxicity | -2%±3% | 0%±5% | -6%±8% | 0.73±0.03 | 1.70±0.13 | 0.26±0.02 | -4%±5% | -2%±4% | 0.36±0.05
MultiETs |  |  |  |  |  |  |  |  |
space_ga | 4%±4% | -5%±1% | 7%±4% | 0.88±0.02 | 2.23±0.10 | 0.24±0.01 | -12%±4% | -14%±5% | 0.40±0.06
cpu_activity | -16%±12% | -5%±2% | -20%±21% | 0.93±0.01 | 1.07±0.20 | 0.20±0.04 | -30%±26% | -54%±22% | 0.79±0.11
naval_propulsion_plant | 58%±5% | -17%±5% | 46%±3% | 0.82±0.05 | 1.38±0.10 | 0.41±0.03 | -36%±23% | -93%±53% | 0.71±0.06
miami_housing | 20%±16% | -10%±5% | 21%±14% | 0.87±0.05 | 1.15±0.26 | 0.14±0.03 | -27%±12% | -37%±12% | 0.67±0.12
kin8nm | 5%±3% | -8%±2% | 12%±4% | 0.87±0.02 | 2.55±0.10 | 0.50±0.02 | -11%±3% | -13%±3% | 0.32±0.03
concrete_compressive_strength | 1%±6% | -14%±4% | 13%±2% | 0.80±0.02 | 1.83±0.06 | 0.37±0.01 | -35%±23% | -50%±28% | 0.49±0.07
cars | -23%±7% | -8%±6% | -35%±16% | 0.83±0.04 | 1.41±0.15 | 0.24±0.02 | -33%±20% | -31%±20% | 0.66±0.10
energy_efficiency | -15%±18% | 0%±10% | -36%±47% | 0.85±0.06 | 0.51±0.05 | 0.15±0.01 | -73%±82% | -84%±72% | 0.94±0.04
california_housing | 21%±7% | -7%±2% | 22%±5% | 0.91±0.02 | 2.19±0.13 | 0.52±0.03 | -14%±3% | -17%±4% | 0.61±0.06
airfoil_self_noise | -7%±20% | -20%±2% | 8%±18% | 0.76±0.03 | 1.96±0.39 | 0.42±0.08 | -65%±60% | -93%±78% | 0.33±0.12
QSAR_fish_toxicity | 10%±8% | -9%±3% | 16%±4% | 0.82±0.04 | 2.12±0.20 | 0.33±0.03 | -8%±9% | -10%±10% | 0.39±0.10
MultiMLPs |  |  |  |  |  |  |  |  |
space_ga | -11%±6% | -15%±2% | 1%±7% | 0.68±0.04 | 1.44±0.09 | 0.16±0.01 | -40%±26% | -34%±17% | 0.21±0.25
cpu_activity | -50%±18% | 22%±16% | -193%±132% | 0.86±0.11 | 0.97±0.42 | 0.18±0.08 | -15%±21% | -23%±19% | 0.75±0.14
naval_propulsion_plant | -52%±21% | -14%±27% | -102%±62% | 0.63±0.20 | 1.17±0.28 | 0.34±0.08 | -108%±74% | -183%±218% | -0.03±0.53
miami_housing | -33%±17% | 13%±11% | -92%±78% | 0.85±0.06 | 0.99±0.32 | 0.12±0.04 | -25%±9% | -28%±12% | 0.67±0.11
kin8nm | -11%±6% | -5%±2% | -8%±6% | 0.77±0.01 | 1.76±0.09 | 0.35±0.02 | -16%±9% | -15%±10% | 0.37±0.04
concrete_compressive_strength | -15%±14% | -13%±31% | -10%±54% | 0.60±0.25 | 1.11±0.59 | 0.22±0.12 | -12%±6% | -12%±7% | 0.54±0.06
cars | -51%±5% | 43%±8% | -211%±38% | 0.93±0.02 | 1.74±0.22 | 0.29±0.04 | -28%±14% | -26%±15% | 0.66±0.11
energy_efficiency | -13%±20% | 9%±17% | -30%±28% | 0.66±0.15 | 0.33±0.15 | 0.09±0.04 | 1%±10% | -1%±18% | 0.95±0.03
california_housing | -17%±7% | 4%±4% | -28%±16% | 0.88±0.03 | 1.84±0.17 | 0.44±0.04 | -5%±3% | -3%±4% | 0.61±0.06
airfoil_self_noise | 2%±18% | -2%±6% | -1%±23% | 0.75±0.08 | 1.15±0.41 | 0.25±0.09 | -15%±9% | -9%±9% | 0.55±0.20
QSAR_fish_toxicity | -17%±5% | 6%±7% | -30%±12% | 0.76±0.04 | 1.82±0.18 | 0.28±0.03 | -13%±7% | -14%±9% | 0.34±0.07

###  5.3 Discussion

The results indicate that specific models, particularly MultiXGBs, perform
strongly with complete data. However, their robustness diminishes with
increasing amounts of missing data. Further, we observe that MultiETs can
produce narrow prediction intervals at the expense of reduced coverage. In
contrast, despite the experiment design prioritizing interval estimation over
point accuracy, MultiMLPs also deliver good point predictions. Precise
choices, such as processing the columns separately, using the same hyper-
parameters for both the complete and incomplete data, and fixing the same
number of latent nodes per input ($m_{k}$), might have further constrained the
models’ ability to tailor their performance to varying dataset complexities.

The stability of SEMF across different seeds, as indicated by the absolute
PICP and NMPIW, contrasts with the significant variability observed in
baseline models, suggesting that SEMF offers a more reliable performance
framework. The method’s sampling operation, akin to cross-validation, helps
combat overfitting, ensuring robust results despite variations in the data.
Nonetheless, most models do not achieve the desired 95% prediction intervals.
This shortfall could be addressed by integrating conformal prediction
techniques with the intervals generated by SEMF, enhancing the reliability and
coverage of prediction intervals across different datasets and missing data
scenarios [34].

##  6 Conclusion

This paper introduces the Supervised Expectation-Maximization Framework
(SEMF), a novel model-agnostic approach for generating prediction intervals in
datasets with missing values. SEMF draws from the EM algorithm for supervised
learning to devise latent representations that produce better prediction
intervals than quantile regression. Due to SEMF’s iterative simulation
technique, training and inference can be done with complete and incomplete
data. A comprehensive set of 330 experimental runs on 11 datasets with three
different model types showed that SEMF is superior to quantile regression,
particularly for models such as XGBoost, which intrinsically lacks latent
representations. This research underscores the potential of SEMF in various
application domains and opens new avenues for further exploration for
supervised latent representation learning and uncertainty estimation.

##  7 Limitations & Future Work

The primary limitation of this study was its reliance on the normality
assumption, which may not fully capture the potential of SEMF across diverse
data distributions. Although in Appendix F we demonstrate that the framework
can learn non-normal patterns, further investigation and exploration of SEMF
under other distributions, such as uniform and log-normal, are needed. The
computational complexity of the approach presents another significant
challenge, as the current implementation can be optimized for large-scale
applications. Additionally, while the CWR metric is useful, it implicitly
assumes that a 1% drop in PICP equates to a 1% reduction in NMPIW, thus
assuming a uniform distribution. Evaluating CWR under various distributional
assumptions would provide a more comprehensive assessment of its implications.
Finally, applying the same hyper-parameters across all datasets without
specific tuning for incomplete data likely contributes to the observed decline
in accuracy and robustness.

Future work presents several intriguing avenues for exploration. A promising
direction is the application of SEMF in multi-modal data settings, where the
distinct $p_{\phi}$ components of the framework could be adapted to process
diverse data types—from images and text to tabular datasets—enabling a more
nuanced and powerful approach to integrating heterogeneous data sources. This
capability positions the framework as a versatile tool for addressing missing
data challenges across various domains and can also help expand it to discrete
and multiple outputs. Another valuable area for development is the exploration
of methods to capture and leverage dependencies among input features, which
could improve the model’s predictive performance and provide deeper insights
into the underlying data structure. These advancements can enhance the broader
appeal of end-to-end approaches like SEMF in the ML community.

## References

  * [1] Michael W. Dusenberry, Dustin Tran, Edward Choi, Jonas Kemp, Jeremy Nixon, Ghassen Jerfel, Katherine Heller, and Andrew M. Dai.  Analyzing the role of model uncertainty for electronic health records.  In Proceedings of the ACM Conference on Health, Inference, and Learning, CHIL ’20, page 204–213, New York, NY, USA, 2020. Association for Computing Machinery.
  * [2] Wojciech Wisniewski, David Lindsay, and Sian Lindsay.  Application of conformal prediction interval estimations to market makers’ net positions.  In Alexander Gammerman, Vladimir Vovk, Zhiyuan Luo, Evgueni Smirnov, and Giovanni Cherubin, editors, Proceedings of the Ninth Symposium on Conformal and Probabilistic Prediction and Applications, volume 128 of Proceedings of Machine Learning Research, pages 285–301. PMLR, 09–11 Sep 2020.
  * [3] Xiaolin Tang, Kai Yang, Hong Wang, Jiahang Wu, Yechen Qin, Wenhao Yu, and Dongpu Cao.  Prediction-uncertainty-aware decision-making for autonomous vehicles.  IEEE Transactions on Intelligent Vehicles, 7(4):849–862, 2022.
  * [4] Zoubin Ghahramani.  Probabilistic machine learning and artificial intelligence.  Nature, 521(7553):452–459, 2015.
  * [5] A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society. Series B (Methodological), 39(1):1–38, 1977.
  * [6] Alp Kucukelbir David M. Blei and Jon D. McAuliffe.  Variational inference: A review for statisticians.  Journal of the American Statistical Association, 112(518):859–877, 2017.
  * [7] Diederik P. Kingma and Max Welling.  Auto-Encoding Variational Bayes.  In 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014.
  * [8] Jakob Gawlikowski, Cedrique Rovile Njieutcheu Tassi, Mohsin Ali, Jongseok Lee, Matthias Humt, Jianxiang Feng, Anna Kruspe, Rudolph Triebel, Peter Jung, Ribana Roscher, et al.  A survey of uncertainty in deep neural networks.  Artificial Intelligence Review, 56(Suppl 1):1513–1589, 2023.
  * [9] CF Jeff Wu.  On the convergence properties of the em algorithm.  The Annals of statistics, pages 95–103, 1983.
  * [10] Zoubin Ghahramani and Michael Jordan.  Supervised learning from incomplete data via an em approach.  In J. Cowan, G. Tesauro, and J. Alspector, editors, Advances in Neural Information Processing Systems, volume 6. Morgan-Kaufmann, 1993.
  * [11] David Williams, Xuejun Liao, Ya Xue, and Lawrence Carin.  Incomplete-data classification using logistic regression.  In Proceedings of the 22nd International Conference on Machine Learning, ICML ’05, page 972–979, New York, NY, USA, 2005. Association for Computing Machinery.
  * [12] Wenchong He and Zhe Jiang.  Semi-supervised learning with the em algorithm: A comparative study between unstructured and structured prediction.  IEEE Transactions on Knowledge and Data Engineering, 34(6):2912–2920, 2022.
  * [13] Robin Louiset, Pietro Gori, Benoit Dufumier, Josselin Houenou, Antoine Grigis, and Edouard Duchesnay.  Ucsl: A machine learning expectation-maximization framework for unsupervised clustering driven by supervised learning.  In Machine Learning and Knowledge Discovery in Databases. Research Track: European Conference, ECML PKDD 2021, Bilbao, Spain, September 13–17, 2021, Proceedings, Part I 21, pages 755–771. Springer, 2021.
  * [14] Alp Kucukelbir David M. Blei and Jon D. McAuliffe.  Variational inference: A review for statisticians.  Journal of the American Statistical Association, 112(518):859–877, 2017.
  * [15] Sivaraman Balakrishnan, Martin J. Wainwright, and Bin Yu.  Statistical guarantees for the EM algorithm: From population to sample-based analysis.  The Annals of Statistics, 45(1):77 – 120, 2017.
  * [16] Durk P Kingma, Shakir Mohamed, Danilo Jimenez Rezende, and Max Welling.  Semi-supervised learning with deep generative models.  Advances in neural information processing systems, 27, 2014.
  * [17] Fuzhen Zhuang, Xiaohu Cheng, Ping Luo, Sinno Jialin Pan, and Qing He.  Supervised representation learning: transfer learning with deep autoencoders.  In Twenty-fourth international joint conference on artificial intelligence, IJCAI’15, page 4119–4125. AAAI Press, 2015.
  * [18] Lei Le, Andrew Patterson, and Martha White.  Supervised autoencoders: Improving generalization performance with unsupervised regularizers.  Advances in neural information processing systems, 31, 2018.
  * [19] Tianchen Ji, Srikanth Vuppala, Girish V. Chowdhary, and K. Driggs-Campbell.  Multi-modal anomaly detection for unstructured and uncertain environments.  In Conference on Robot Learning, 2020.
  * [20] Yilin Zhuang, Zhuobin Zhou, Burak Alakent, and Mehmet Mercangöz.  Semi-supervised variational autoencoders for regression: Application to soft sensors.  In 2023 IEEE 21st International Conference on Industrial Informatics (INDIN), pages 1–8, 2023.
  * [21] Zhekun Luo, Devin Guillory, Baifeng Shi, Wei Ke, Fang Wan, Trevor Darrell, and Huijuan Xu.  Weakly-supervised action localization with expectation-maximization multi-instance learning.  In Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXIX 16, pages 729–745. Springer, 2020.
  * [22] Moucheng Xu, Yukun Zhou, Chen Jin, Marius de Groot, Daniel C Alexander, Neil P Oxtoby, Yipeng Hu, and Joseph Jacob.  Expectation maximization pseudo labels.  Medical Image Analysis, page 103125, 2024.
  * [23] Christopher Williams and Carl Rasmussen.  Gaussian processes for regression.  Advances in neural information processing systems, 8, 1995.
  * [24] James Hensman, Alexander Matthews, and Zoubin Ghahramani.  Scalable variational gaussian process classification.  In Artificial Intelligence and Statistics, pages 351–360. PMLR, 2015.
  * [25] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In Maria Florina Balcan and Kilian Q. Weinberger, editors, Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1050–1059, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * [26] Leo Breiman.  Random forests.  Machine learning, 45:5–32, 2001.
  * [27] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  Advances in neural information processing systems, 30, 2017.
  * [28] Roger Koenker and Gilbert Bassett.  Regression quantiles.  Econometrica, 46(1):33–50, 1978.
  * [29] Roger Koenker and Kevin F. Hallock.  Quantile regression.  The Journal of Economic Perspectives, 15(4):143–156, 2001.
  * [30] Vladimir Vovk, Alexander Gammerman, and Glenn Shafer.  Algorithmic learning in a random world, volume 29.  Springer, 2005.
  * [31] Vladimir Vovk, Alexander Gammerman, and Glenn Shafer.  Conformal prediction: General case and regression.  In Algorithmic Learning in a Random World, pages 19–69. Springer, 2022.
  * [32] Ingo Steinwart and Andreas Christmann.  Estimating conditional quantiles with the help of the pinball loss.  Bernoulli, 17(1), February 2011.
  * [33] Nicolai Meinshausen and Greg Ridgeway.  Quantile regression forests.  Journal of machine learning research, 7(6), 2006.
  * [34] Yaniv Romano, Evan Patterson, and Emmanuel Candes.  Conformalized quantile regression.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [35] Jinsung Yoon, Yao Zhang, James Jordon, and Mihaela van der Schaar.  Vime: Extending the success of self- and semi-supervised learning to tabular domain.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 11033–11043. Curran Associates, Inc., 2020.
  * [36] Anil Jadhav, Dhanya Pramod, and Krishnan Ramanathan.  Comparison of performance of data imputation methods for numeric dataset.  Applied Artificial Intelligence, 33(10):913–933, 2019.
  * [37] Joseph G Ibrahim, Hongtu Zhu, and Niansheng Tang.  Model selection criteria for missing-data problems using the em algorithm.  Journal of the American Statistical Association, 103(484):1648–1658, 2008.
  * [38] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [39] Stef van Buuren and Karin Groothuis-Oudshoorn.  mice: Multivariate imputation by chained equations in r.  Journal of Statistical Software, 45(3):1–67, 2011.
  * [40] S. F. Buck.  A method of estimation of missing values in multivariate data suitable for use with an electronic computer.  Journal of the Royal Statistical Society. Series B (Methodological), 22(2):302–306, 1960.
  * [41] Joseph L Schafer.  Analysis of incomplete multivariate data.  CRC press, 1997.
  * [42] Siwei Wang, Miaomiao Li, Ning Hu, En Zhu, Jingtao Hu, Xinwang Liu, and Jianping Yin.  K-means clustering with incomplete data.  IEEE Access, 7:69162–69171, 2019.
  * [43] Ricardo Cardoso Pereira, Miriam Seoane Santos, Pedro Pereira Rodrigues, and Pedro Henriques Abreu.  Reviewing autoencoders for missing data imputation: Technical trends, applications and outcomes.  Journal of Artificial Intelligence Research, 69:1255–1285, 2020.
  * [44] Marek Śmieja, Łukasz Struski, Jacek Tabor, Bartosz Zieliński, and Przemysław Spurek.  Processing of missing data by neural networks.  Advances in neural information processing systems, 31, 2018.
  * [45] Sebastian Felix Fischer, Matthias Feurer, and Bernd Bischl.  OpenML-CTR23 – a curated tabular regression benchmarking suite.  In AutoML Conference 2023 (Workshop), 2023.
  * [46] Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.
  * [47] Pierre Geurts, Damien Ernst, and Louis Wehenkel.  Extremely randomized trees.  Machine learning, 63:3–42, 2006.
  * [48] Natasa Tagasovska and David Lopez-Paz.  Single-model uncertainties for deep learning.  Advances in Neural Information Processing Systems, 32, 2019.
  * [49] Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter.  Self-normalizing neural networks.  Advances in neural information processing systems, 30, 2017.
  * [50] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785–794, 2016.
  * [51] Reid A. Johnson.  quantile-forest: A python package for quantile regression forests.  Journal of Open Source Software, 9(93):5976, 2024.
  * [52] Tim Pearce, Alexandra Brintrup, Mohamed Zaki, and Andy Neely.  High-quality prediction intervals for deep learning: A distribution-free, ensembled approach.  In International conference on machine learning, pages 4075–4084. PMLR, 2018.
  * [53] Ting Zhou, Yuxin Jie, Yingjie Wei, Yanyi Zhang, and Hui Chen.  A real-time prediction interval correction method with an unscented kalman filter for settlement monitoring of a power station dam.  Scientific Reports, 13(1):4055, 2023.
  * [54] Kaggle.  Moneyball dataset.  https://www.kaggle.com/datasets/wduckett/moneyball-mlb-stats-19622012, 2017.  Accessed: 2024-02-02.
  * [55] R. Kelley Pace and Ronald Barry.  Quick computation of spatial autoregressive estimators.  Geographical Analysis, 29(3):232–247, 1997.
  * [56] C. Rasmussen, R. Neal, G. Hinton, D. van Camp, M. Revow, Z. Ghahramani, R. Kustra, and R. Tibshirani.  Computer activity dataset.  http://www.cs.toronto.edu/~delve/data/datasets.html, 1996.  Accessed: 2024-02-02.
  * [57] Andrea Coraddu, Luca Oneto, Aessandro Ghio, Stefano Savio, Davide Anguita, and Massimo Figari.  Machine learning approaches for improving condition-based maintenance of naval propulsion plants.  Proceedings of the Institution of Mechanical Engineers, Part M: Journal of Engineering for the Maritime Environment, 230(1):136–153, 2016.
  * [58] Kaggle.  Miami housing dataset.  https://www.kaggle.com/datasets/deepcontractor/miami-housing-dataset, 2022.  Accessed: 2024-02-02.
  * [59] Z. Ghahramani.  The kin datasets.  https://www.cs.toronto.edu/~delve/data/kin/desc.html, 1996.  Accessed: 2024-02-02.
  * [60] I.-C. Yeh.  Modeling of strength of high-performance concrete using artificial neural networks.  Cement and Concrete Research, 28(12):1797–1808, 1998.
  * [61] Shonda Kuiper.  Introduction to multiple regression: How much is your car worth?  Journal of Statistics Education, 16(3), 2008.
  * [62] Athanasios Tsanas and Angeliki Xifara.  Accurate quantitative estimation of energy performance of residential buildings using statistical machine learning tools.  Energy and Buildings, 49:560–567, 2012.
  * [63] R. Kelley Pace and Ronald Barry.  Sparse spatial autoregressions.  Statistics & Probability Letters, 33(3):291–297, 1997.
  * [64] Thomas Brooks, Dennis Pope, and Michael Marcolini.  Airfoil self-noise and prediction.  NASA Technical Report 1218, NASA, 1989.
  * [65] R. Todeschini M. Cassotti, D. Ballabio and V. Consonni.  A similarity-based qsar model for predicting acute toxicity towards the fathead minnow (pimephales promelas).  SAR and QSAR in Environmental Research, 26(3):217–243, 2015.  PMID: 25780951.
  * [66] Lukas Biewald.  Experiment tracking with weights and biases, 2020.  Software available from wandb.com.

##  Appendix A SEMF algorithm

Algorithm 1 SEMF Training: two input sources where $x_{1}$ can be missing

1:$y,x_{1},x_{2},R$

2:$\theta,\phi_{1},\phi_{2},\xi$

3:Initialize $\theta,\phi_{1},\phi_{2},\xi$

4:Initialize $D_{y},D_{z_{1}},D_{z_{2}},D_{x},D_{i}$ to $\emptyset$

5:Split $I=\\{1,\ldots,N\\}$ into $L$ batches $\\{b_{1},\ldots,b_{L}\\}$

6:for $\ell=1$, …, $L$ do

7: for all $i$ in $b_{\ell}$ do

8: if $x_{1,i}$ is absent then

9: for $r=1$, …, $R$ do

10: Simulate $[j_{i,r},x_{1,i,r}]\sim p_{\xi}(\cdot|x_{2,i})$

11: Simulate $z_{1,i,r}\sim p_{\phi_{1}}(\cdot|x_{1,i,r})$

12: Simulate $z_{2,i,r}\sim p_{\phi_{2}}(\cdot|x_{2,i})$

13: Set $z_{i,r}=[z_{1,i,r},z_{2,i,r}]$

14: end for

15: else

16: for $r=1$, …, $R$ do

17: Simulate $z_{1,i,r}\sim p_{\phi_{1}}(\cdot|x_{1,i})$

18: Simulate $z_{2,i,r}\sim p_{\phi_{2}}(\cdot|x_{2,i})$

19: Set $z_{i,r}=[z_{1,i,r},z_{2,i,r}]$

20: end for

21: end if

22: for $r=1$, …, $R$ do

23: Compute
$w_{i,r}=\frac{p_{\theta}(y_{i}|z_{i,r})}{\sum_{t=1}^{R}p_{\theta}(y_{i}|z_{i,t})}$

24: Update $D_{y}\leftarrow D_{y}\cup[y_{i}|z_{i,r}|w_{i,r}]$

25: Update $D_{z_{2}}\leftarrow D_{z_{2}}\cup[z_{2,i,r}|x_{2,i}|w_{i,r}]$

26: if $x_{1,i}$ is absent then

27: Update $D_{z_{1}}\leftarrow D_{z_{1}}\cup[z_{1,i,r}|x_{1,i,r}|w_{i,r}]$

28: Update $D_{x}\leftarrow D_{x}\cup[j_{i,r}|x_{2,i}|w_{i,r}]$

29: else

30: Update $D_{z_{1}}\leftarrow D_{z_{1}}\cup[z_{1,i,r}|x_{1,i}|w_{i,r}]$

31: end if

32: end for

33: Update $\theta\leftarrow Q_{y}(\theta,D_{y})$

34: Update $\phi_{1}\leftarrow Q_{1}(\phi_{1},D_{z_{1}})$

35: Update $\phi_{2}\leftarrow Q_{2}(\phi_{2},D_{z_{2}})$

36: Update $\xi\leftarrow Q_{x}(\xi,D_{x})$

37: end for

38:end for

39:Check convergence; Go to step 4 if not

Algorithm 2 SEMF Inference

1:$\theta^{*},\phi_{1}^{*},\phi_{2}^{*},\xi^{*},x_{1},x_{2},R$

2:$z_{i,r}$

3:for $r=1,\ldots,R$ do

4: if $x_{1,i}$ is absent then

5: Simulate $[j_{i,r},x_{1,i,r}]\sim p_{\xi^{*}}(\cdot|x_{2,i})$

6: Simulate $z_{1,i,r}\sim p_{\phi_{1}^{*}}(\cdot|x_{1,i,r})$

7: Simulate $z_{2,i,r}\sim p_{\phi_{2}^{*}}(\cdot|x_{2,i})$

8: Set $z_{i,r}=[z_{1,i,r},z_{2,i,r}]$

9: else

10: Simulate $z_{1,i,r}\sim p_{\phi_{1}^{*}}(\cdot|x_{1,i})$

11: Simulate $z_{2,i,r}\sim p_{\phi_{2}^{*}}(\cdot|x_{2,i})$

12: Set $z_{i,r}=[z_{1,i,r},z_{2,i,r}]$

13: end if

14:end for

##  Appendix B Datasets for tabular benchmark

OpenML-CTR23 [45] datasets are selected in the following manner. The first
criterion is to exclude datasets exceeding 30,000 instances or 30 features to
maintain computational tractability. Moreover, we exclude the moneyball data
[54] to control for missing values and any datasets with non-numeric features,
such as those with temporal or ordinal data not encoded numerically. We then
categorize the datasets based on size: small for those with less than ten
features, medium for 10 to 19 features, and large for 20 to 29 features. We
apply a similar size classification based on the number of instances,
considering datasets with more than 10,000 instances as large. To avoid
computational constraints, we exclude datasets that were large in both
features and instances, ensuring a varied yet manageable set for our
experiments. This leads us to the final list of 11 datasets listed in Table 4.

Table 4: Summary of benchmark tabular datasets retained from [45]

Dataset Name | N Samples | N Features | OpenML Data ID | Y [Min:Max] | Source
---|---|---|---|---|---
space_ga | 3,107 | 7 | 45402 | [-3.06:0.1] | [55]
cpu_activity | 8,192 | 22 | 44978 | [0:99] | [56]
naval_propulsion_plant | 11,934 | 15 | 44969 | [0.95:1.0] | [57]
miami_housing | 13,932 | 16 | 44983 | [72,000:2,650,000] | [58]
kin8nm | 8,192 | 9 | 44980 | [0.04:1.46] | [59]
concrete_compressive_strength | 1,030 | 9 | 44959 | [2.33:82.6] | [60]
cars | 804 | 18 | 44994 | [8,639:70,756] | [61]
energy_efficiency | 768 | 9 | 44960 | [6.01:43.1] | [62]
california_housing | 20,640 | 9 | 44977 | [14,999:500,001] | [63]
airfoil_self_noise | 1,503 | 6 | 44957 | [103.38:140.98] | [64]
QSAR_fish_toxicity | 908 | 7 | 44970 | [0.053:9.612] | [65]

##  Appendix C Optimal set of hyper-parameters

The hyper-parameter tuning for SEMF is implemented and monitored using Weights
& Biases [66]. A random search is done in the hyper-parameter space for a
maximum of 500 iterations on all 11 datasets, focusing on tuning the models
only on the complete datasets. Key hyper-parameters are varied across a
predefined set to balance accuracy and computational efficiency. The following
grid is used for hyper-parameter tuning: the number of importance sampling
operations $R\in\\{5,10,25,50,100\\}$ (100 is omitted for MultiMLPS), nodes
per latent dimension $m_{k}\in\\{1,5,10,20,30\\}$, and standard deviations
$\sigma_{m_{k}}\in\\{0.001,0.01,0.1,1.0\\}$. Early stopping steps (patience)
are set to five or ten, and $R_{\text{infer}}$ is explored at [30, 50, 70].
The option to run the models in parallel must be consistently enabled. Table 5
shows the optimal set of hyper-parameters. This table includes common hyper-
parameters for complete and 50% datasets and another part showing
$\xi_{\mathrm{nodes}}$, which is tuned manually and only relevant to the
missing data.

Table 5: Hyper-parameters for MultiXGBs, MultiETs, and MultiMLPs used for both
complete and missing data.

Dataset | Complete and Missing | Missing |
---|---|---|---
R | $m_{k}$ | $\sigma_{k}$ | Patience | $R_{\mathrm{infer}}$ | $\xi_{\mathrm{nodes}}$ |
MultiXGBs |  |  |  |  |  |  |
space_ga | 10 | 30 | 1.0 | 5 | 70 | 100 |
cpu_activity | 5 | 30 | 1.0 | 5 | 70 | 50 |
naval_propulsion_plant | 5 | 30 | 0.01 | 5 | 50 | 100 |
miami_housing | 5 | 10 | 0.1 | 5 | 50 | 50 |
kin8nm | 5 | 30 | 1.0 | 10 | 70 | 100 |
concrete_compressive_strength | 25 | 30 | 1.0 | 5 | 70 | 100 |
cars | 50 | 10 | 1.0 | 10 | 70 | 100 |
energy_efficiency | 5 | 1 | 0.01 | 10 | 70 | 100 |
california_housing | 5 | 10 | 0.1 | 10 | 50 | 100 |
airfoil_self_noise | 25 | 1 | 0.01 | 10 | 70 | 100 |
QSAR_fish_toxicity | 50 | 30 | 1.0 | 5 | 70 | 100 |
MultiETs |  |  |  |  |  |  |
space_ga | 10 | 30 | 1.0 | 10 | 70 | 100 |
cpu_activity | 5 | 30 | 1.0 | 10 | 70 | 50 |
naval_propulsion_plant | 5 | 30 | 0.01 | 10 | 50 | 100 |
miami_housing | 10 | 10 | 0.1 | 10 | 50 | 50 |
kin8nm | 5 | 30 | 1.0 | 10 | 70 | 100 |
concrete_compressive_strength | 25 | 30 | 1.0 | 10 | 70 | 100 |
cars | 100 | 5 | 0.1 | 10 | 100 | 100 |
energy_efficiency | 5 | 1 | 0.01 | 10 | 70 | 100 |
california_housing | 5 | 10 | 0.1 | 5 | 50 | 100 |
airfoil_self_noise | 25 | 1 | 0.01 | 10 | 70 | 100 |
QSAR_fish_toxicity | 50 | 30 | 1.0 | 10 | 70 | 100 |
MultiMLPs |  |  |  |  |  |  |
space_ga | 25 | 10 | 0.001 | 10 | 50 | 100 |
cpu_activity | 5 | 20 | 0.001 | 5 | 50 | 50 |
naval_propulsion_plant | 5 | 20 | 0.001 | 5 | 50 | 100 |
miami_housing | 5 | 20 | 0.01 | 5 | 50 | 50 |
kin8nm | 5 | 20 | 0.001 | 5 | 50 | 100 |
concrete_compressive_strength | 5 | 30 | 0.001 | 10 | 50 | 100 |
cars | 5 | 30 | 0.1 | 5 | 50 | 100 |
energy_efficiency | 50 | 30 | 0.1 | 10 | 50 | 100 |
california_housing | 5 | 20 | 0.01 | 5 | 50 | 100 |
airfoil_self_noise | 25 | 10 | 0.01 | 10 | 50 | 100 |
QSAR_fish_toxicity | 50 | 30 | 1.0 | 10 | 70 | 100 |

MultiXGBs and MultiMLPs benefit from early stopping to reduce computation time
in complete and incomplete cases. Similarly, the baseline models for these
instances use the same hyper-parameters for early stopping. Further, the
number of epochs in the case of MultiMLPs is set as 1000, except for
energy_efficiency and QSAR_fish_toxicity, where this is changed to 5000. Any
model-specific hyperparameter we did not specify in this paper remains at the
implementation’s default value (e.g., the number of leaves in XGBoost [50]).
Along with the supplementary code, we provide two additional CSV files: one
for the optimal hyperparameters of SEMF models and the other for the results
and hyperparameters of all 330 runs.

Additional conditions are applied only for experiments with missing data.
Datasets california_housing, cpu_activity, miami_housing, and
naval_propulsion_plant—have a patience of five to expedite the training
process. Additionally, for california_housing and cpu_activity, the
$R_{\text{infer}}$ value is set to 30, while for all the other datasets, it is
set to 50. We do this to ensure efficient computation, speed, and memory usage
(especially for the GPU).

For training MultiXGBs and MultiETs, the computations are performed in
parallel using CPU cores (Intel® Core™ i9-13900KF). For MultiMLPs, they are
done on a GPU (NVIDIA® GeForce RTX™ 4090). The GPU is also consistently used
for the missing data simulator and training $p_{\xi}$. All the computations
are done on a machine with 32 GB of memory. The code provides further details
on hardware and reproducibility.

##  Appendix D Metrics for prediction intervals

###  D.1 Common metrics

The most common metrics for evaluating prediction intervals [52, 53] are:

  * •

Prediction Interval Coverage Probability (PICP): This metric assesses the
proportion of times the true value of the target variable falls within the
constructed prediction intervals. For a set of test examples
${(x_{1},y_{1}),\ldots,(x_{N},y_{N})}$, a given level of confidence $\alpha$,
and their corresponding prediction intervals ${I_{1},\ldots,I_{N}}$, the PICP
is calculated as:

$\text{PICP}=\frac{1}{N}\sum_{i=1}^{N}\mathds{1}(y_{i}\in[L_{i},U_{i}]),$ (21)

where $U_{i}$ and $L_{i}$ are the upper and lower bounds of the predicted
values for the $i$-th instance. $y_{i}$ is the actual value of the $i$-th test
example, and $\mathds{1}$ is the indicator function, which equals 1 if $y_{i}$
is in the interval $[L_{i},U_{i}]$ and 0 otherwise. $0\leq\text{{PICP}}\leq 1$
where PICP closer to 1 and higher than the confidence level $\alpha$ is
favored.

  * •

Mean Prediction Interval Width (MPIW): The average width is computed as

$\text{MPIW}=\frac{1}{N}\sum_{i=1}^{N}(U_{i}-L_{i}),$ (22)

which shows the sharpness or uncertainty, where $0\leq\text{{MPIW}}<\infty$
and MPIW close to 0 is preferred.

  * •

Normalized Mean Prediction Interval Width (NMPIW): Since MPIW varies by
dataset, it can be normalized by the range of the target variable

$\text{NMPIW}=\frac{\text{MPIW}}{\max(y)-\min(y)}$ (23)

where $\max(y)$ and $\min(y)$ are the maximum and minimum values of the target
variable, respectively. The interpretation remains the same as MPIW.

###  D.2 Impact of relative metrics for modeling

As our primary focus is on interval prediction, configurations demonstrating
the most significant improvements in $\Delta$CWR and $\Delta$PICP are
prioritized when selecting the optimal hyper-parameters. Furthermore, both
$\Delta$PICP and $\Delta$CWR must be positive, indicating that we must at
least have the same reliability of the baseline (PICP) with better or same
interval ratios (CWR). In instances where no configuration meets the initial
improvement criteria for both metrics, we relax the requirement for positive
$\Delta$PICP to accept values greater than -5% and subsequently -10%, allowing
us to consider configurations where SEMF significantly improves CWR, even if
the PICP improvement is less marked but remains within an acceptable range for
drawing comparisons.

##  Appendix E Full results

Table 6: Test results for MultiXGBs with complete data at 95% quantiles for
seeds 0, 10, 20, 30, 40, with rows ordered by seed (ascending). Performance
over the baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
space_ga | -1% | 6% | -7% | 0.88 | 1.85 | 0.20 | -10% | -13% | 0.60
space_ga | -4% | 5% | -9% | 0.89 | 1.85 | 0.20 | -8% | -10% | 0.59
space_ga | 1% | 9% | -8% | 0.90 | 1.76 | 0.19 | -10% | -10% | 0.62
space_ga | -7% | 5% | -12% | 0.90 | 1.92 | 0.21 | -8% | -10% | 0.59
space_ga | 8% | 3% | 4% | 0.87 | 1.61 | 0.17 | -9% | -8% | 0.60
cpu_activity | 27% | 7% | 16% | 0.88 | 0.37 | 0.07 | 21% | 0% | 0.98
cpu_activity | 33% | 11% | 17% | 0.89 | 0.37 | 0.07 | 23% | -2% | 0.98
cpu_activity | 30% | 10% | 15% | 0.88 | 0.36 | 0.07 | 20% | -5% | 0.98
cpu_activity | 16% | 9% | 7% | 0.89 | 0.38 | 0.07 | 20% | -3% | 0.98
cpu_activity | 18% | 9% | 8% | 0.90 | 0.38 | 0.07 | 22% | 2% | 0.98
naval_propulsion_plant | 176% | 8% | 61% | 0.96 | 0.39 | 0.12 | -20% | -14% | 0.99
naval_propulsion_plant | 151% | 7% | 58% | 0.96 | 0.40 | 0.12 | -9% | -15% | 0.99
naval_propulsion_plant | 136% | 8% | 54% | 0.96 | 0.40 | 0.12 | -15% | -12% | 0.99
naval_propulsion_plant | 156% | 7% | 58% | 0.95 | 0.38 | 0.11 | -12% | -9% | 0.99
naval_propulsion_plant | 161% | 9% | 58% | 0.96 | 0.38 | 0.11 | -15% | -11% | 0.99
miami_housing | 50% | -5% | 37% | 0.82 | 0.46 | 0.06 | 3% | -1% | 0.91
miami_housing | 65% | -7% | 44% | 0.83 | 0.45 | 0.06 | 4% | 1% | 0.91
miami_housing | 62% | -10% | 44% | 0.81 | 0.44 | 0.06 | 7% | 7% | 0.90
miami_housing | 62% | -7% | 43% | 0.82 | 0.44 | 0.05 | 9% | 5% | 0.91
miami_housing | 65% | -9% | 45% | 0.81 | 0.44 | 0.05 | -4% | 0% | 0.90
kin8nm | 12% | 4% | 8% | 0.89 | 1.88 | 0.37 | -17% | -20% | 0.64
kin8nm | 11% | 5% | 6% | 0.89 | 1.91 | 0.37 | -21% | -22% | 0.63
kin8nm | 11% | 5% | 6% | 0.90 | 1.92 | 0.38 | -19% | -21% | 0.64
kin8nm | 5% | 7% | -3% | 0.91 | 2.06 | 0.40 | -20% | -23% | 0.62
kin8nm | 8% | 3% | 5% | 0.88 | 1.91 | 0.37 | -21% | -24% | 0.62
concrete_compressive_strength | 7% | 34% | -26% | 0.91 | 1.26 | 0.26 | -34% | -52% | 0.86
concrete_compressive_strength | 13% | 22% | -7% | 0.89 | 1.25 | 0.25 | -29% | -52% | 0.84
concrete_compressive_strength | 13% | 18% | -5% | 0.91 | 1.20 | 0.24 | -21% | -34% | 0.86
concrete_compressive_strength | 20% | 31% | -9% | 0.93 | 1.27 | 0.26 | -16% | -22% | 0.86
concrete_compressive_strength | 14% | 26% | -10% | 0.90 | 1.25 | 0.25 | -28% | -39% | 0.83
cars | 35% | 14% | 16% | 0.89 | 0.67 | 0.11 | -5% | -1% | 0.95
cars | 16% | 10% | 5% | 0.85 | 0.67 | 0.11 | -4% | -4% | 0.95
cars | 13% | 20% | -6% | 0.91 | 0.78 | 0.13 | -2% | 0% | 0.95
cars | 12% | 16% | -4% | 0.88 | 0.68 | 0.11 | 2% | 5% | 0.96
cars | 16% | 18% | -2% | 0.91 | 0.71 | 0.12 | -7% | -4% | 0.95
energy_efficiency | 255% | 24% | 64% | 0.88 | 0.14 | 0.04 | -50% | -29% | 1.00
energy_efficiency | 320% | 19% | 71% | 0.91 | 0.14 | 0.04 | -3% | -13% | 1.00
energy_efficiency | 257% | 17% | 67% | 0.90 | 0.16 | 0.05 | -21% | -31% | 1.00
energy_efficiency | 322% | 16% | 73% | 0.85 | 0.12 | 0.04 | 17% | 16% | 1.00
energy_efficiency | 78% | 34% | 26% | 0.89 | 0.22 | 0.06 | -17% | -22% | 1.00
california_housing | 24% | 1% | 18% | 0.88 | 1.19 | 0.28 | -2% | -3% | 0.81
california_housing | 28% | 1% | 21% | 0.88 | 1.16 | 0.28 | 0% | -1% | 0.82
california_housing | 33% | 0% | 25% | 0.89 | 1.19 | 0.28 | 0% | 0% | 0.81
california_housing | 30% | 1% | 22% | 0.87 | 1.16 | 0.28 | -1% | -1% | 0.82
california_housing | 38% | -2% | 29% | 0.87 | 1.15 | 0.27 | -1% | -2% | 0.81
airfoil_self_noise | 68% | 6% | 37% | 0.83 | 0.80 | 0.17 | -41% | -43% | 0.89
airfoil_self_noise | 134% | 5% | 55% | 0.81 | 0.59 | 0.13 | -18% | -18% | 0.93
airfoil_self_noise | 13% | -3% | 15% | 0.72 | 1.10 | 0.23 | -109% | -138% | 0.78
airfoil_self_noise | 38% | 20% | 13% | 0.91 | 1.12 | 0.24 | -73% | -76% | 0.85
airfoil_self_noise | 51% | -6% | 38% | 0.76 | 0.84 | 0.18 | -82% | -93% | 0.86
QSAR_fish_toxicity | 11% | 17% | -5% | 0.77 | 1.48 | 0.23 | 3% | 3% | 0.53
QSAR_fish_toxicity | 11% | 9% | 3% | 0.75 | 1.42 | 0.22 | 8% | 7% | 0.58
QSAR_fish_toxicity | 4% | 8% | -5% | 0.77 | 1.48 | 0.23 | 4% | 3% | 0.57
QSAR_fish_toxicity | 10% | 9% | 1% | 0.77 | 1.45 | 0.23 | -2% | -1% | 0.54
QSAR_fish_toxicity | 13% | 15% | -1% | 0.74 | 1.48 | 0.23 | 0% | -6% | 0.52

Table 7: Test results for MultiETs with complete data at 95% quantiles for
seeds 0, 10, 20, 30, 40, with rows ordered by seed (ascending). Performance
over the baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
space_ga | 3% | 0% | 3% | 0.92 | 2.15 | 0.23 | -19% | -24% | 0.51
space_ga | 8% | -1% | 9% | 0.91 | 2.04 | 0.22 | -15% | -18% | 0.54
space_ga | 10% | -1% | 10% | 0.91 | 2.04 | 0.22 | -16% | -17% | 0.54
space_ga | 11% | -2% | 11% | 0.91 | 2.01 | 0.22 | -14% | -17% | 0.55
space_ga | 12% | -3% | 13% | 0.90 | 1.96 | 0.21 | -14% | -15% | 0.55
cpu_activity | 24% | -8% | 26% | 0.91 | 0.48 | 0.09 | -15% | -22% | 0.98
cpu_activity | 24% | -8% | 25% | 0.91 | 0.49 | 0.09 | -17% | -22% | 0.97
cpu_activity | 30% | -9% | 30% | 0.90 | 0.45 | 0.09 | -13% | -18% | 0.98
cpu_activity | 29% | -8% | 28% | 0.90 | 0.47 | 0.09 | -13% | -18% | 0.98
cpu_activity | 27% | -8% | 27% | 0.91 | 0.48 | 0.09 | -14% | -20% | 0.98
naval_propulsion_plant | 265% | -9% | 75% | 0.91 | 0.58 | 0.17 | -268% | -371% | 0.97
naval_propulsion_plant | 205% | -9% | 70% | 0.90 | 0.69 | 0.20 | -332% | -419% | 0.96
naval_propulsion_plant | 178% | -9% | 67% | 0.91 | 0.74 | 0.22 | -382% | -464% | 0.95
naval_propulsion_plant | 198% | -9% | 69% | 0.91 | 0.69 | 0.20 | -336% | -422% | 0.96
naval_propulsion_plant | 248% | -10% | 74% | 0.90 | 0.58 | 0.17 | -283% | -368% | 0.97
miami_housing | 70% | -12% | 48% | 0.86 | 0.66 | 0.08 | -12% | -22% | 0.90
miami_housing | 63% | -12% | 46% | 0.86 | 0.68 | 0.08 | -9% | -21% | 0.90
miami_housing | 68% | -13% | 48% | 0.86 | 0.65 | 0.08 | -9% | -20% | 0.90
miami_housing | 71% | -13% | 49% | 0.86 | 0.63 | 0.08 | -7% | -18% | 0.90
miami_housing | 72% | -12% | 49% | 0.86 | 0.65 | 0.08 | -12% | -21% | 0.89
kin8nm | 12% | -9% | 19% | 0.88 | 2.20 | 0.43 | -32% | -36% | 0.50
kin8nm | 9% | -9% | 17% | 0.88 | 2.28 | 0.45 | -36% | -40% | 0.46
kin8nm | 11% | -10% | 18% | 0.88 | 2.24 | 0.44 | -34% | -37% | 0.49
kin8nm | 9% | -9% | 16% | 0.89 | 2.29 | 0.45 | -34% | -36% | 0.49
kin8nm | 7% | -9% | 15% | 0.88 | 2.33 | 0.46 | -37% | -40% | 0.47
concrete_compressive_strength | 25% | -15% | 32% | 0.79 | 1.02 | 0.21 | -59% | -85% | 0.78
concrete_compressive_strength | 9% | -12% | 19% | 0.81 | 1.19 | 0.24 | -71% | -96% | 0.76
concrete_compressive_strength | 4% | -11% | 14% | 0.83 | 1.22 | 0.25 | -62% | -92% | 0.77
concrete_compressive_strength | -2% | -9% | 8% | 0.83 | 1.31 | 0.27 | -77% | -110% | 0.74
concrete_compressive_strength | 11% | -12% | 20% | 0.81 | 1.14 | 0.23 | -64% | -88% | 0.77
cars | 3% | 15% | -11% | 0.83 | 0.54 | 0.09 | 3% | 2% | 0.95
cars | 16% | 18% | -2% | 0.85 | 0.52 | 0.09 | 3% | 0% | 0.95
cars | 21% | 13% | 7% | 0.80 | 0.47 | 0.08 | 4% | -1% | 0.95
cars | 6% | 17% | -10% | 0.84 | 0.56 | 0.10 | -2% | 0% | 0.94
cars | 10% | 17% | -7% | 0.85 | 0.54 | 0.09 | 8% | 5% | 0.95
energy_efficiency | -3% | -3% | 0% | 0.63 | 0.07 | 0.02 | 3% | 0% | 1.00
energy_efficiency | -3% | 1% | -5% | 0.66 | 0.08 | 0.02 | 7% | 3% | 1.00
energy_efficiency | -14% | 9% | -29% | 0.70 | 0.09 | 0.03 | 0% | 0% | 1.00
energy_efficiency | 2% | -11% | 14% | 0.59 | 0.07 | 0.02 | -3% | -2% | 1.00
energy_efficiency | -5% | 1% | -5% | 0.65 | 0.08 | 0.02 | 6% | 0% | 1.00
california_housing | 33% | -7% | 30% | 0.91 | 1.71 | 0.41 | -14% | -19% | 0.72
california_housing | 28% | -7% | 27% | 0.91 | 1.74 | 0.41 | -16% | -22% | 0.70
california_housing | 34% | -7% | 31% | 0.90 | 1.67 | 0.40 | -14% | -20% | 0.72
california_housing | 29% | -7% | 28% | 0.91 | 1.72 | 0.41 | -16% | -22% | 0.71
california_housing | 31% | -8% | 30% | 0.90 | 1.69 | 0.40 | -15% | -20% | 0.71
airfoil_self_noise | 8% | -21% | 27% | 0.77 | 1.20 | 0.26 | -168% | -189% | 0.74
airfoil_self_noise | -16% | -14% | -3% | 0.84 | 1.68 | 0.36 | -170% | -214% | 0.67
airfoil_self_noise | 20% | -8% | 24% | 0.89 | 1.22 | 0.26 | -95% | -105% | 0.84
airfoil_self_noise | 62% | -18% | 49% | 0.81 | 0.82 | 0.18 | -57% | -81% | 0.89
airfoil_self_noise | 18% | -13% | 26% | 0.86 | 1.19 | 0.25 | -93% | -118% | 0.85
QSAR_fish_toxicity | 11% | -10% | 19% | 0.81 | 1.82 | 0.28 | -6% | -13% | 0.53
QSAR_fish_toxicity | 13% | -10% | 21% | 0.82 | 1.79 | 0.28 | -5% | -11% | 0.54
QSAR_fish_toxicity | 14% | -10% | 21% | 0.82 | 1.82 | 0.28 | -8% | -12% | 0.51
QSAR_fish_toxicity | 11% | -8% | 17% | 0.82 | 1.79 | 0.28 | -5% | -11% | 0.53
QSAR_fish_toxicity | 9% | -8% | 16% | 0.84 | 1.84 | 0.29 | -6% | -10% | 0.52

Table 8: Test results for MultiMLPs with complete data at 95% quantiles for
seeds 0, 10, 20, 30, 40, with rows ordered by seed (ascending). Performance
over the baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
space_ga | 7% | -4% | 10% | 0.79 | 1.22 | 0.13 | -1% | -1% | 0.75
space_ga | 8% | 2% | 5% | 0.82 | 1.25 | 0.14 | 0% | 0% | 0.76
space_ga | 7% | 4% | 3% | 0.83 | 1.27 | 0.14 | -3% | -3% | 0.74
space_ga | 8% | -1% | 8% | 0.79 | 1.24 | 0.14 | 0% | -1% | 0.74
space_ga | 9% | -2% | 10% | 0.81 | 1.27 | 0.14 | 1% | 0% | 0.76
cpu_activity | 8% | -11% | 16% | 0.72 | 0.25 | 0.05 | 6% | 5% | 0.98
cpu_activity | 14% | -15% | 25% | 0.68 | 0.23 | 0.04 | 8% | 5% | 0.98
cpu_activity | 8% | -8% | 15% | 0.73 | 0.24 | 0.04 | 8% | 7% | 0.98
cpu_activity | 5% | -3% | 8% | 0.74 | 0.25 | 0.05 | 6% | 5% | 0.98
cpu_activity | 6% | -8% | 12% | 0.70 | 0.24 | 0.04 | 4% | 2% | 0.98
naval_propulsion_plant | 50% | 0% | 34% | 0.92 | 0.19 | 0.06 | -44% | -48% | 1.00
naval_propulsion_plant | 20% | 0% | 17% | 0.93 | 0.20 | 0.06 | -35% | -13% | 1.00
naval_propulsion_plant | 3% | 1% | 3% | 0.92 | 0.26 | 0.08 | -51% | -45% | 0.99
naval_propulsion_plant | 3% | 3% | 0% | 0.94 | 0.23 | 0.07 | -86% | -63% | 1.00
naval_propulsion_plant | 29% | -6% | 27% | 0.91 | 0.23 | 0.07 | -9% | -11% | 1.00
miami_housing | 13% | -2% | 13% | 0.78 | 0.37 | 0.05 | -2% | 7% | 0.91
miami_housing | 10% | 3% | 6% | 0.81 | 0.40 | 0.05 | -6% | 1% | 0.91
miami_housing | 14% | 1% | 10% | 0.82 | 0.42 | 0.05 | -2% | 3% | 0.92
miami_housing | 8% | 1% | 7% | 0.82 | 0.43 | 0.05 | -1% | 4% | 0.91
miami_housing | 3% | 6% | -4% | 0.83 | 0.44 | 0.06 | -4% | 6% | 0.91
kin8nm | -5% | 9% | -15% | 0.81 | 0.69 | 0.13 | 6% | 5% | 0.93
kin8nm | 15% | 1% | 12% | 0.78 | 0.62 | 0.12 | 7% | 7% | 0.93
kin8nm | 2% | 4% | -2% | 0.82 | 0.67 | 0.13 | 10% | 10% | 0.94
kin8nm | 2% | 14% | -12% | 0.83 | 0.69 | 0.14 | 10% | 6% | 0.94
kin8nm | -3% | 12% | -15% | 0.83 | 0.71 | 0.14 | 4% | 5% | 0.93
concrete_compressive_strength | 89% | -14% | 54% | 0.62 | 0.30 | 0.06 | 4% | 10% | 0.90
concrete_compressive_strength | 82% | -21% | 57% | 0.55 | 0.29 | 0.06 | 9% | 15% | 0.91
concrete_compressive_strength | 120% | -17% | 62% | 0.60 | 0.28 | 0.06 | 18% | 19% | 0.92
concrete_compressive_strength | 100% | -33% | 66% | 0.52 | 0.27 | 0.06 | 9% | 12% | 0.91
concrete_compressive_strength | 101% | -30% | 65% | 0.52 | 0.28 | 0.06 | 18% | 20% | 0.91
cars | 10% | 3% | 7% | 0.77 | 0.38 | 0.06 | -4% | 5% | 0.95
cars | -23% | 28% | -65% | 0.88 | 0.54 | 0.09 | 1% | 1% | 0.96
cars | 16% | -4% | 17% | 0.65 | 0.29 | 0.05 | -6% | -2% | 0.95
cars | 2% | -6% | 8% | 0.68 | 0.34 | 0.06 | -2% | -1% | 0.95
cars | -5% | 8% | -13% | 0.81 | 0.45 | 0.08 | -5% | -2% | 0.95
energy_efficiency | 93% | 59% | 17% | 0.70 | 0.07 | 0.02 | 37% | 38% | 1.00
energy_efficiency | 83% | 27% | 29% | 0.61 | 0.06 | 0.02 | 28% | 32% | 1.00
energy_efficiency | 109% | 28% | 39% | 0.64 | 0.06 | 0.02 | 35% | 32% | 1.00
energy_efficiency | 65% | 26% | 24% | 0.58 | 0.07 | 0.02 | 31% | 31% | 1.00
energy_efficiency | 60% | -15% | 47% | 0.58 | 0.06 | 0.02 | 31% | 32% | 1.00
california_housing | -13% | 9% | -25% | 0.90 | 1.33 | 0.32 | 7% | 8% | 0.81
california_housing | -15% | 6% | -24% | 0.89 | 1.27 | 0.30 | 7% | 9% | 0.82
california_housing | -11% | 9% | -23% | 0.90 | 1.29 | 0.31 | 8% | 11% | 0.82
california_housing | -14% | 8% | -26% | 0.89 | 1.31 | 0.31 | 7% | 8% | 0.81
california_housing | -18% | 9% | -33% | 0.90 | 1.32 | 0.32 | 7% | 9% | 0.82
airfoil_self_noise | 140% | -10% | 62% | 0.74 | 0.34 | 0.07 | 18% | 18% | 0.97
airfoil_self_noise | 42% | 9% | 23% | 0.73 | 0.34 | 0.07 | 11% | 9% | 0.97
airfoil_self_noise | 43% | -1% | 31% | 0.77 | 0.36 | 0.08 | 22% | 22% | 0.97
airfoil_self_noise | 39% | 3% | 26% | 0.73 | 0.35 | 0.07 | 17% | 12% | 0.97
airfoil_self_noise | 105% | -2% | 52% | 0.77 | 0.38 | 0.08 | 21% | 19% | 0.97
QSAR_fish_toxicity | -6% | 8% | -14% | 0.72 | 1.40 | 0.22 | -4% | 6% | 0.52
QSAR_fish_toxicity | 0% | 17% | -17% | 0.82 | 1.50 | 0.23 | 6% | 11% | 0.54
QSAR_fish_toxicity | -5% | 8% | -13% | 0.77 | 1.57 | 0.24 | 0% | -1% | 0.57
QSAR_fish_toxicity | -8% | 13% | -23% | 0.78 | 1.51 | 0.23 | 10% | 9% | 0.61
QSAR_fish_toxicity | -12% | 10% | -25% | 0.73 | 1.49 | 0.23 | 5% | 11% | 0.51

Table 9: Test results for MultiXGBs with 50% missing data at 95% quantiles for
seeds 0, 10, 20, 30, 40, with rows ordered by seed (ascending). Performance
over the baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
space_ga | -12% | 5% | -20% | 0.87 | 2.25 | 0.24 | -11% | -13% | 0.35
space_ga | -9% | 6% | -18% | 0.87 | 2.15 | 0.23 | -12% | -11% | 0.44
space_ga | 3% | 6% | -5% | 0.85 | 1.90 | 0.21 | -1% | -6% | 0.46
space_ga | -8% | 4% | -14% | 0.89 | 2.06 | 0.22 | -5% | -11% | 0.54
space_ga | 0% | 1% | -4% | 0.82 | 1.88 | 0.20 | -6% | -6% | 0.45
cpu_activity | 2% | 2% | -1% | 0.79 | 0.58 | 0.11 | -30% | -30% | 0.73
cpu_activity | 9% | 7% | -12% | 0.87 | 0.61 | 0.12 | -7% | -23% | 0.94
cpu_activity | -41% | 27% | -127% | 0.96 | 1.17 | 0.22 | -1% | -26% | 0.81
cpu_activity | -8% | 14% | -34% | 0.94 | 0.78 | 0.15 | 36% | 5% | 0.88
cpu_activity | 12% | 8% | -12% | 0.91 | 0.66 | 0.12 | -2% | -11% | 0.81
naval_propulsion_plant | 82% | -9% | 47% | 0.76 | 0.64 | 0.19 | -30% | -42% | 0.83
naval_propulsion_plant | 121% | -13% | 59% | 0.76 | 0.65 | 0.19 | 5% | -6% | 0.84
naval_propulsion_plant | 127% | -13% | 58% | 0.73 | 0.66 | 0.20 | 6% | 6% | 0.76
naval_propulsion_plant | 134% | -21% | 63% | 0.64 | 0.52 | 0.15 | -7% | -11% | 0.70
naval_propulsion_plant | 58% | -25% | 53% | 0.62 | 0.64 | 0.19 | -45% | -90% | 0.73
miami_housing | 16% | -6% | 13% | 0.84 | 0.86 | 0.11 | -4% | -4% | 0.72
miami_housing | -9% | 4% | -22% | 0.90 | 0.94 | 0.12 | -14% | -20% | 0.79
miami_housing | 25% | -11% | 17% | 0.72 | 0.60 | 0.07 | -35% | -35% | 0.54
miami_housing | 23% | -4% | 21% | 0.85 | 0.83 | 0.10 | -7% | -12% | 0.80
miami_housing | 25% | -9% | 21% | 0.78 | 0.63 | 0.08 | -10% | -19% | 0.76
kin8nm | -3% | 2% | -12% | 0.88 | 2.41 | 0.47 | 0% | 1% | 0.40
kin8nm | -6% | 3% | -14% | 0.89 | 2.47 | 0.48 | -6% | -8% | 0.39
kin8nm | -4% | 2% | -10% | 0.88 | 2.48 | 0.48 | -5% | -6% | 0.37
kin8nm | -9% | 7% | -23% | 0.90 | 2.57 | 0.50 | -5% | -7% | 0.44
kin8nm | -3% | 3% | -10% | 0.87 | 2.38 | 0.47 | 1% | 2% | 0.41
concrete_compressive_strength | -2% | 0% | -6% | 0.71 | 1.51 | 0.30 | -19% | -22% | 0.48
concrete_compressive_strength | 2% | 11% | -15% | 0.81 | 1.44 | 0.29 | -4% | -3% | 0.69
concrete_compressive_strength | 4% | 11% | -12% | 0.76 | 1.44 | 0.29 | -15% | -24% | 0.50
concrete_compressive_strength | 1% | 8% | -9% | 0.77 | 1.47 | 0.30 | -4% | -9% | 0.64
concrete_compressive_strength | 7% | 8% | -7% | 0.80 | 1.49 | 0.30 | -37% | -44% | 0.57
cars | -3% | 24% | -34% | 0.85 | 1.34 | 0.23 | -9% | -16% | 0.75
cars | -35% | 9% | -98% | 0.82 | 1.78 | 0.30 | -68% | -80% | 0.55
cars | -6% | 29% | -52% | 0.88 | 1.54 | 0.26 | -8% | -14% | 0.76
cars | -29% | 32% | -90% | 0.81 | 1.73 | 0.29 | -28% | -39% | 0.40
cars | -5% | 34% | -40% | 0.92 | 1.49 | 0.25 | -37% | -25% | 0.81
energy_efficiency | 20% | 15% | -13% | 0.93 | 0.71 | 0.21 | -37% | -39% | 0.94
energy_efficiency | 35% | 3% | 23% | 0.88 | 0.58 | 0.17 | -15% | -24% | 0.97
energy_efficiency | 37% | 7% | 19% | 0.92 | 0.63 | 0.18 | -16% | -44% | 0.96
energy_efficiency | -21% | 18% | -81% | 0.95 | 1.01 | 0.29 | -155% | -137% | 0.93
energy_efficiency | -22% | 14% | -66% | 0.91 | 0.92 | 0.27 | -705% | -520% | 0.76
california_housing | 30% | -8% | 28% | 0.82 | 1.52 | 0.36 | -4% | -7% | 0.61
california_housing | 26% | -2% | 21% | 0.89 | 1.61 | 0.39 | 0% | 1% | 0.73
california_housing | 34% | -1% | 26% | 0.90 | 1.66 | 0.40 | -11% | -13% | 0.73
california_housing | 13% | 0% | 11% | 0.90 | 1.85 | 0.44 | -3% | -3% | 0.65
california_housing | 27% | 0% | 21% | 0.91 | 1.70 | 0.41 | -1% | -3% | 0.73
airfoil_self_noise | -8% | 4% | -14% | 0.77 | 1.70 | 0.36 | -14% | -8% | 0.52
airfoil_self_noise | -18% | -2% | -20% | 0.80 | 1.77 | 0.38 | -94% | -110% | 0.50
airfoil_self_noise | -12% | -7% | -11% | 0.70 | 1.86 | 0.40 | -43% | -48% | 0.19
airfoil_self_noise | 34% | -5% | 29% | 0.69 | 1.02 | 0.22 | -19% | -29% | 0.38
airfoil_self_noise | -6% | -6% | -5% | 0.71 | 1.58 | 0.34 | -14% | -22% | 0.43
QSAR_fish_toxicity | -3% | -4% | -1% | 0.68 | 1.56 | 0.24 | -4% | -5% | 0.38
QSAR_fish_toxicity | -1% | -4% | 2% | 0.75 | 1.83 | 0.28 | -9% | -3% | 0.35
QSAR_fish_toxicity | 3% | 9% | -8% | 0.75 | 1.66 | 0.26 | 4% | 4% | 0.43
QSAR_fish_toxicity | -6% | 1% | -19% | 0.72 | 1.61 | 0.25 | -7% | -3% | 0.35
QSAR_fish_toxicity | -1% | -3% | -2% | 0.74 | 1.84 | 0.29 | -5% | -3% | 0.30

Table 10: Test results for MultiETs with 50% missing data at 95% quantiles for
seeds 0, 10, 20, 30, 40, with rows ordered by seed (ascending). Performance
over the baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
space_ga | 7% | -6% | 12% | 0.86 | 2.17 | 0.23 | -17% | -18% | 0.30
space_ga | 3% | -5% | 8% | 0.87 | 2.14 | 0.23 | -14% | -15% | 0.41
space_ga | 2% | -5% | 5% | 0.88 | 2.23 | 0.24 | -9% | -10% | 0.41
space_ga | 7% | -5% | 9% | 0.86 | 2.21 | 0.24 | -7% | -7% | 0.40
space_ga | -1% | -3% | 0% | 0.91 | 2.39 | 0.26 | -14% | -18% | 0.47
cpu_activity | -8% | -5% | -8% | 0.93 | 0.90 | 0.17 | -36% | -46% | 0.93
cpu_activity | -8% | -6% | -4% | 0.92 | 0.89 | 0.17 | 2% | -23% | 0.87
cpu_activity | -30% | -3% | -43% | 0.94 | 1.33 | 0.25 | -36% | -80% | 0.73
cpu_activity | -27% | -3% | -41% | 0.94 | 1.21 | 0.23 | -13% | -49% | 0.74
cpu_activity | -6% | -6% | -3% | 0.91 | 1.00 | 0.19 | -65% | -70% | 0.67
naval_propulsion_plant | 55% | -22% | 48% | 0.77 | 1.29 | 0.38 | -18% | -50% | 0.64
naval_propulsion_plant | 55% | -10% | 42% | 0.89 | 1.52 | 0.45 | -66% | -158% | 0.79
naval_propulsion_plant | 54% | -15% | 44% | 0.84 | 1.44 | 0.42 | -29% | -73% | 0.74
naval_propulsion_plant | 67% | -16% | 50% | 0.83 | 1.29 | 0.38 | -56% | -142% | 0.74
naval_propulsion_plant | 57% | -20% | 48% | 0.78 | 1.38 | 0.41 | -14% | -42% | 0.66
miami_housing | 20% | -6% | 22% | 0.92 | 1.33 | 0.17 | -13% | -26% | 0.78
miami_housing | 32% | -16% | 31% | 0.79 | 0.87 | 0.11 | -41% | -50% | 0.48
miami_housing | -4% | -5% | -1% | 0.92 | 1.39 | 0.17 | -38% | -51% | 0.67
miami_housing | 37% | -14% | 36% | 0.84 | 0.88 | 0.11 | -21% | -34% | 0.79
miami_housing | 16% | -7% | 16% | 0.89 | 1.30 | 0.16 | -20% | -26% | 0.65
kin8nm | 5% | -9% | 13% | 0.86 | 2.48 | 0.49 | -10% | -12% | 0.32
kin8nm | 4% | -8% | 9% | 0.87 | 2.58 | 0.51 | -10% | -12% | 0.33
kin8nm | 5% | -8% | 12% | 0.87 | 2.62 | 0.51 | -13% | -15% | 0.30
kin8nm | 8% | -11% | 17% | 0.84 | 2.41 | 0.47 | -8% | -8% | 0.30
kin8nm | 1% | -5% | 6% | 0.90 | 2.66 | 0.52 | -14% | -18% | 0.37
concrete_compressive_strength | 4% | -11% | 15% | 0.79 | 1.88 | 0.38 | -24% | -36% | 0.46
concrete_compressive_strength | -3% | -17% | 12% | 0.78 | 1.74 | 0.35 | -26% | -36% | 0.42
concrete_compressive_strength | 4% | -12% | 13% | 0.81 | 1.83 | 0.37 | -30% | -47% | 0.56
concrete_compressive_strength | 7% | -11% | 15% | 0.83 | 1.81 | 0.37 | -19% | -31% | 0.58
concrete_compressive_strength | -7% | -18% | 9% | 0.80 | 1.87 | 0.38 | -76% | -99% | 0.45
cars | -18% | 1% | -39% | 0.89 | 1.29 | 0.22 | -16% | -15% | 0.81
cars | -18% | -11% | -11% | 0.84 | 1.39 | 0.23 | -23% | -30% | 0.72
cars | -33% | -6% | -54% | 0.84 | 1.63 | 0.28 | -59% | -59% | 0.60
cars | -28% | -12% | -41% | 0.80 | 1.27 | 0.21 | -18% | -10% | 0.59
cars | -19% | -14% | -29% | 0.79 | 1.48 | 0.25 | -49% | -42% | 0.59
energy_efficiency | -11% | 13% | -35% | 0.91 | 0.50 | 0.14 | -7% | -13% | 0.98
energy_efficiency | 0% | -7% | 0% | 0.83 | 0.49 | 0.14 | -8% | -20% | 0.95
energy_efficiency | -9% | -14% | 0% | 0.76 | 0.47 | 0.14 | -112% | -137% | 0.93
energy_efficiency | -11% | 4% | -30% | 0.87 | 0.59 | 0.17 | -38% | -76% | 0.96
energy_efficiency | -46% | 2% | -114% | 0.89 | 0.52 | 0.15 | -198% | -176% | 0.86
california_housing | 30% | -7% | 28% | 0.92 | 2.09 | 0.50 | -11% | -13% | 0.64
california_housing | 28% | -9% | 28% | 0.87 | 2.02 | 0.48 | -14% | -18% | 0.50
california_housing | 13% | -6% | 17% | 0.92 | 2.24 | 0.53 | -19% | -23% | 0.64
california_housing | 19% | -5% | 20% | 0.93 | 2.23 | 0.53 | -12% | -16% | 0.64
california_housing | 16% | -6% | 18% | 0.93 | 2.35 | 0.56 | -11% | -16% | 0.60
airfoil_self_noise | -7% | -18% | 10% | 0.77 | 2.02 | 0.43 | -27% | -51% | 0.36
airfoil_self_noise | -26% | -18% | -11% | 0.80 | 2.21 | 0.47 | -171% | -229% | 0.30
airfoil_self_noise | 25% | -20% | 36% | 0.74 | 1.28 | 0.27 | -28% | -38% | 0.41
airfoil_self_noise | -22% | -24% | -2% | 0.73 | 2.26 | 0.48 | -54% | -86% | 0.13
airfoil_self_noise | -7% | -18% | 10% | 0.79 | 2.02 | 0.43 | -44% | -58% | 0.43
QSAR_fish_toxicity | 11% | -8% | 15% | 0.85 | 2.37 | 0.37 | -11% | -10% | 0.40
QSAR_fish_toxicity | 13% | -8% | 18% | 0.82 | 2.02 | 0.31 | -1% | -5% | 0.47
QSAR_fish_toxicity | 19% | -7% | 22% | 0.84 | 1.90 | 0.29 | -3% | -3% | 0.48
QSAR_fish_toxicity | -3% | -15% | 12% | 0.76 | 2.02 | 0.32 | -24% | -27% | 0.25
QSAR_fish_toxicity | 8% | -9% | 14% | 0.83 | 2.29 | 0.36 | -3% | -5% | 0.36

Table 11: Test results for MultiMLPs with 50% missing data at 95% quantiles
for seeds 0, 10, 20, 30, 40, with rows ordered by seed (ascending).
Performance over the baseline is highlighted in bold.

| Interval Predictions | Point Predictions
---|---|---
| Relative | Absolute | Relative | Absolute
Dataset |  $\Delta$CWR |  $\Delta$PICP |  $\Delta$NMPIW | PICP | MPIW | NMPIW |  $\Delta$RMSE |  $\Delta$MAE | $\text{R}^{2}$
space_ga | -12% | -16% | 4% | 0.68 | 1.50 | 0.16 | -41% | -40% | 0.15
space_ga | -21% | -17% | -11% | 0.64 | 1.41 | 0.15 | -85% | -61% | -0.21
space_ga | -5% | -14% | 6% | 0.65 | 1.30 | 0.14 | -21% | -20% | 0.37
space_ga | -9% | -13% | 4% | 0.70 | 1.50 | 0.16 | -28% | -25% | 0.35
space_ga | -9% | -13% | 1% | 0.73 | 1.49 | 0.16 | -25% | -26% | 0.38
cpu_activity | -38% | 14% | -89% | 0.87 | 0.75 | 0.14 | 14% | 1% | 0.92
cpu_activity | -60% | 32% | -240% | 0.94 | 1.08 | 0.20 | -20% | -17% | 0.78
cpu_activity | -67% | 36% | -355% | 0.95 | 1.50 | 0.28 | -30% | -34% | 0.72
cpu_activity | -24% | -2% | -28% | 0.68 | 0.39 | 0.07 | -3% | -16% | 0.78
cpu_activity | -63% | 28% | -251% | 0.85 | 1.13 | 0.21 | -37% | -49% | 0.54
naval_propulsion_plant | -41% | 0% | -70% | 0.63 | 0.85 | 0.25 | -62% | -68% | -0.21
naval_propulsion_plant | -53% | -13% | -87% | 0.82 | 1.60 | 0.47 | -141% | -163% | 0.50
naval_propulsion_plant | -48% | -3% | -85% | 0.72 | 1.23 | 0.36 | -99% | -94% | 0.16
naval_propulsion_plant | -88% | -61% | -210% | 0.29 | 1.10 | 0.33 | -215% | -562% | -0.85
naval_propulsion_plant | -33% | 4% | -55% | 0.70 | 1.04 | 0.31 | -24% | -27% | 0.27
miami_housing | -48% | 22% | -148% | 0.93 | 1.33 | 0.17 | -21% | -20% | 0.74
miami_housing | -54% | 27% | -202% | 0.89 | 1.33 | 0.17 | -41% | -48% | 0.53
miami_housing | -19% | 11% | -40% | 0.84 | 0.75 | 0.09 | -18% | -19% | 0.75
miami_housing | -17% | 6% | -32% | 0.82 | 0.69 | 0.09 | -23% | -22% | 0.77
miami_housing | -28% | 1% | -39% | 0.76 | 0.82 | 0.10 | -24% | -32% | 0.57
kin8nm | -7% | -5% | -3% | 0.75 | 1.67 | 0.33 | -12% | -10% | 0.39
kin8nm | -11% | -6% | -8% | 0.77 | 1.80 | 0.35 | -20% | -22% | 0.36
kin8nm | -19% | -4% | -19% | 0.77 | 1.88 | 0.37 | -20% | -19% | 0.31
kin8nm | -4% | -2% | -5% | 0.77 | 1.68 | 0.33 | -3% | 0% | 0.43
kin8nm | -15% | -9% | -8% | 0.78 | 1.79 | 0.35 | -26% | -24% | 0.36
concrete_compressive_strength | -15% | 5% | -26% | 0.74 | 1.53 | 0.31 | -9% | -6% | 0.51
concrete_compressive_strength | -5% | -38% | 34% | 0.40 | 0.57 | 0.12 | -10% | -10% | 0.46
concrete_compressive_strength | -32% | 21% | -80% | 0.84 | 1.69 | 0.34 | -5% | -9% | 0.60
concrete_compressive_strength | 2% | -53% | 54% | 0.27 | 0.39 | 0.08 | -20% | -23% | 0.52
concrete_compressive_strength | -25% | 0% | -33% | 0.76 | 1.37 | 0.28 | -15% | -13% | 0.59
cars | -50% | 38% | -176% | 0.93 | 1.62 | 0.27 | -36% | -38% | 0.75
cars | -44% | 49% | -169% | 0.93 | 1.55 | 0.26 | -3% | -5% | 0.75
cars | -53% | 32% | -232% | 0.94 | 2.04 | 0.34 | -34% | -33% | 0.70
cars | -49% | 44% | -223% | 0.89 | 1.59 | 0.27 | -32% | -15% | 0.48
cars | -57% | 51% | -257% | 0.93 | 1.90 | 0.32 | -32% | -37% | 0.63
energy_efficiency | -35% | 4% | -68% | 0.81 | 0.43 | 0.13 | -8% | 1% | 0.97
energy_efficiency | -12% | 25% | -47% | 0.60 | 0.27 | 0.08 | 8% | 2% | 0.92
energy_efficiency | -19% | -14% | -5% | 0.57 | 0.31 | 0.09 | -1% | -32% | 0.93
energy_efficiency | -18% | 7% | -29% | 0.83 | 0.50 | 0.14 | -8% | 6% | 0.96
energy_efficiency | 20% | 26% | -3% | 0.51 | 0.12 | 0.03 | 14% | 18% | 0.98
california_housing | -16% | 1% | -20% | 0.88 | 1.76 | 0.42 | -7% | -3% | 0.63
california_housing | -14% | 4% | -25% | 0.83 | 1.66 | 0.40 | -4% | -3% | 0.54
california_housing | -9% | 0% | -10% | 0.89 | 1.79 | 0.43 | -1% | 0% | 0.68
california_housing | -20% | 5% | -30% | 0.90 | 1.91 | 0.45 | -2% | -1% | 0.65
california_housing | -27% | 10% | -52% | 0.89 | 2.09 | 0.50 | -8% | -9% | 0.55
airfoil_self_noise | 3% | -1% | 0% | 0.76 | 1.22 | 0.26 | -3% | 1% | 0.54
airfoil_self_noise | 27% | -7% | 26% | 0.77 | 0.84 | 0.18 | -20% | -11% | 0.84
airfoil_self_noise | 8% | -7% | 14% | 0.62 | 0.63 | 0.14 | -22% | -13% | 0.31
airfoil_self_noise | -9% | 1% | -12% | 0.76 | 1.44 | 0.31 | -23% | -20% | 0.44
airfoil_self_noise | -20% | 6% | -33% | 0.83 | 1.62 | 0.35 | -9% | -1% | 0.62
QSAR_fish_toxicity | -17% | -1% | -22% | 0.75 | 1.89 | 0.29 | -21% | -28% | 0.28
QSAR_fish_toxicity | -24% | 9% | -45% | 0.79 | 1.97 | 0.31 | -8% | -9% | 0.44
QSAR_fish_toxicity | -19% | 0% | -23% | 0.70 | 1.56 | 0.24 | -6% | -8% | 0.39
QSAR_fish_toxicity | -17% | 16% | -41% | 0.76 | 1.72 | 0.27 | -21% | -17% | 0.28
QSAR_fish_toxicity | -9% | 8% | -18% | 0.81 | 1.98 | 0.31 | -10% | -9% | 0.32

##  Appendix F Example of learning with non-normality

To illustrate how SEMF adapts to non-normal outcomes, we provide an example
from the naval_propulsion_plant dataset [57]. Figure 2 shows the distribution
of the ground-truth $y$ variable which in this case is
gt_compressor_decay_state_coefficient. The values are uniformly distributed,
and we only standardize the values without changing the shape of the
distribution.

Figure 2: Distribution of the standardized outcome ($y$) variable for the
naval_propulsion_plant dataset which shows that $y$ is uniformly distributed
prior to any training.

After training our SEMF model under the normality assumption with the ideal
hyper-parameters (and a seed of 0), sampling from a normal distribution for
the $z$ dimension, we infer on some randomly sampled test instances that
provide us with the prediction intervals in Figure 3. The ‘SEMF intervals’ can
be compared with XGBoost quantile regression, constituting our baseline. This
figure shows that SEMF’s predicted intervals are better than the baseline.
This plot alone does not tell us much about the predicted output distribution.
Therefore, we provide Figure 4. The last plot shows that for a handful of the
instances, the predicted values can take any shape and are not necessarily
normal.

Figure 3: Predicted intervals and true values on 100 randomly selected test
samples for SEMF MultiXGBs and XGBoost quantile regression with complete data
at 95% quantiles (according to Eq 13).

Figure 4: Violin plot for the first ten instances of the test set where both
the ground truth and 50 values inferred by SEMF (according to Eq 12) are
generated for each instance. We have added jitter here, so the points do not
perfectly align along the x-axis.