1 Introduction
Rare event modeling (REM) deals with situations with an extremely low prevalence of events but influential consequences. Common examples of high social-economic value events include vehicle collisions in traffic safety [10], autonomous vehicle testing [34], life-threatening rare diseases in healthcare, and fraudulent financial activities [32]. Accurate and robust modeling of rare events is critical for the identification of the relevant risk factors so as to predict, and hopefully prevent, associated negative outcomes via intervention. As an example, it is estimated that an adult driver in the U.S. will, on average, travel approximately 6.8 million miles before experiencing a single traffic crash [10]. Given that traffic crashes are relatively rare events compared to normal driving, the precise identification and prediction of such incidents are of utmost importance for reducing fatalities and enhancing overall safety [39].
Characterized by severe event class imbalance and the lack of minority labels, rare event modeling falls outside the comfort zone of standard statistical approaches [46]. Without explicit statistical adjustments, the imbalance drives a learning agent to bias toward the majority class; at the same time, the absence of adequate minority examples causes unprotected models to wrongly capture spurious features that do not generalize. Without adequate generalization ability, the implemented REM model would substantially diminish effectiveness in detecting rare events, which are usually associated with risks that people aim to avoid. In the context of traffic crashes, inaccurate predictions by the REM model could potentially result in severe accidents occurring.
To establish reliability for REMs, the implemented model should have sufficient robustness to handle future out-of-sample predictions for the extreme event, i.e. the data which has not been exposed to during the training of REM. [48] point out that learning from less represented events is vulnerable to data perturbations. These generated misjudgments might lead to safety-critical situations. For instance, if practitioners employ the REM model to assess the safety functions of autonomous vehicles, including decisions about whether to execute maneuvers to prevent traffic crashes based on crash predictions, a less robust REM model could potentially permit self-driving cars to take unforeseen actions, possibly resulting in collisions [5].
Pioneered by the work of [23], there has been extensive research to overcome the potential pitfalls of rare-event modeling in the classical statistical regression framework. Prominent examples include coefficient bias correction [23] and penalized maximum likelihood estimation [14]. However, difficulties arise when applying such correction techniques in real-world settings: (i) the estimation efficiency for rare event modeling is dictated by the number of events of interest given fixed dimensions, and high-dimensional inputs render the conclusions unreliable and ungeneralizable [46, 43]; ($ii$) heavy tails are commonly expected in rare event data, and most existing solutions fail to account for such characteristics [48].
Statistical re-sampling and re-weighting are the two most popular sample adjustment strategies to mitigate severe event class unbalancing issues. Re-sampling approaches alter the learner’s exposure frequency by over-sampling the minority class and down-sampling the majority class. The re-weighting approach adjusts the relative importance by assigning more weight to the less-represented examples. In line with these practices, recent work represented by [4] and [29] also considers modifications to derive class-sensitive loss that properly penalizes the minority event class, allowing the model to learn efficient representation without getting trapped in the complications associated with sample adjustments.
The re-sampling and re-weighting statistical approaches have limited capability to handle the situation when inputs are outside the norm. Re-sampling typically involves either over-sampling or under-sampling [3]. The former is often subjected to the loss of estimation efficiency [46], the latter is often associated with compromised generalization [4]. The re-weighting schemes are often criticized for the numerical instability [35]. Recent works by [4, 29] modify the hinge loss and entropy loss to capitalize the minority class prediction. However, research points out this might result in overfitting issues due to training bias and label noise [35].
An alternative way to overcome the difficulties of rare event modeling is to solicit inductive bias, thereby imposing structural constraints to suppress the spurious associations and identify robust predictive features. Many such models are framed under the few-shot learning setup, capitalizing on its superior ability to generalize in the small sample regime. These models often operate under relatively strong assumptions that the majority event classes sufficiently characterize the feature encoder and meta-predictors, leaving only a handful of parameters to tune for the minority classes, thereby boosting efficiency and generalization. Sharing an aggregating knowledge across different learning tasks defined in majority classes contributes significantly to the success of few-shot learning on generalization [13]. Few-shot learning relies on the assumption that all classes in the data share common features and leverages knowledge transfer from majority classes to enhance learning for minority classes. However, in reality, this assumption may not always hold true, as different classes can have varying data generation processes. For example, traffic crashes and normal driving segments often display distinct characteristics [16]. The failure to meet the assumptions of few-shot learning methods can potentially lead to a degradation in their performance [11].
Adding to the difficulty is the challenge of discerning spurious correlations as rare events do not always associate with stable, recognizable patterns. Instead, feature irregularities are common for rare events, which makes it unattainable to confidently model the rare event behavior. In such scenarios, rare event modeling can be formulated by “establishing the norm” using the majority of examples and to test if input is out of the ordinary, i.e., the anomaly detection approach with abundant tools developed based on various heuristics [8, 37]. Unsupervised anomaly detection approaches often involve density modeling, which has been a long-standing challenge for high-dimensional inputs [8, 37]
Despite the varying degree of empirical success of existing rare event modeling schemes, a few limitations have been largely overlooked. First, assumptions made by different modeling schemes are often at odds, implying performance depends on whether underlying assumptions have been satisfied. Second, representation has not been sufficiently valued in existing research, potentially taking a toll on the minority class generalization. In light of the above, it is imperative to develop robust and interpretable learning strategies that can fully recognize the unique characteristics of rare-event data.
In recognition of the limitations of current rare-event modeling discussed above, we present a novel learning framework, named Variational Inference for Rare Event Modeling (VI-REM), that explicitly addresses some of the weaknesses of existing schemes. Our main contributions include (i) formulation of a variational representation learning scheme based on information theoretic model, allowing disentangled extreme representations for rare-event prediction; ($ii$) development of a robust and interpretable prediction arm that joins the strength of a generalized additive model and an isotonic neural net; ($iii$) presentation of theoretical insights that justify the empirical gains through multiple numerical experiments of VI-REM.
2 Background
Denote the input features as $\mathbf{x}\in {\mathbb{R}^{d}}$, the latent variables as $\mathbf{z}\in {\mathbb{R}^{k}}$, and class label $y\in \{0,1\}$ for which $y=1$ represents the rare-event of interest with probability of occurrence p, and $y=0$ as Normality. Training instances are given by n paired samples $\mathcal{D}={\{{\mathbf{x}^{(i)}},{y^{(i)}}\}_{i=1}^{n}}$ with sample size ${n_{1}}$ in class Event and ${n_{0}}=n-{n_{1}}$ in class Normality. The goal is to predict the conditional likelihood of an event given the features, i.e., $\Pr (y|\mathbf{x})$ with accuracy and robustness. The rare-event problem arises when $p\ll 1$.
2.1 Information Theoretic Generative Model
The generative Bayesian model has gained popularity due to its scalability, interpretability, and generalizability [24]. Under the full Bayesian setup, the generative process for the data likelihood of y given x can be obtained by marginalization over latent variables z:
A common assumption in such generative process is z contains sufficient information of x for y, i.e $\Pr (\text{y}|\mathbf{x},\mathbf{z})=\Pr (\text{y}|\mathbf{z})$. Similar assumption can be found in [41]. With this, Equation 2.1 can be simplified as
For REM which belongs to supervised learning tasks, the core of the generative Bayesian model is to reconstruct the conditional probability of response variable y given predictors x, which is $\Pr (y|\mathbf{x})$, with high accuracy [46, 43].
(2.1)
\[ \Pr (y|\mathbf{x})=\frac{\textstyle\int \Pr (\text{y},\mathbf{x},\mathbf{z})d\mathbf{z}}{\Pr (\mathbf{x})}=\int \Pr (\text{y}|\mathbf{x},\mathbf{z})\Pr (\mathbf{z}|\mathbf{x})d\mathbf{z}\]For REM, maximizing over Equation 2.2 is challenging because of the computational intractability of the integration as well as the difficulty to find the latent z given high-dimensional x.
To address the computational intractability issue for latent z, the Variational Information Bottleneck (VIB) [1] utilizes the information-theoretic methods to find a meaningful representation z from x by maximizing the following equations:
where $I(\mathbf{z},y)$ measures the mutual information between z and y, and the parameter $\beta \gt 0$ is a hyper parameter.
Under the generative Bayesian setup, maximizing Equation 2.4 is identical to minimize the following loss function based on Variational Inference (VI) [24]:
where $\Pr (\mathbf{z})$ denotes a prior distribution on z, ${p_{\theta }}(\mathbf{y}|\mathbf{z})$ denotes a conditional distribution of y given z parameterized by θ, and ${q_{\phi }}(\mathbf{z}|\mathbf{x})$ is a family of distributions parameterized by ϕ to approximate the true posterior density [24].
(2.4)
\[ \begin{aligned}{}{\text{Loss}_{\text{VIB}}}=& -{\mathcal{D}_{\mathrm{KL}}}\Big({q_{\phi }}(\mathbf{z}|\mathbf{x})\| \Pr (\mathbf{z})\Big)\\ {} & +\beta {\mathbb{E}_{{q_{\phi }}(\mathbf{z}|\mathbf{x})}}\Big[\log {p_{\theta }}(\mathbf{y}|\mathbf{z})\Big]\end{aligned}\]The first term ${\mathcal{D}_{\mathrm{KL}}}({q_{\phi }}(\mathbf{z}|\mathbf{x})\| \Pr (\mathbf{z}))$ is the KL divergence between $\Pr (\mathbf{z})$ and ${q_{\phi }}(\mathbf{z}|\mathbf{x})$. It shows the complexity of the posterior distribution. Smaller values mean that ${q_{\phi }}(\mathbf{z}|\mathbf{x})$ is similar to $\Pr (\mathbf{z})$, which implies a simpler posterior. The second term ${\mathbb{E}_{{q_{\phi }}(\mathbf{z}|\mathbf{x})}}[\log {p_{\theta }}(\mathbf{y}|\mathbf{z})]$ shows how well latent z can reconstruct label y, i.e. the task performance.
The parameter β is used to keep a balance between the performance on the task and the complexity of the latent representation z [50]. If $\beta =1$, then Variational Information Bottleneck in Equation 2.4 is the same as maximizing the evidence lower bound (ELBO) of $\Pr (y|\mathbf{x})$, given latent variable z [41]; if $\beta \ne 1$ Variational Information Bottleneck will be the β-Variational Auto Encoder (VAE) [21] in the version of supervised learning.
For supervised learning tasks, maximizing over the VIB objective can improve model’s generalizability [41], because the ${\mathcal{D}_{\mathrm{KL}}}\Big({q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x},\mathbf{y})\| \Pr (\mathbf{z})\Big)$ regularizes the Kullback–Leibler (KL) divergence between approximated posterior ${q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x},\mathbf{y})$ and prior distribution $\Pr (\mathbf{z})$ and helps to prevent ${q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x},\mathbf{y})$ being overfitted because of simplicity in prior distribution. The generative Bayesian model also shows superior robustness over the deterministic models, especially in handling testing data with considerable noises [1].
2.2 Challenges and Motivation
When applying the generative Bayesian model for REM, there are three challenges from statistical perspectives. The first challenge is setting an appropriate prior $\Pr (\mathbf{z})$ for the less represented events. Within the current state of the art (SOTA) approaches, most work uses the isotropic Gaussian [41, 24] or non-parametric based approaches [45]. However, neither of them could be a good fit in REM, where the events of interest are always associated with heavy tail and rarity of observations.
The second challenge in REM comes from the model’s generalization ability to unseen data [48], as the events of interest do not always associate with robust patterns. The disentanglement learning approach proves its effectiveness in representing information with stability, especially for the patterns that are less stable [2]. Promising work includes generalized additive models (GAM) [19] and β-Variational Auto Encoder (VAE) [21]. A feasible approach for disentanglement is to assume sets of statistically independent variables for modeling [7]. In this work, we follow the work in [21] to assume sets of statistically independent variables assuring the model’s disentanglement in REM. This setup helps to decouple the representations into individual features that each capture a unique aspect of data.
The third challenge arises from how we can build an interpretable prediction model. Isotonic regression, also known as monotonic regression, is an important class of constrained regression techniques that works under the assumption of monotonicity [18, 47]. Research has verified its appropriateness in preventing overfitting and enhancing model interpretability [18]. Specifically, for a scalar predictor x, the relationship of x toward $y=g(x)$ is assumed to be non-decreasing with regard to x; i.e., $g({x_{1}})\le g({x_{2}})$ if ${x_{1}}\le {x_{2}}$. Such relation is ubiquitous in practice, and monotonicity allows extra flexibility compared to linear relations yet retains a similar level of interpretability of the predictive features and modeling robustness. In this paper, we leverage the recent introduction of the unconstrained monotonic neural net (UMNN) [47] to approximate flexible mapping functions in the generative Bayesian model.
3 Methodology
We formulate VI-REM under the latent variable model framework, in which the event class labels can be predicted from the latent z. The key idea of the VI-REM is to amortize the difficulty of direct prediction of rare-events to the representation learning stage. The premise is the extreme latent representations lead to extreme events. The approach partitions the representation space into normal and extreme regions, where the prevalence of target events in the latter far exceeds those in the former. This alleviates the majority bias issue that afflicts conventional schemes, as the event distribution is more balanced for the extreme region.
3.1 Extreme Prior Distribution
The Extreme Value Theory (EVT) studies the statistical behavior of extremes [9]. To model the extreme behavior of statistical features, maximum model within the domain of EVT describes the maximal value in a sequence of independent observations. Modeling extremes poses a challenge due to the limited number of observations. However, extreme value theory (EVT) asserts that, under mild conditions [9], extreme values actually belong to the Generalized Extreme Value (GEV) family of distributions. For a sequence of independent and identically distributed random variables, denoted as ${X_{1}},{X_{2}},\dots ,{X_{N}}$, the asymptotic distribution of the maximum value $\max ({X_{1}},{X_{2}},\dots ,{X_{N}})$ follows:
where $\text{GEV}({\mu _{1}},{\sigma _{1}},\xi )$ denotes a GEV distribution with ${\mu _{1}}$ as the location, ${\sigma _{1}}$ as the scale and ξ as the extreme value index. Equation 3.1 is known as the Fisher–Tippett–Gnedenko theorem [9], allowing for accurate modeling of extreme values with very limited examples.
(3.1)
\[ \underset{N\to \infty }{\lim }\Pr \big(\max \left\{{X_{1}},\dots ,{X_{N}}\right\}\le x\big)=\text{GEV}({\mu _{1}},{\sigma _{1}},\xi )\]The GEV has three forms depending on the parameter ξ: Fréchet distribution ($\xi \gt 0$), Gumbel distribution($\xi =0$), and reversed Weibull distribution ($\xi \lt 0$). The Fréchet and Weibull distributions have bounded support, e.g. $x\gt {\mu _{1}}-\frac{{\sigma _{1}}}{\xi }$ for Fréchet distribution. To cover the whole space, we use the Gumbel distribution, denoted as $\mathcal{G}({\mu _{1}},{\sigma _{1}})$, as the prior because the support for Gumbel is $x\in \mathbb{R}$.
In VI-REM, we assume extreme latent in the representation space lead to the happening of Event cases. As these extreme representation usually deviate from the center of normality, we leverage maximum model in Equation 3.1 to represent these extreme cases. We extend standard Gaussian prior [24, 41] via incorporating a Gumbel based prior to accommodate heavy tails. We model the normal regular representations with a Gaussian distribution, and the extreme representation with a Gumbel distribution. We assume the representation space of latent variable of z can be partitioned into the into normal and extreme regions. So that the probability density can be:
where $\mathcal{N}({\mu _{2}},{\sigma _{2}})$ is a Gaussian distribution with parameter ${\mu _{2}}$, ${\sigma _{2}}$. For simplicity, we set $\lambda =\Pr (z\in \textit{Normality})$. We define the extreme prior (EP) as a weighted mixture of a Gumbel distribution and Gaussian distribution in equation 3.2, and its associated probability density function is formulated as:
(3.2)
\[ \begin{aligned}{}\Pr (z)=\Pr (z\in & \textit{Normality})\times \mathcal{N}({\mu _{2}},{\sigma _{2}})\\ {} & +\Pr (z\notin \textit{Normality})\times \mathcal{G}({\mu _{1}},{\sigma _{1}})\end{aligned}\](3.3)
\[ \begin{aligned}{}{p_{\boldsymbol{\kappa }}^{\text{EP}}}(z)=& \frac{\lambda }{\sqrt{2\pi {\sigma _{2}^{2}}}}\exp \Bigg\{-\frac{{(z-{\mu _{2}})^{2}}}{2{\sigma _{2}^{2}}}\Bigg\}\\ {} & +\frac{1-\lambda }{{\sigma _{1}}}\exp \Bigg\{-\frac{z-{\mu _{1}}}{{\sigma _{1}}}-\exp \Big(-\frac{z-{\mu _{1}}}{{\sigma _{1}}}\Big)\Bigg\}\end{aligned}\]Figure 1
Illustration of the Extreme Prior distribution. ${z_{1}}$ and ${z_{2}}$ are simulated independently in the figure. (a) The scatter plot for different EP distributions and the red points come from the Gumbel distribution part and the green points come from the Normal distribution part. For case 1, both of the latent variables ${z_{1}},{z_{2}}\sim \operatorname{EP}(10,2,0.995,0,2)$. (b) We set both the latent factors ${z_{1}},{z_{2}}\sim \operatorname{EP}(10,0.2,0.995,0,2)$ and as ${z_{1}}$, ${z_{2}}$ comes from the same distribution, we only show the results for the probability density estimation for ${z_{1}}$. We consider the EP distribution, kernel density estimation (a non-parametric approach), and Normal distribution to fit the distribution of ${z_{1}}$. The black line represents the estimated entire distribution functions, and the blue displays the distribution function for the tailed part.
To induce more flexibility into the model, we set the parameters $\boldsymbol{\kappa }=(\lambda ,{\mu _{1}},{\sigma _{1}},{\mu _{2}},{\sigma _{2}})$ of EP prior are all learnable during the optimization process. The parameters ${\mu _{1}}$, ${\mu _{2}}$ controls the center of Normal and Gumbel distribution. The parameter ${\sigma _{1}}$ controls the shape of distribution tail. Figure 1a displays the scatter plot for different EP distributions. Theorem 3.1 states that as long as ${\sigma _{1}},{\sigma _{2}}\gt 0$, the proposed EP prior is a proper prior. If the tailed behavior does not actually lead to rarity, the EP can fall back to a standard Gaussian distribution.
Figure 1b compares the performance of different estimators in modeling the tail part of the distributions. The fitted Normal distribution considerably underestimates the tail portion, while non-parametric kernel density estimation tends to overfit the observed samples. In scenarios involving heavy-tailed distributions, which are commonly observed in risk-related analyses [9], traditional parametric methods (such as fitted Normal distribution) and non-parametric methods (like kernel density estimation) may struggle to accurately capture the tail probabilities due to the limited number of rare events. However, extreme value theory-based estimation (EP distribution) excels in handling heavy-tailed scenarios effectively. It not only addresses the challenges posed by limited rare events but also demonstrates substantial generalization ability in accurately describing the tails of distributions.
EP distribution can be generalized to multi-dimensional latent components by applying the above setup to each individual dimension. We assume mutual independence among multi-dimensional latents $\mathbf{z}\in {\mathbb{R}^{k}}$, so as the EP distribution for z with parameter $\boldsymbol{\omega }\in {\mathbb{R}^{5k}}$ can be extended as where ${\mathbf{z}_{i}}\in \mathbb{R}$:
where ${\boldsymbol{\kappa }_{i}}=({\lambda _{i}},{\mu _{1,i}},{\sigma _{1,i}},{\mu _{2,i}},{\sigma _{2,i}})$ and $\boldsymbol{\omega }=({\boldsymbol{\kappa }_{1}},\dots ,{\boldsymbol{\kappa }_{k}})$.
(3.5)
\[ {p_{\boldsymbol{\omega }}^{\text{EP}}}(\mathbf{z})={\prod \limits_{i=1}^{k}}{p_{{\boldsymbol{\kappa }_{i}}}^{\text{EP}}}({\mathbf{z}_{i}})\]The decomposition above assumes that the prior distribution is independent, which implies a possible disentanglement of the latent factor z [21]. One can also consider dependent priors for z. For example, a Gaussian process-based prior $\Pr (\mathbf{z})=\mathcal{N}(\mathbf{0},\mathcal{K}zz)$ with a kernel matrix $\mathcal{K}zz$ [6] to account for dependence in the prior distribution, or to introduce a hierarchical structure to the prior with $\Pr (\mathbf{z})=\Pr ({\mathbf{z}_{1}}){\textstyle\prod _{i=2}^{k}}\Pr ({\mathbf{z}_{i}}\mid {\mathbf{z}_{i-1}})$ [42]. However, more complex priors will increase the number of parameters to be optimized, thus to increase the overfitting risk in REM [45]. We focus on the independent prior in this work, and defer more flexible priors to future research.
3.2 Monotonic Additive Neural Network
To facilitate interpretability and generalization for REM, it is beneficial to impose constraints on ${p_{\boldsymbol{\theta }}}(y|\mathbf{z})$. We combine two prominent techniques: the generalized additive model (GAM) [19] and isotonic regression for ${p_{\boldsymbol{\theta }}}(y|\mathbf{z})$, and we name it as Monotonic Additive Neural Network (MANN). Let y come from a Bernoulli(p) distribution and ${\psi _{\boldsymbol{\theta }}}(y|\mathbf{z})$ be the Logit function and the ${p_{\boldsymbol{\theta }}}(y|\mathbf{z})$ can be the following:
We assume sets of statistically independent variables [7] in MANN to assure model’s disentanglement. Each ${f_{j}}$ takes input from the j-th dimension of the latent to model the generation of y. This additive decomposition encourages disentangled representation. The Logit function ${\psi _{\boldsymbol{\theta }}}(y|\mathbf{z})$ is based on GAM:
where each of the individual function is in Equation 3.8, $j=1,..k$ and k is the dimension of latent variables.
(3.7)
\[ {\psi _{\boldsymbol{\theta }}}(y|\mathbf{z})=\log \left(\frac{{p_{\boldsymbol{\theta }}}(y|\mathbf{z})}{1-{p_{\boldsymbol{\theta }}}(y|\mathbf{z})}\right)=\alpha (\boldsymbol{\theta })+{\sum \limits_{j=1}^{k}}{f_{j}}({\mathbf{z}_{j}})\]It is easy to see all ${f_{j}}$ are monotonic, with the direction dictated by the sign of ${\beta _{j}}$. We model ${h_{j}}(t;\boldsymbol{\theta })$ above using MLP based deep neural networks, and compute ${f_{j}}({\mathbf{z}_{j}})$ via Clenshaw-Curtis numerical integration [47].
We use a toy examples to illustrate the advantages of monotonicity in dealing with REM. We set y from a Bernoulli distribution and the conditional probability follows
where the random variables ${z_{1}}\sim \operatorname{EP}(2,0.5,0.99,-2,0.5)$ and ${z_{2}}\sim \operatorname{EP}(1,0.5,0.99,-1,0.5)$. Through this, the probability for $\Pr (y=1)=0.02$. We set a $2\times 500$ fully connected layers for both MANN and MLP. From Figure 2, the MANN models shows better fittings for the low data density region.
Figure 2
The marginal relationship of ${z_{1}}$, ${z_{2}}$ with the Logit value of condition probability, that is $\log \Big(\frac{\Pr (y|{z_{1}})}{1-\Pr (y|{z_{1}})}\Big)$ and $\log \Big(\frac{\Pr (y|{z_{2}})}{1-\Pr (y|{z_{2}})}\Big)$ for ${z_{1}}$ and ${z_{2}}$. The black points represent the observed samples of the Gumbel part.
3.3 Variational Inference Learning
Our VI-REM model combines the MANN and EP prior distribution. We establish VI-REM on the Multilayer Perceptron (MLP), which is flexible to handle different learning tasks. We adopt the setup from VIB [1], the loss function to be minimized for VI-REM is:
where β is empirically set to be ${n_{0}}/{n_{1}}$ for models to be better optimized. ${p_{\boldsymbol{\omega }}^{\text{EP}}}(\mathbf{z})$ is the EP distribution with parameters $\boldsymbol{\omega }$.
(3.9)
\[ \begin{aligned}{}\mathcal{L}(\boldsymbol{\theta },\boldsymbol{\phi },\boldsymbol{\omega })={\mathcal{D}_{\mathrm{KL}}}& \Big({q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})||{p_{\boldsymbol{\omega }}^{\text{EP}}}(\mathbf{z})\Big)\\ {} -& \beta {\mathbb{E}_{{q_{\phi }}(\mathbf{z}|\mathbf{x})}}\Big(\log {p_{\boldsymbol{\theta }}}(y|\mathbf{z})\Big)\end{aligned}\]The ${q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})$ is a family of distributions to approximate the true posterior density [24]. Practitioners usually refer ${q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})$ as Encoder structure as it reduces the original high-dimensional features x into lower dimensions z. The ${q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})$ can be set as isotropic Gaussian distributions through:
\[ {q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})=\mathcal{N}\Big(\mathbf{z}\mid {\mu _{\phi }}(\mathbf{x}),{\sigma _{\phi }^{2}}(\mathbf{x})\times \mathbf{I}\Big)\]
where I is the identity matrix. The parameters ${\mu _{\phi }}(\mathbf{x}),{\sigma _{\phi }^{2}}(\mathbf{x})={g_{\boldsymbol{\phi }}}(\mathbf{x})$ and ${g_{\phi }}(\mathbf{x})$ is established based on MLP parameterized by $\boldsymbol{\phi }$. This is a standard choice in literature of VI [25]. The usage of isotropic Gaussian distributions imply the individual components of the posterior ${q_{\phi }}({\mathbf{z}_{1}}|\mathbf{x})$, ..., ${q_{\phi }}({\mathbf{z}_{K}}|\mathbf{x})$ are mutually independent [25]. For more complex cases with dependent posterior, we can use distributions with tractable likelihood, such as normalizing flows [36] or inverse autoregressive flows [26]. Our results show that isotropic Gaussian distributions perform well. We use them in the following analysis. $\log {p_{\boldsymbol{\theta }}}(y|\mathbf{z})$ is based on MANN parameterized by $\boldsymbol{\theta }$:
Updating the parameters $\boldsymbol{\theta }$, $\boldsymbol{\phi }$, $\boldsymbol{\omega }$ in VI-REM is based on stochastic gradient descent. Updating the parameter $\boldsymbol{\phi }$ is the same as the strategies in [41, 24]. As for parameters $\boldsymbol{\omega }\in {\mathbb{R}^{5k}}$ for ${p_{\boldsymbol{\omega }}^{\text{EP}}}(\mathbf{z})$ and we calculate the derivatives with Equation 3.5:
The parameters ${\lambda _{1}},\dots ,{\lambda _{k}}\in (0,1)$ determine the weight for the Gumbel distribution within the EP prior. Their derivatives can be calculated as follows:
where $i=1,\dots ,k$, and k is the dimension of the latent variable z; ${\mathbf{z}_{i}}$ is the ith component of latent z; ${\lambda _{i}}$ is the parameter λ for latent ${\mathbf{z}_{i}}$; ${h_{1,i}}=\frac{1}{\sqrt{2\pi {\sigma _{2,i}^{2}}}}\exp \Big(\frac{-{({\mathbf{z}_{i}}-{\mu _{2,i}})^{2}}}{2{\sigma _{2,i}^{2}}}\Big)$; and ${h_{2,i}}=\frac{1}{{\sigma _{1,i}}}\exp \Big(\frac{-({\mathbf{z}_{i}}-{\mu _{1,i}})}{{\sigma _{1,i}}}-\exp \big(\frac{-({\mathbf{z}_{i}}-{\mu _{1,i}})}{{\sigma _{1,i}}}\big)\Big)$.
(3.11)
\[ \begin{aligned}{}{\nabla _{\boldsymbol{\omega }}}\mathcal{L}(\boldsymbol{\theta },\boldsymbol{\phi },\boldsymbol{\omega })& ={\nabla _{\boldsymbol{\omega }}}{\mathcal{D}_{\mathrm{KL}}}[{q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})||{p_{\boldsymbol{\omega }}}(\mathbf{z})]\\ {} & =-{\mathbb{E}_{{q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})}}\Big({\sum \limits_{i=1}^{k}}{\nabla _{\boldsymbol{\omega }}}\log {p_{{\boldsymbol{\kappa }_{i}}}}({\mathbf{z}_{i}})\Big)\end{aligned}\](3.12)
\[ \begin{aligned}{}{\nabla _{{\lambda _{i}}}}\mathcal{L}(\boldsymbol{\theta },\boldsymbol{\phi },\boldsymbol{\omega })& =-{\mathbb{E}_{{q_{\boldsymbol{\phi }}}({\mathbf{z}_{i}}|\mathbf{x})}}\Big(\frac{{\nabla _{{\lambda _{i}}}}{p_{{\boldsymbol{\kappa }_{i}}}}({\mathbf{z}_{i}})}{{p_{{\boldsymbol{\kappa }_{i}}}}({\mathbf{z}_{i}})}\Big)\\ {} & =-{\mathbb{E}_{{q_{\boldsymbol{\phi }}}({\mathbf{z}_{i}}|\mathbf{x})}}\Big\{\frac{{h_{1,i}}-{h_{2,i}}}{{\lambda _{i}}{h_{1,i}}+(1-{\lambda _{i}}){h_{2,i}}}\Big\}\end{aligned}\]During the training process, parameter ${\lambda _{i}}$ might collapse to zero due to sparse parameter regime. To avoid this issues, we set parameter ${\lambda _{i}}$ through a logistic function as ${\lambda _{i}}=\frac{1}{1+\exp (-{\Lambda _{i}})}$. This transformation enables us to update the weight for ${\lambda _{i}}\in (0,1)$ by updating the real-valued parameter ${\Lambda _{i}}\in \mathbb{R}$. The logistic function helps prevent the parameter ${\lambda _{i}}$ from collapsing to zero. To obtain ${\lambda _{i}}$, we utilize the chain rule on Equation 3.12 to compute the derivatives of ${\nabla _{{\Lambda _{i}}}}\mathcal{L}(\boldsymbol{\theta },\boldsymbol{\phi },\boldsymbol{\omega })$. The effectiveness of this approach is confirmed through numerical studies in Sections 5 and 6.
The calculation of derivatives $\boldsymbol{\theta }$ involves the derivatives of each individual function ${f_{j}}({\mathbf{z}_{j}})$ in Equation 3.8. Based on Leibniz integral rule, we can have:
Based on Equation 3.13, the derivatives of the Logit function ${\psi _{\boldsymbol{\theta }}}(y|\mathbf{z})$ with respect to $\boldsymbol{\theta }$ can be obtained through:
(3.14)
\[ \begin{aligned}{}{\nabla _{\boldsymbol{\theta }}}{\psi _{\boldsymbol{\theta }}}(y|\mathbf{z})={\nabla _{\boldsymbol{\theta }}}& \alpha (\boldsymbol{\theta })+{\sum \limits_{j=1}^{k}}{\nabla _{\boldsymbol{\theta }}}{\beta _{j}}(\boldsymbol{\theta })\\ {} & +{\sum \limits_{j=1}^{k}}{\int _{{s_{0}}}^{{\mathbf{z}_{j}}}}{\nabla _{\boldsymbol{\theta }}}\exp \Big({h_{j}}(t;\boldsymbol{\theta })\Big)\mathrm{d}t\end{aligned}\]With $\mathbf{z}\sim {q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})$, the derivatives w.r.t parameters $\boldsymbol{\theta }$ are:
With the optimized parameters $\boldsymbol{\theta }$, $\boldsymbol{\phi }$ through the introduced learning process, we can generate the predictions for $\Pr (y|\mathbf{x})$ through a Monte Carlo estimation where L is number of generated samples, where ${\mathbf{z}^{(l)}}\sim {q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})$:
(3.15)
\[ \begin{aligned}{}& {\nabla _{\boldsymbol{\theta }}}\mathcal{L}(\boldsymbol{\theta },\boldsymbol{\phi },\boldsymbol{\omega })={\mathbb{E}_{{q_{\phi }}(\mathbf{z}|\mathbf{x})}}\Big[{\nabla _{\boldsymbol{\theta }}}\log {p_{\theta }}(y|\mathbf{z})\Big]\\ {} & ={\mathbb{E}_{{q_{\phi }}(\mathbf{z}|\mathbf{x})}}\Big[\frac{\exp (-{\psi _{\theta }}(y|\mathbf{z}))}{1+\exp (-{\psi _{\theta }}(y|\mathbf{z}))}\times {\nabla _{\boldsymbol{\theta }}}{\psi _{\boldsymbol{\theta }}}(y|\mathbf{z})\Big]\end{aligned}\]3.4 Theoretical Justifications of VI-REM
We theoretically discuss modeling the tail in REM can improve generalization power. The generalization power of $\Pr (y|\mathbf{x})$ involve ${q_{\phi }}(\mathbf{z}|\mathbf{x})$ and ${p_{\theta }}(y|\mathbf{z})$ as in Equation 3.16. The asymptotic behavior of ${q_{\phi }}(\mathbf{z}|\mathbf{x})$ through VI has been extensively discussed in [49]. This paper focuses on ${p_{\theta }}(y|\mathbf{z})$.
We use generalization gap to measure model’s generalization ability, which is defined as the difference between the model’s performance on training data and on unseen testing data. Denote ${\mathrm{P}_{n}}\big(g(\mathbf{z})\big)$ and $\Pr \big(g(\mathbf{z})\big)$ respectively as the empirical measure and underlying ground truth for some function $g(\mathbf{z})$. Define $\mathcal{F}$ as the set of all monotonic functions, let f be the candidate function belonging to function space $\mathcal{F}$, and ${f_{0}}\in \mathcal{F}$ is the ground truth. We use ${\mathrm{L}_{f}}(\mathbf{z},y)$ to denote the loss measured by function f taking value at z for the response y. We identify ${\Omega _{\eta }}\triangleq \{\| \mathbf{z}\| \lt \eta \}$ as an η-ball where the bulk of the distribution mass resides, i.e., $\Pr ({\Omega _{\eta }})\ge 1-\epsilon $, $\epsilon \ll 1$, and consequently ${\Omega _{\eta }^{C}}\triangleq {\mathbb{R}^{d}}\setminus \Omega $ contains the extremes.
We point out the theoretical generalization gap with the introduction of cutoff points in Theorem 3.2. The following regularized assumptions are needed to establish Theorem 3.2:
Theorem 3.2.
Under the assumptions in C.1–C.3 described above, the following generalization gap holds
where $\delta =|\mathcal{F}|{e^{-2n{\phi ^{2}}}}$ and $|\mathcal{F}|$ is the cardinality of space $\mathcal{F}$. ϕ, γ are some constant, ${I_{\mathbf{z}\in {\Omega _{\eta }}:||\mathbf{z}||\le \eta }}$ is parameter λ in EP distribution, and h is the empirical value of ${I_{\mathbf{z}\in {\Omega _{\eta }}:||\mathbf{z}||\le \eta }}$.
(3.17)
\[ \begin{aligned}{}\mathbb{E}\Big(& \underset{f\in \mathcal{F}}{\sup }\big|{\mathrm{P}_{n}}\big({\mathrm{L}_{f}}(\mathbf{z},y)\big)-\Pr \big({\mathrm{L}_{{f_{0}}}}(\mathbf{z},y)\big)\big|\Big)\\ {} & \le \sqrt{\frac{\log |\mathcal{F}|+\log (1/\delta )}{2n}}+\gamma {h^{1/2}}{\eta ^{1/2}}\underset{f\in \mathcal{F}}{\sup }\| f-{f_{0}}{\| _{\mathbf{z}\in {\Omega _{\eta }}}^{1/2}}\\ {} & +\gamma \underset{\mathbf{z}\in {\Omega _{\eta }^{C}}}{\sup }\| \mathbf{z}{\| ^{1/2}}{(1-h)^{1/2}}\underset{f\in \mathcal{F}}{\sup }\| f-{f_{0}}{\| _{\mathbf{z}\in {\Omega _{\eta }^{C}}}^{1/2}}\end{aligned}\]Through Theorem 3.2, we decompose the generalization gap into estimation error, the approximation error for the non-tailed part, and the approximation error for the tailed part over a certain cutoff point as stated in the following Remark.
Remark 1.
Assume ${f_{1}},{f_{2}}\in \mathcal{F}$ are two functions belonging the space $\mathcal{F}$. We set ${f_{1}}$ as the $\log {p_{\boldsymbol{\theta }}}(y|\mathbf{z})$ which is obtained by maximizing Equation 3.9. Similarly, we set ${f_{2}}$ as $\log {p_{{\boldsymbol{\theta }^{\prime }}}}(y|{\mathbf{z}^{\prime }})$is obtained by optimizing
where r is some constant. Then based on conditions C.1–C.3 and Theorem 3.2, the generalization gap over function ${f_{1}}$, ${f_{2}}$ can have the following:
\[ \beta {\mathbb{E}_{{q_{{\boldsymbol{\phi }^{\prime }}}}({\mathbf{z}^{\prime }}|\mathbf{x})}}\Big(\log {p_{{\boldsymbol{\theta }^{\prime }}}}(y|{\mathbf{z}^{\prime }})\Big)-{\mathcal{D}_{\mathrm{KL}}}\Big({q_{{\boldsymbol{\phi }^{\prime }}}}({\mathbf{z}^{\prime }}|\mathbf{x})\mid {p^{\textit{Gauss}}}({\mathbf{z}^{\prime }})\Big)\]
where ${p^{\textit{Gauss}}}(\mathbf{z})$ is a isotropic Gaussian distribution. Denote ${\Delta _{{\Omega _{\eta }^{C}}}}=||{f_{2}}-{f_{0}}|{|_{\mathbf{z}\in {\Omega _{\eta }^{C}}}^{1/2}}-||{f_{1}}-{f_{0}}|{|_{\mathbf{z}\in {\Omega _{\eta }^{C}}}^{1/2}}$ and ${\Delta _{{\Omega _{\eta }}}}=||{f_{1}}-{f_{0}}|{|_{\mathbf{z}\in {\Omega _{\eta }}}^{1/2}}-||{f_{2}}-{f_{0}}|{|_{\mathbf{z}\in {\Omega _{\eta }}}^{1/2}}$. As long as the follows holds:
(3.18)
\[ {\Delta _{{\Omega _{\eta }^{C}}}}/{\Delta _{{\Omega _{\eta }}}}\gt \frac{{\eta ^{1/2}}}{r{\sup _{\mathbf{z}\in {\Omega _{\eta }^{C}}}}\| \mathbf{z}{\| ^{1/2}}}{\Big(\frac{h}{1-h}\Big)^{1/2}}\](3.19)
\[ \begin{aligned}{}\mathbb{E}\Big(& \underset{{f_{1}}\in \mathcal{F}}{\sup }\big|{\mathrm{P}_{n}}\big({\mathrm{L}_{{f_{1}}}}(\mathbf{z},y)\big)-\Pr \big({\mathrm{L}_{{f_{0}}}}(\mathbf{z},y)\big)\big|\Big)\\ {} & \lt \mathbb{E}\Big(\underset{{f_{2}}\in \mathcal{F}}{\sup }\big|{\mathrm{P}_{n}}\big({\mathrm{L}_{{f_{2}}}}({\mathbf{z}^{\prime }},y)\big)-\Pr \big({\mathrm{L}_{{f_{0}}}}({\mathbf{z}^{\prime }},y)\big)\big|\Big)\end{aligned}\]Remark 1 supports that the proposed extreme prior can improve model generalization. The misrepresented tailed part in REM amplifies the approximation error, thus restricting the model’s generalizability. Intuitively from Figure 1b, utilizing the heavy-tailed prior distribution remedies the misrepresentation problem, thus helps to tighten the generalization bound.
4 Model Implementations
4.1 Model Interpretation
The interpretation of VI-REM, i.e., the relationship between x and label y, needs to be made through latent factor z. The MANN component in recognizing the ${p_{\boldsymbol{\theta }}}(y|\mathbf{z})$ makes the relationship between latent factors z and response y easy to assess. However, interpreting the relationship between x and z is challenging due to the black-box neural network. We frame the model interpretability as a variable selection process, i.e., identifying the most relevant x for each of the latent factors.
The core of this interpretation approach is the evaluation of the statistical dependence between x and z. Recent studies suggest that distance correlation [44] is a powerful non-parametric tool to measure the statistical correlation between two random variables. The MANN structure guarantees that the latent factors z in VI-REM are largely statistically independent. Based on this fact, we can use the Distance Correlation [44] to select the critical x in determining ${\mathbf{z}_{1}}$, $...{\mathbf{z}_{k}}$.
Considering the high dimensionality of the input features x, we adopt the Distance Correlation Sure Independence Screening (DC-SIS) framework [27], which is the SOTA approach to use Distance Correlation to select high-dimensional features. The essence of a DC-SIS is to rank variables based on the Distance Correlation. Under our setup, we can use the Distance Correlation between ${\mathbf{z}_{i}}$ and each of the input features ${\mathbf{x}_{1}}$, ..., ${\mathbf{x}_{d}}$ to rank the importance of input features and select the most relevant ones. However, when directly using the DC-SIS to the proposed VI-REM model, the $\mathcal{O}({n^{2}})$ computational complexity for Distance Correlation is intolerable for a large sample size [28], especially when hundreds of thousands of training data. Therefore, we adopted a Divide and Conquer algorithm [28] to overcome the computation challenge.
For latent factor ${\mathbf{z}_{j}}$, we equally spilt the n samples into G groups with $s=n/G$ data points in each subgroup. For each subgroup $g\in [1,\dots G]$ with s samples, we measure the Distance Correlation of ${\mathbf{x}_{l}^{g}}$ and ${\mathbf{z}_{j}^{g}}$, where ${\mathbf{x}_{l}^{g}}$ represents the input feature ${\mathbf{x}_{l}}$ within subgroup g and ${\mathbf{z}_{j}^{g}}$ represents the latent factor ${\mathbf{z}_{j}}$ within subgroup g. For simplicity, denote $\mathbf{u}={\mathbf{x}_{l}^{g}}$ and $\mathbf{v}={\mathbf{z}_{j}^{g}}$, ${d_{u}}$, ${d_{v}}$ as the dimension for u and v, and ${\mathbf{u}^{(i)}}$, ${\mathbf{v}^{(i)}}$ as the i th sample of u and v correspondingly. The detailed formula for estimation of distance correlation is as follows:
where
and ${\widehat{S}_{1}^{g}}$, ${\widehat{S}_{2}^{g}}$, ${\widehat{S}_{3}^{g}}$ for ${\widehat{\operatorname{dcov}}^{2}}(\mathbf{u},\mathbf{v})$ are as follow:
(4.1)
\[ {\widehat{\operatorname{DC}}_{l}^{g}}=\widehat{\operatorname{dcorr}}({\mathbf{x}_{l}^{g}},{\mathbf{z}_{j}^{g}})=\frac{\widehat{\operatorname{dcov}}(\mathbf{u},\mathbf{v})}{\sqrt{\widehat{\operatorname{dcov}}(\mathbf{u},\mathbf{u})\widehat{\operatorname{dcov}}(\mathbf{v},\mathbf{v})}}\](4.2)
\[ {\widehat{\operatorname{dcov}}^{2}}(\mathbf{u},\mathbf{v})={\widehat{S}_{1}^{g}}+{\widehat{S}_{2}^{g}}-2{\widehat{S}_{3}^{g}}\](4.3)
\[ \begin{array}{l}\displaystyle {\widehat{S}_{1}^{g}}=\displaystyle \frac{1}{{s^{2}}}{\displaystyle \sum \limits_{i=1}^{s}}{\displaystyle \sum \limits_{j=1}^{s}}{\left\| {\mathbf{u}^{(i)}}-{\mathbf{u}^{(j)}}\right\| _{{d_{u}}}}{\left\| {\mathbf{v}^{(i)}}-{\mathbf{v}^{(j}}\right\| _{{d_{v}}}}\\ {} \displaystyle {\widehat{S}_{3}^{g}}=\displaystyle \frac{1}{{s^{3}}}{\displaystyle \sum \limits_{i=1}^{s}}{\displaystyle \sum \limits_{j=1}^{s}}{\displaystyle \sum \limits_{l=1}^{s}}{\left\| {\mathbf{u}^{(i)}}-{\mathbf{u}^{(l)}}\right\| _{{d_{u}}}}{\left\| {\mathbf{v}^{(j)}}-{\mathbf{v}^{(l)}}\right\| _{{d_{v}}}}\\ {} \displaystyle {\widehat{S}_{2}^{g}}=\displaystyle \frac{1}{{s^{2}}}{\displaystyle \sum \limits_{i=1}^{s}}{\displaystyle \sum \limits_{j=1}^{s}}{\left\| {\mathbf{u}^{(i)}}-{\mathbf{u}^{(j)}}\right\| _{{d_{u}}}}\displaystyle \frac{1}{{s^{2}}}{\displaystyle \sum \limits_{i=1}^{s}}{\displaystyle \sum \limits_{j=1}^{s}}{\left\| {\mathbf{v}^{(i)}}-{\mathbf{v}^{(j)}}\right\| _{{d_{v}}}}\end{array}\]For each subgroup, we conduct the DC-SIS [27] to select the set with the most contribution of input features x toward the latent factor ${\mathbf{z}_{j}}$. In practice for each of the subgroup g, the selection process for a given size $q\lt s=n/G$ can be:
where $l\in [1,\dots ,d]$ and d is the dimension of input features x.
(4.4)
\[ {\widehat{\mathcal{D}}_{g}}=\left\{l:{\widehat{\operatorname{DC}}_{l}^{g}}\hspace{2.5pt}\text{is among top}\hspace{2.5pt}q\hspace{2.5pt}\text{largest of}\hspace{2.5pt}d\hspace{2.5pt}\right\},\]After conducting these G (in total we have G subgroups) independent calculations and selections, the final component of ${\mathbf{z}_{j}}$ denoted by ${\widehat{\mathcal{D}}_{{\mathbf{z}_{j}}}}$ is established by: ${\widehat{\mathcal{D}}_{{\mathbf{z}_{j}}}}={\cap _{g\in \{1,\dots G\}}}{\widehat{\mathcal{D}}_{g}}$. In practice, we can simply uses the most common selected relevant features ${\mathbf{x}_{1}}$, ..., ${\mathbf{x}_{d}}$ in each of the subgroup procedures. Repeat the above process for all ${\mathbf{z}_{j}}$ and calculate the ${\widehat{\mathcal{D}}_{{\mathbf{z}_{j}}}}$; i.e., the identified component of latent factors.
Theorem 4.1.
Denote ${\mathcal{D}_{\mathbf{{z_{j}}}}^{\ast }}$ as the true components associated with latent $\mathbf{{z_{j}}}$, and ${\widehat{\mathcal{D}}_{\mathbf{{z_{j}}}}}$ as the set identified by the previously mentioned model’s interpretation strategy. Then:
where $G={n^{1-\alpha }}$, $\alpha \gt 0$, d is the dimension of input features x, ${c_{1}}$, ${c_{2}}$, u, v are some constant and $u+v\lt 1/2$.
(4.5)
\[ \begin{aligned}{}\Pr \left({\mathcal{D}_{\mathbf{{z_{j}}}}^{\ast }}\subseteq {\widehat{\mathcal{D}}_{\mathbf{{z_{j}}}}}\right)\ge \Big\{1& -\mathcal{O}\{d\exp [-{c_{1}}{n^{\alpha -2\alpha (u+v)}}]\\ {} & +d{n^{\alpha }}\exp (-{c_{2}}{n^{v}})\}{\Big\}^{G}}\to 1\end{aligned}\]Theorem 4.1 shows that as long as the sample size n is getting sufficiently large, the true risk modulators for the learnt latents have high confidence to be selected by the given procedure latent factor identification.
4.2 Robustness Evaluations
A robust prediction model should be less vulnerable to attacks and remain valid when testing data is out of domain [5]. Building a robust model is critical for deploying deep learning models in the application of safety-critical situation [48]. Robustness measures the model’s ability to handle data perturbations and helps understand the boundary of successful applications. Mathematically for a model f and its predictions over x $f(x)$, robustness evaluations aims to find how the model f can handle $f({x^{\prime }})$, where ${x^{\prime }}$ has minor perturbations to x. if $f({x^{\prime }})$ does not show much difference with $f(x)$, then model f is robust and trustworthy to deploy [48].
Researchers have proposed various approaches under different application scenarios to generate ${x^{\prime }}$ to evaluate robustness. Regarding natural language processing tasks, [51] summarizes the current well-established methods and suggests mixing the words from different documents [22, 51] could be a simple and efficient approach. For computer vision studies, image blurring and image rotation are all useful in practice. In this work, we use the rotation operator to investigate whether the model can efficiently extract invariant features from rare events.
4.3 Extension to General Setups
The previous VI-REM model mostly focuses on binary classification, and we generalize it to multi-label classification in risk analysis. Suppose we have $\mathrm{L}+1$ classified labels, and we treat one dominant class as class Normality, the rest L minority classes as Event. For example, traffic crashes are rare and can be classified from the most severe fatal crashes to injury crashes to minor property-damage only crashes. To mimic the general settings in risk analysis, we assume the severity of the L minority class ranges from low to high.
Without generality, we set 1, ..., L as labels for the L minority classes, and class 1 as the least severe Event and class L is the most severe one. As the response variable, the severity of an event increases from class 1 toward class L and belongs to the ordinal variables, we applied the setup from the Ordered Logit Model (OLM). Under the OLM setup, we want to model $\Pr (y\ge m|\mathbf{x})$ given different event severity where $m=1,,,.L$. So that under the VI-REM, we modify the Logit function in Equation 3.6 through:
and the Logit function in the OLM based VI-REM follows:
where the coefficients ${\boldsymbol{\alpha }_{\boldsymbol{m}}}(\boldsymbol{\theta })$ are intercepts for different event categories and $0\lt {\boldsymbol{\alpha }_{\mathbf{1}}}(\boldsymbol{\theta })\lt \cdots \lt {\boldsymbol{\alpha }_{\boldsymbol{L}}}(\boldsymbol{\theta })$. The individual function ${f_{j}}({\mathbf{z}_{j}})$ function follows the same as Equation 3.8.
(4.6)
\[ {\psi _{\boldsymbol{\theta }}}\big(y\ge m|\mathbf{z}\big)=\log \left(\frac{{p_{\boldsymbol{\theta }}}(y\ge m|\mathbf{z})}{1-{p_{\boldsymbol{\theta }}}(y\ge m|\mathbf{z})}\right)\](4.7)
\[ {\psi _{\boldsymbol{\theta }}}\big(y\ge m|\mathbf{z}\big)={\boldsymbol{\alpha }_{\boldsymbol{m}}}(\boldsymbol{\theta })+{\sum \limits_{j=1}^{k}}{f_{j}}({\mathbf{z}_{j}})\]4.4 Experiment Setup
We considered the following machine learning approaches for comparison: (1) Sampling MLP, short as S-MLP, is a re-sampling scheme based on Multilayer Perceptron (MLP) over-samples from the less-represented class [3]. The hyperparameter r in the S-MLP is used to control the over-sample ratio, which is equal to the expected ${n_{1}}$ over ${n_{0}}$ after over-sampling. We use grid search cross-validation to select the sampling ratio $r\in (0.01,0.05,0.1,0.2)$. (2) Focal MLP, short as F-MLP, uses a re-weighting based loss that adapts the cross entropy loss putting more penalization on the less-represented class [29] based on MLPs. We use grid search cross-validation to select the hyperparameters $\alpha \in (0.25,0.5,0.75)$, $\gamma \in (2,5)$. (3) MAML is a few-shot learning scheme from [13] and it focuses on improving the generalization performance for less represented data samples and we use 5-shot learning.
We also consider the Deep-SVDD [37], short as D-SVDD, which is a unsupervised anomaly detection approach. D-SVDD will train a MLP based neural network to represent the majority of the data. The data points falling outside the majority will be considered as outliers. We also consider GBDT, a tree based classifier optimizing objective function through a multiple boosting process, and it is a good fit for REM [31].
In terms of the evaluation metrics, we use the F1 score to evaluate the model’s classification ability, and we set 0.5 as the decision threshold. Another adopted evaluation metric is the Precision-Recall AUC (PR-AUC), which focuses on the trade-off between precision and recall ability [20]. Besides the F1 and PR-AUC score for classification ability, we also care about the model’s ability with uncertainty prediction as the model results provide support for decision-making [40]. We use the negative log likelihood (NLL), Brier score [40]. Considering the extreme imbalanced label in REM, we modify the Brier score as the Brier Skill Score (BSS) in the following:
where $\Pr ({y^{(i)}}|{\mathbf{x}^{(i)}})$ is predicted probability of label ${y^{(i)}}$ and $\operatorname{BSS}\in [0,1]$. Similar as the BSS, we modify the NLL with the relative NLL (short as R-NLL) given the imbalanced label:
where $\operatorname{NLL}\Big(\Pr ({y^{(i)}}|{\mathbf{x}^{(i)}}),{y^{(i)}}\Big)$ measure the NLL between $\Pr ({y^{(i)}}|{\mathbf{x}^{(i)}})$ and label ${y^{(i)}}$.
(4.8)
\[ \operatorname{BSS}=1-\frac{\frac{1}{n}{\textstyle\textstyle\sum _{i=1}^{n}}{\Big(\Pr ({y^{(i)}}|{\mathbf{x}^{(i)}})-{y^{(i)}}\Big)^{2}}}{\frac{1}{n}{\textstyle\textstyle\sum _{i=1}^{n}}{\Big(\frac{1}{n}{\textstyle\textstyle\sum _{i}^{n}}{y^{(i)}}-{y^{(i)}}\Big)^{2}}}\](4.9)
\[ R‐NLL=\frac{{\textstyle\textstyle\sum _{i=1}^{n}}\operatorname{NLL}\Big(\Pr ({y^{(i)}}|{\mathbf{x}^{(i)}}),{y^{(i)}}\Big)}{{\textstyle\textstyle\sum _{i=1}^{n}}\operatorname{NLL}\Big(\frac{1}{n}{\textstyle\textstyle\sum _{i}^{n}}{y^{(i)}},{y^{(i)}}\Big)}\]5 Benchmark Dataset
We investigated VI-REM’s generalization performance on unstructured image and natural language processing tasks and evaluated robustness. We conducted the experiments via Pytorch on an NVIDIA GeForce RTX 4080 Laptop GPU. To ensure fair comparisons and convenience of visualization, we set the number of latents to 2 in all experiments. As the two considered benchmark datasets are well established in practice, we design an ablation study to understand the functionality of each component in VI-REM. For model VI-REM-V1, we alternate the EP distribution with the Gaussian distribution. Model VI-REM-V2 further replaces the MANN component with the fully connected MLPs. Model VI-REM-V3 alternates the weight $\beta ={n_{0}}/{n_{1}}$ in Equation 3.9 with a small number. For simplicity, we set the weight $\beta =5$ in VI-REM-V3.
5.1 Benchmark Dataset: Imbalanced Digit Recognition
Figure 3
Learned representations via VI-REM under different tasks on original testing. Red points represent class Event, green points represent class Normality. Red points represent the class Event, green points represent the class Normality.
We applied the VI-REM to MNIST1 dataset. The MNIST dataset covers handwritten digit numbers 0 to 9, corresponding to class 0 to class 9 as the label. Following a similar setup in [3], in each task, one of the ten classes is set as Event, while the samples from the remaining nine classes represent Normality. In the model training part, 100 images are randomly selected from the Event class, and nearly 600,000 images are selected from Normality class. As the selection of Event cases involves randomness, we repeat the experiment ten times for each task. In total, we conduct 100 individual runs of experiments. For model testing part, we create balanced testings consist of 5,000 images from the class Event and 5,000 images from the class Normality. Regarding robustness evaluations, we applied the random rotation operator to the original testing data. Each of the original testing images will rotate at a different angle, and the rotated angle is a random value from a uniform distribution $[-{60^{\circ }},{60^{\circ }}]$.
For VI-REM, the experiment sets $500\times 500\times 2$ fully connected (FC) MLPs for ${q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})$ and $500\times 500\times 500$ for MANN in recognizing ${p_{\boldsymbol{\theta }}}(y|\mathbf{z})$. We used Leaky ReLU for the activation function for both ${q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{x})$ and ${p_{\boldsymbol{\theta }}}(y|\mathbf{z})$. For S-MLP, F-MLP and MAML, we set the neural network to be $500\times 500\times \times 2\times 500\times 500\times 500\times 2$ FC layers; for D-SVDD, we set the neural network to be $1000\times 1000\times 1000\times 2$ FC layers. For the GBDT model, we set the maximum depth of GBDT to be 10 and number of estimators in GBDT to be 50.
Table 1
Digit Number Detection Overall Performance Comparison.
| R-NLL↑ | BSS↑ | PR-AUC↑ | F1 Score↑ | ||
| Original | GBDT | 0.28 | −0.49 | 0.75 | 0.24 |
| MAML | 1.32 | 0.47 | 0.98 | 0.82 | |
| S-MLP | 0.74 | 0.31 | 0.98 | 0.78 | |
| F-MLP | 0.83 | 0.40 | 0.98 | 0.81 | |
| D-SVDD | 0.38 | −0.80 | 0.60 | 0.01 | |
| VI-REM | 1.13 | 0.51 | 0.96 | 0.86 | |
| VI-REM-V1 | 1.25 | 0.54 | 0.95 | 0.86 | |
| VI-REM-V2 | 0.63 | −0.22 | 0.78 | 0.45 | |
| VI-REM-V3 | 0.55 | −0.43 | 0.70 | 0.32 | |
| Rotation | GBDT | 0.24 | −0.76 | 0.65 | 0.24 |
| MAML | 0.23 | −0.82 | 0.69 | 0.12 | |
| S-MLP | 0.31 | −0.50 | 0.77 | 0.33 | |
| F-MLP | 0.25 | −0.70 | 0.78 | 0.22 | |
| D-SVDD | 0.38 | −0.81 | 0.54 | 0.01 | |
| VI-REM | 0.87 | 0.21 | 0.89 | 0.71 | |
| VI-REM-V1 | 0.34 | −0.45 | 0.72 | 0.36 | |
| VI-REM-V2 | 0.26 | −0.75 | 0.65 | 0.17 | |
| VI-REM-V3 | 0.23 | −0.92 | 0.57 | 0.05 |
↑ for higher is better and ↓ for lower is better.
Figure 3 shows the learned representations via VI-REM under different tasks under the original testing. It’s clear that all the class Event deviate from the distribution and mostly lie on the tail part. This pattern is consistent with the EP prior from the VI-REM model. Table 1 reports the average of the 100 experiment results for both original testings and robustness evaluations. The proposed VI-REM model consistently outperformed the competing models regarding model accuracy and robustness. Under random rotation tasks, the VI-REM shows significantly improved model robustness compared to other SOTAs. Through ablation study, the MANN structure largely increases the model’s generalization ability under original testings, while the EP distribution contributes much larger in robustness.
Table 2
Digit Number Detection Overall Performance Comparison with Different Number of Latent Variables.
| Dim | R-NLL↑ | BSS↑ | PR-AUC↑ | F1 Score↑ | |
| Original | 2 | 1.13 | 0.51 | 0.96 | 0.86 |
| 3 | 0.33 | 0.31 | 0.96 | 0.77 | |
| 5 | 0.48 | 0.56 | 0.96 | 0.87 | |
| Rotation | 2 | 0.87 | 0.21 | 0.89 | 0.71 |
| 3 | 0.13 | −0.42 | 0.77 | 0.38 | |
| 5 | 0.13 | −0.35 | 0.77 | 0.43 |
↑ for higher is better and ↓ for lower is better.
The dimension of latent variable is an important hyper parameter in the VI-REM. To examine how the latent variable dimension affects the performance of our model, we vary it from 2 to 3 and 5. The results are reported in Table 2, and the column “Dim” displays the dimension for the latent variables. We find that setting the dimension of the latent variable can generate satisfying results on both original testing and rotation testing.
5.2 Benchmark Dataset: Imbalanced News Detection
This experiment applied the VI-REM to the standard natural language processing (NLP) benchmark dataset AG News Classification.2 The AG News dataset covers various news articles from “Sports”, “World”, “Business”, and “Sci/Tech” these four different categories. For implementing the applications, we used pre-trained fastText [33] and SWEM [38] to extract raw text representations. For the original input text data, we firstly apply the pre-trained fastText model to transform raw word input into the word embedding, where we set the dimensionality to 300. Based on the extracted word embedding, we used the average operator from SWEM [38] to combine word embedding and generate a $300\times 1$ vector for an individual news article. In terms of robustness evaluations, we randomly select 20% words from class Event and class Normality, and exchange them. Figure 4 gives an illustrative example of the data processing and robustness evaluations.
We created an extremely imbalanced text classification task. For the training task, we treated 30,000 “Sports” news articles as Normality and 100 randomly selected news articles from other categories as Event types. Within each task, we repeat the random samplings for Event types ten times to alleviate randomness. In task 1, we use the news articles from the “World” category as class Event while in task 2 and 3, we set “Business” and “Sci/Tech” separately as Event types. For out-of-sample model performance evaluation, 1,900 news articles for each category made up the testing sample. For VI-REM, we set a $500\times 500\times 2$ FC neural layer for ${q_{\phi }}(\mathbf{z}|\mathbf{x})$ and $500\times 500\times 1$ for the individual function ${f_{j}}({\mathbf{z}_{j}})$ of MANN in ${p_{\theta }}(y|\mathbf{z})$. For S-MLP, F-MLP and MAML, we set the neural network as $500\times 500\times 2\times 500\times 500\times 1$ FC layers; for D-SVDD, we set a $1000\times 1000\times 1000\times 2$ network.
Table 3 reports average performance under the original testings and robustness evaluations, and VI-REM shows better generalization ability and more robust model performance in dealing with data perturbations. The MANN structure largely increases the model’s generalization ability through the ablation study in Table 3. We also check how the model performances will change when setting a different number of latent dimension in Table 4.
Table 3
News Detection Model performance Comparisons.
| Method | R-NLL ↑ | BSS ↑ | PR-AUC↑ | F1-Score ↑ | |
| Original | GBDT | 0.06 | $-0.61$ | 0.81 | 0.33 |
| MAML | 0.60 | 0.44 | 0.99 | 0.82 | |
| F-MLP | 0.34 | 0.43 | 0.99 | 0.83 | |
| S-MLP | 0.35 | 0.43 | 0.98 | 0.83 | |
| D-SVDD | 0.45 | $-0.69$ | 0.60 | 0.02 | |
| VI-REM | 1.64 | 0.58 | 0.98 | 0.87 | |
| VI-REM-V1 | 1.43 | 0.52 | 0.98 | 0.84 | |
| VI-REM-V2 | 0.68 | 0.20 | 0.95 | 0.63 | |
| VI-REM-V3 | 0.50 | −0.03 | 0.92 | 0.50 | |
| Mixture | GBDT | 0.05 | $-0.95$ | 0.66 | 0.06 |
| MAML | 0.34 | 0.02 | 0.97 | 0.65 | |
| F-MLP | 0.11 | $-0.49$ | 0.86 | 0.39 | |
| S-MLP | 0.12 | −0.44 | 0.85 | 0.43 | |
| D-SVDD | 0.42 | $-0.75$ | 0.57 | 0.01 | |
| VI-REM | 0.92 | 0.24 | 0.95 | 0.74 | |
| VI-REM-V1 | 0.80 | 0.15 | 0.95 | 0.69 | |
| VI-REM-V2 | 0.42 | −0.18 | 0.87 | 0.47 | |
| VI-REM-V3 | 0.32 | −0.38 | 0.85 | 0.32 |
↑ for higher is better and ↓ for lower is better.
Table 4
News Detection Overall Performance Comparison with Different Number of Latent Variables.
| Dim | R-NLL↑ | BSS↑ | PR-AUC↑ | F1 Score↑ | |
| Original | 2 | 1.64 | 0.58 | 0.98 | 0.87 |
| 3 | 1.06 | 0.42 | 0.96 | 0.85 | |
| 5 | 0.91 | 0.06 | 0.99 | 0.27 | |
| Mixture | 2 | 0.92 | 0.24 | 0.95 | 0.74 |
| 3 | 0.65 | 0.04 | 0.91 | 0.70 | |
| 5 | 0.69 | −0.17 | 0.95 | 0.19 |
↑ for higher is better and ↓ for lower is better.
6 Real Data Applications
Two real data applications are considered. First application is to detect the fraudulent credit card transactions. The second application identifies traffic crash events based on driving data.
6.1 Application Case 1: Fraud Detection
We implemented the VI-REM model on the fraud detection applications3 where fraudulent credit card transactions are the rare events of interest. The dataset to which the VI-REM was applied consists of about 284,000 financial transactions associated with 28 predictive features, among which only 492 examples are labeled as Event resulting in an extremely low prevalence of $0.2\% $. We adopted a 5:5 train test split. To further evaluate model robustness in handling massive data noise, we created an artificially noisy setting by injecting Gaussian noise $\mathcal{N}(0,1)$ to the testing data, labeled as “Noisy 1” in Table 5. Further, we impose Gaussian noise $\mathcal{N}(0,2)$ to testing data and we refer it as “Noisy 2” in Table 5. The neural network structure is the same as Imbalanced News Detection.
Figure 5a and Figure 5b display the learnt representations for the fraud detection dataset in the original testing and noisy testings. Figure 5a and Figure 5b confirm that for the VI-REM model, the learnt representation for the Event cases lie on the tailed part while Normality mostly lie in the center.
Figure 5
The latent representation by different model in testing Set for Fraud Detection. PCA is short as Principal Component Analysis. For VI-REM model, the latents comes from $\mathbf{z}\sim {q_{\boldsymbol{\phi }}}(\mathbf{z}|\mathbf{z})$. The red points represent the fraud transactions, the green points represent the normal ones.
Table 5
Fraud Detection Model Comparison.
| Method | R-NLL ↑ | BSS ↑ | PR-AUC↑ | F1-Score ↑ | |
| Original | GBDT | 0.54 | 0.42 | 0.64 | 0.71 |
| MAML | 1.06 | $-0.31$ | 0.71 | 0.49 | |
| F-MLP | 2.23 | 0.63 | 0.79 | 0.79 | |
| S-MLP | 0.89 | 0.56 | 0.76 | 0.77 | |
| D-SVDD | 1.09 | $-0.01$ | 0.03 | 0.00 | |
| VI-REM | 3.72 | 0.68 | 0.80 | 0.81 | |
| Noisy 1 | GBDT | 0.07 | $-5.49$ | 0.36 | 0.17 |
| MAML | 0.54 | $-2.16$ | 0.67 | 0.26 | |
| F-MLP | 0.76 | $-0.61$ | 0.68 | 0.44 | |
| S-MLP | 0.31 | $-1.56$ | 0.46 | 0.32 | |
| D-SVDD | 0.99 | $-0.01$ | 0.02 | 0.00 | |
| VI-REM | 2.83 | 0.57 | 0.72 | 0.73 | |
| Noisy 2 | GBDT | 0.02 | $-30.2$ | 0.33 | 0.04 |
| MAML | 0.30 | $-4.46$ | 0.53 | 0.17 | |
| F-MLP | 0.28 | $-2.82$ | 0.41 | 0.22 | |
| S-MLP | 0.13 | $-5.10$ | 0.21 | 0.14 | |
| D-SVDD | 0.55 | $-0.26$ | 0.01 | 0.00 | |
| VI-REM | 1.32 | 0.16 | 0.53 | 0.56 |
↑ for higher is better and ↓ for lower is better.
Table 5 reports results for both original and noisy settings. The proposed VI-REM model consistently outperformed the competing models in terms of accuracy and robustness. Under noisy testing tasks, which sometimes can be associated with data quality issues in real-world settings, the VI-REM shows significantly improved model robustness compared to all SOTAs. VI-REM’s supreme performance on robustness to data noises indicates it can learn ubiquitous discriminative information. We notice that under “Noisy 1” testings, our model VI-REM still has reliable prediction performance while the other competitors are like a random guesses.
Next, we investigate the relationship between latents z and label y. Figure 6a plots the decision boundary and associated predicted probability for determining the Event cases under original testings. It is clear that the larger distance with the distribution center, the more likely it will be class Event. Figure 6b details the relationship between latent factors ${\mathbf{z}_{1}}$, ${\mathbf{z}_{2}}$ and label y. Latent ${\mathbf{z}_{1}}$ has a negative relationship with the occurrences of the events while ${\mathbf{z}_{2}}$ has a positive relationship.
We identify the relationship between input features x and latent z. We consider the 28 individual features and the interaction of any two of the input features. In total, $\left(\genfrac{}{}{0.0pt}{}{28}{1}\right)+\left(\genfrac{}{}{0.0pt}{}{28}{2}\right)=406$ combinations are considered for model interpretations. We implement the strategy discussed in Section 4.1 where we $q=5$ for each round of selections. Table 6 reports the major components of latent factor ${\mathbf{z}_{1}}$ and ${\mathbf{z}_{2}}$ respectively, where ${V_{22}}$ represents the 22th variable of the input features. The most contributed part of ${\mathbf{z}_{1}}$ and ${\mathbf{z}_{2}}$ do not overlap, confirming the disentanglement. We also fit a linear regression between ${\mathbf{z}_{1}}$ and ${\mathbf{z}_{2}}$ and their selected major components, and it is reported as column “Coefficients” in Table 6. For example, the fitted coefficient between the combination of ${V_{16}}$, ${V_{24}}$ and latent ${\mathbf{z}_{1}}$ is 0.024, indicating they have a positive relationship.
Table 6
Model Interpretations for Latent Factor ${\mathbf{z}_{1}}$ and ${\mathbf{z}_{2}}$ in Application Case 1 and Application Case 2.
| Selected Variable | Coefficient | ||
| Case 1 | Latent Variable ${\mathbf{z}_{1}}$ | ${V_{22}}\times {V_{26}}$ | −0.001 |
| Latent Variable ${\mathbf{z}_{1}}$ | ${V_{16}}\times {V_{24}}$ | 0.024 | |
| Latent Variable ${\mathbf{z}_{1}}$ | ${V_{16}}\times {V_{18}}$ | −0.041 | |
| Latent Variable ${\mathbf{z}_{1}}$ | ${V_{13}}\times {V_{24}}$ | −0.006 | |
| Latent Variable ${\mathbf{z}_{1}}$ | ${V_{16}}\times {V_{23}}$ | 0.011 | |
| Case 1 | Latent Variable ${\mathbf{z}_{2}}$ | ${V_{3}}\times {V_{13}}$ | 0.010 |
| Latent Variable ${\mathbf{z}_{2}}$ | ${V_{7}}\times {V_{13}}$ | −0.002 | |
| Latent Variable ${\mathbf{z}_{2}}$ | ${V_{13}}\times {V_{19}}$ | 0.040 | |
| Latent Variable ${\mathbf{z}_{2}}$ | ${V_{3}}\times {V_{16}}$ | 0.008 | |
| Latent Variable ${\mathbf{z}_{2}}$ | ${V_{13}}\times {V_{20}}$ | 0.023 | |
| Case 2 | Latent Variable ${\mathbf{z}_{1}}$ | ${\text{acc-x}_{26}}\times {\text{acc-y}_{26}}$ | −0.747 |
| Latent Variable ${\mathbf{z}_{1}}$ | ${\text{acc-z}_{26}}\times {\text{acc-z}_{29}}$ | −2.042 | |
| Latent Variable ${\mathbf{z}_{1}}$ | ${\text{acc-z}_{26}}\times {\text{acc-z}_{28}}$ | 0.133 | |
| Latent Variable ${\mathbf{z}_{1}}$ | ${\text{acc-z}_{25}}\times {\text{acc-z}_{27}}$ | −0.973 | |
| Latent Variable ${\mathbf{z}_{1}}$ | ${\text{acc-z}_{26}}\times {\text{acc-z}_{30}}$ | 1.896 |
Figure 6
Model interpretations for VI-REM in fraud detection. (a)The decision boundary and associated predicted probability for determining the Event cases under original testings. (b) Latent factors plot. Part (1) and part (3) of displays the probability density functions of latent variable 1 and latent variable 2, where the red lines represent the Event class and the green lines represent the Normality class.Part (3) and (4) presents the marginal relationship of the latent factors ${\mathbf{z}_{1}}$, ${\mathbf{z}_{2}}$ with the Logit value of condition probability, that is $\log \Big(\frac{\Pr (y|{\mathbf{z}_{1}})}{1-\Pr (y|{\mathbf{z}_{1}})}\Big)$ and $\log \Big(\frac{\Pr (y|{\mathbf{z}_{2}})}{1-\Pr (y|{\mathbf{z}_{1}})}\Big)$ for ${\mathbf{z}_{1}}$ and ${\mathbf{z}_{2}}$.
6.2 Application Case 2: Traffic Crash Identification
Traffic crashes are rare events, with an average rate of one crash 6.8 per million vehicle miles traveled in the US [10]. The naturalistic driving study (NDS) provides an unprecedented opportunity to evaluate crash risk [10]. NDSs are characterized by continuously recording driving information, such as three-dimensional Inertial Measurement Unit acceleration and multi-channel video recordings. Accurately identifying crashes with robustness benefits further understanding of driving behavior with the overall goal of reducing traffic accidents [10], identifying high risky segments at road network [52] or improving the validation efficiency for autonomous vehicles [34].
This experiment applied the VI-REM to the Second Strategic Highway Research Program (SHRP2) NDS [10], which is largest NDS to-date, with more than 1 million hours of continuous driving data. SHRP2 Insight website4 gives a more detailed introduction. We applied the VI-REM to identify traffic crashes based on the high-frequency longitudinal, lateral and vertical tri-axial acceleration for the vehicles.
SHRP2 crash events have four different severity levels [10], ranging from level 4 to level 1. Level 4 (L4) crash events are low-risk tire strikes, level 3 (L3) crash events are minor crashes, level 2 (L2) crash events are police reportable crashes with at least 1,500 dollars worth of damage, and level 1 (L1) are fatal crash events. Our data consists of 100 level 1, 150 level 2, 578 level 3, and 588 level 4 crash events. After a random selection process, we also collected nearly 120,000 safe driving segments from SHRP2. These selected safe driving segments have similar maximum and minimum speeds with crashes. Table 7 summarizes the detail of different events and the last column calculates the imbalance ratio.
Table 7
Summary Description of Level 1 to Level 4 Crash Events and Normal Drivings.
| Description | Number | Ratio | |
| L1 Crash | Fatal Crash | 100 | 1200:1 |
| L2 Crash | Police Reportable Crashes | 150 | 800:1 |
| L3 Crash | Minor Crash and MC | 578 | 208:1 |
| L4 Crash | Low-risk Tire Strikes | 588 | 204:1 |
| Normal Driving | Routine Safe Driving | 119,996 | 1:1 |
We applied the ordinal-based VI-REM model to the SHRP2 data, as the label of events has different severity levels. To standardize the data, we put the most volatile time in the middle of each driving segment. Each driving segment has 5.1 seconds and as the collection frequency for acceleration is 10 HZ, so the length for each segment is $5.1\times 10=51$. As the acceleration has three directions, the dimensionality for the input features will be $51\times 3=153$. Figure 7a shows several driving segments. Our experiment adopted a 5/5 split for training and testing, and the neural network structure for each model is the same as Imbalanced Digit Recognition.
Figure 7
Several plots for identifying traffic crash events in SHRP2 data (a) Several driving segments are associated with the three-dimensional accelerations under different event categories. (b) Part (1) and Part (2) represent the kernel density estimation of latent ${\mathbf{z}_{1}}$ and ${\mathbf{z}_{2}}$ given different event categories. Part (3) and Part (4) display the associated probability of different event categories given different values on ${\mathbf{z}_{1}}$ and ${\mathbf{z}_{2}}$.
Our work transforms the assessment of an ordinal model into assessing the performance of a classification model based on different selections on the severity level. For example:
where $Y=({y_{1}},\dots ,{y_{N}})$ and $\widehat{\Pr }(Y\ge j)$ is the predicted probability from model in Equation 4.7. Similar transformations will also be applied to the other metrics. Table 8 compares the performances of different models for identifying different severity traffic crashes, where we skip the Deep-SVDD model as it is difficult to extend to ordinal data conditions. Table 8 shows VI-REM’s supreme performances over other competitive models, especially for severe level 1, level 2 and level 3 crash events. Overall for identifying crash events (level 4 up to level 1), our proposed VI-REM can have 0.77 in F1 score, while the most competitive F1 score from benchmark model is 0.68 by Sampling MLP.
(6.1)
\[ \begin{aligned}{}& \text{PR-AUC}\hspace{2.5pt}\Big(\widehat{\Pr }(Y\ge j),\mathcal{I}(Y\ge j)\Big)\hspace{1em}\\ {} & \text{and}\hspace{1em}j\in [\text{L1, L2, L3, L4}]\end{aligned}\]Analytical modeling for identifying traffic crashes is sensitive to kinematic noise (such as signal connection issues or different driving environments), where model performance can largely deteriorate as noise increases [30]. Considering this, where model robustness is fundamentally important for decision making, we created an artificially noisy setting where Gaussian noise $\mathcal{N}(0,0.05)$ was injected to the input accelerations. The model performance metrics are presented in Table 9. As the table shows, the proposed VI-REM significantly outperforms alternative models based on the evaluation metrics. Compared with benchmarks, the VI-REM shows much more robust disincentive power under noisy testing conditions.
As for the model interpretations, Figure 7b shows the mutual relationship between z and label y. ${\mathbf{z}_{1}}$ has larger discriminative over different crash events, and larger values in ${\mathbf{z}_{1}}$ are associated with more considerable severity. We implement the same strategy in application case 1 to identify the relationship between input features x and latent z. We report the identified major components of ${\mathbf{z}_{1}}$ in Table 6, as latent ${\mathbf{z}_{2}}$ does not have strong distinguishing ability. In Table 6, ${\text{acc-x}_{26}}$ represents the longitudinal acceleration at point 26, and $\text{acc-y}$, $\text{acc-z}$ represent the lateral and vertical acceleration respectively. All the major contributing input features are around the middle of the driving segment through the model interpretations. This phenomenon is consistent with our experiment setup. We also find that vertical accelerations could be the most influential factor for identifying traffic crashes and determining their associated severity level. This observation is consistent with various empirical studies [30]
Table 8
Performance Comparisons in Traffic Crash Identifications.
| Method | R-NLL ↑ | BSS ↑ | PR-AUC↑ | F1-Score ↑ | |
| VI-REM | $\text{Y}\ge \text{L4}$ | 1.58 | 0.58 | 0.70 | 0.77 |
| $\text{Y}\ge \text{L3}$ | 2.50 | 0.52 | 0.69 | 0.68 | |
| $\text{Y}\ge \text{L2}$ | 3.32 | 0.59 | 0.77 | 0.74 | |
| $\text{Y}\ge \text{L1}$ | 2.71 | 0.56 | 0.78 | 0.71 | |
| GBDT | $\text{Y}\ge \text{L4}$ | 0.65 | $-0.23$ | 0.49 | 0.52 |
| $\text{Y}\ge \text{L3}$ | 0.48 | $-0.86$ | 0.27 | 0.31 | |
| $\text{Y}\ge \text{L2}$ | 0.62 | $-0.61$ | 0.38 | 0.36 | |
| $\text{Y}\ge \text{L1}$ | 0.94 | $-0.21$ | 0.50 | 0.42 | |
| MAML | $\text{Y}\ge \text{L4}$ | 1.91 | 0.42 | 0.64 | 0.56 |
| $\text{Y}\ge \text{L3}$ | 2.16 | 0.44 | 0.67 | 0.64 | |
| $\text{Y}\ge \text{L2}$ | 1.67 | 0.40 | 0.65 | 0.61 | |
| $\text{Y}\ge \text{L1}$ | 1.42 | 0.30 | 0.53 | 0.49 | |
| S-MLP | $\text{Y}\ge \text{L4}$ | 1.43 | 0.26 | 0.57 | 0.65 |
| $\text{Y}\ge \text{L3}$ | 1.36 | 0.08 | 0.63 | 0.66 | |
| $\text{Y}\ge \text{L2}$ | 1.68 | 0.24 | 0.54 | 0.54 | |
| $\text{Y}\ge \text{L1}$ | 1.44 | 0.09 | 0.57 | 0.51 | |
| F-MLP | $\text{Y}\ge \text{L4}$ | 0.71 | $-0.01$ | 0.70 | 0.68 |
| $\text{Y}\ge \text{L3}$ | 1.25 | 0.40 | 0.65 | 0.60 | |
| $\text{Y}\ge \text{L2}$ | 1.15 | 0.55 | 0.72 | 0.68 | |
| $\text{Y}\ge \text{L1}$ | 0.66 | 0.49 | 0.62 | 0.58 |
↑ for higher is better and ↓ for lower is better.
Table 9
Robustness Evaluations in Traffic Crash Identifications.
| Method | R-NLL ↑ | BSS ↑ | PR-AUC↑ | F1-Score ↑ | |
| VI-REM | $\text{Y}\ge \text{L4}$ | 1.48 | 0.48 | 0.69 | 0.73 |
| $\text{Y}\ge \text{L3}$ | 2.43 | 0.47 | 0.68 | 0.66 | |
| $\text{Y}\ge \text{L2}$ | 3.28 | 0.59 | 0.77 | 0.74 | |
| $\text{Y}\ge \text{L1}$ | 2.57 | 0.57 | 0.78 | 0.72 | |
| GBDT | $\text{Y}\ge \text{L4}$ | 0.38 | $-1.09$ | 0.42 | 0.38 |
| $\text{Y}\ge \text{L3}$ | 0.38 | $-1.42$ | 0.24 | 0.23 | |
| $\text{Y}\ge \text{L2}$ | 0.38 | $-0.84$ | 0.31 | 0.31 | |
| $\text{Y}\ge \text{L1}$ | 0.67 | $-0.64$ | 0.43 | 0.38 | |
| MAML | $\text{Y}\ge \text{L1}$ | 1.85 | 0.40 | 0.58 | 0.56 |
| $\text{Y}\ge \text{L2}$ | 1.99 | 0.39 | 0.64 | 0.62 | |
| $\text{Y}\ge \text{L3}$ | 1.60 | 0.38 | 0.61 | 0.60 | |
| $\text{Y}\ge \text{L4}$ | 1.43 | 0.30 | 0.52 | 0.49 | |
| S-MLP | $\text{Y}\ge \text{L4}$ | 0.87 | 0.05 | 0.58 | 0.61 |
| $\text{Y}\ge \text{L3}$ | 0.66 | $-0.65$ | 0.62 | 0.63 | |
| $\text{Y}\ge \text{L2}$ | 0.87 | $-0.42$ | 0.45 | 0.48 | |
| $\text{Y}\ge \text{L1}$ | 0.63 | $-1.34$ | 0.45 | 0.29 | |
| F-MLP | $\text{Y}\ge \text{L4}$ | 0.55 | $-0.53$ | 0.64 | 0.60 |
| $\text{Y}\ge \text{L3}$ | 1.00 | 0.27 | 0.62 | 0.60 | |
| $\text{Y}\ge \text{L2}$ | 1.14 | 0.55 | 0.71 | 0.68 | |
| $\text{Y}\ge \text{L1}$ | 0.65 | 0.48 | 0.61 | 0.61 |
↑ for higher is better and ↓ for lower is better.
6.3 Discussion on the Running Time
Before concluding our paper, we present the average running time of our method VI-REM and other competitive methods across different tasks in Table 10. Column “CV” represents the task for imbalanced digit recognition in Section 5.1, column “NLP” represents the task for imbalanced detection in Section 5.2, and column “Fraud” represents the task for fraud detection in Section 6.1, and “Crash” represents the task in Section 6.2. Our experiments show that MAML is the most robust among the competitive methods. In terms of efficiency, VI-REM is faster than MAML, especially under the MLP and CV task. We also find that most of the running time for the VI-REM model comes from the MANN part. Designing a more efficient regularized neural networks could be an interesting topic to explore in the future. The GBDT method performs fast under CV and Fraud task, but may experience longer iteration (with an average of 2113 seconds) to converge when dealing with complex text data. Under the OLM setup in Section 6.2, the VI-REM is even faster than S-MLP.
Table 10
Overall Computation Time across Different Methods on Different Tasks. The unit is $\times 100$ Second.
| Model | CV | NLP | Fraud | Crash |
| GBDT | 1.05 | 21.13 | 1.01 | 6.44 |
| MAML | 4.82 | 1.53 | 2.78 | 3.79 |
| F-MLP | 0.30 | 0.11 | 0.35 | 0.32 |
| S-MLP | 0.51 | 0.14 | 0.32 | 3.53 |
| D-SVDD | 0.25 | 0.19 | 0.32 | - |
| VI-REM | 3.17 | 0.80 | 2.40 | 3.04 |
| VI-REM-V1 | 2.95 | 0.77 | 2.10 | 2.89 |
| VI-REM-V2 | 0.28 | 0.20 | 0.16 | 0.27 |
7 Summary and Discussion
This study proposes a Variational Inference toward extreme approach for rare event modeling (short as VI-REM). Our model induces uncertainty to representation learning through a novel EP distribution and addresses the over-fitting issue through MANN, which leads to a more robust model for rare events. We provide theoretical properties on the efficacy of the proposed model. Extensive empirical experiments are conducted to confirm the model performance over various application scenarios with promising results, especially with respect to model generalization, robustness and interpretability.
The VI-REM framework greatly improves the generalizability and interpretability of rare event modeling, two challenging issues associated with limited numbers of events. The application of VI-REM could accurately depict the risk associated with rare-events and allow the general public and decision-makers to set realistic expectations for rare events. The identification of adverse events at an early stage may allow mitigation of the damage and loss associated with the events. The features identified through the model are crucial for researchers and practitioners to identify the causes of rare events and take proper countermeasures to prevent and reduce the occurrence of future adverse events.