\section{Introduction}

Neutrinoless double beta decay consists in the transformation of an even-even nucleus ($A$,~$Z$) into a  lighter isobar ($A$,~$Z+2$) accompanied by the emission of 2 electrons and no other particle [1]. This process is forbidden within the Standard Model since its discovery would confirm that lepton number is not a symmetry of nature. In addition, it will imply that neutrinos are the only fermions that coincide with their antiparticle, i.e. Majorana particles. This process is energetically allowed for 35 nuclei ($^{76}$Ge, $^{82}$Se, $^{100}$Mo, $^{116}$Cd, $^{130}$Te, $^{136}$Xe ...) [2] and its signature is a monochromatic peak at the $Q$-value ($Q_{\beta\beta}$) of the transition corresponding to the sum of the kinetic energies carried by the two emitted electrons. The most promising and desirable isotopes to study are those with a high $Q_{\beta\beta}$ above the natural $\gamma$ radioactivity endline ($^{208}$Tl, ~ 2.615 MeV), since this will imply less background in the region of interest (ROI) and higher decay probability (since the phase-space of the process increases with higher $Q_{\beta\beta}$). Another important consideration is the natural abundance for the $0\nu\beta\beta$ isotope. High isotopic abundance (I.A.) will translate into easiness of enrichment of the $\beta\beta$ isotope and a lower cost as well. The current most stringent half-life limits on $0\nu\beta\beta$ decay are of the order of 10$^{25}$–10$^{26}$ y [3]. This makes any experiment aiming at searching for this rare decay extremely challenging.
An experiment has to be performed with radiopure large-mass high-resolution detectors and negligible radioactivity from anything except the nucleus under study, because otherwise the signal will be buried in the background.

\section{Bolometric technique and choice}
Bolometers are very powerful nuclear detectors to study $0\nu\beta\beta$ [4]. When a particle releases energy in the absorber (crystal) it produces lattice vibrations that generate a temperature rise in the crystal measured by a thermometer glued on the crystal surface. Working at cryogenic temperature is crucial to lower the heat capacity of the insulating crytal ({\it C~$\propto$~T$^3$}), because otherwise the signal amplitude will be dimmed. Bolometers offer wide choices of absorber materials making it possible to improve the sensitivity of the experiment through the maximization of the efficiency and the minimization of the background by choosing a compound compatible with an active background-rejection technique (pulse shape discrimination or detection of scintillation [5] or Cherenkov light [6]). The CROSS choice of the isotopes and the compound crystals is based on years of dedicated R$\&$Ds. Li$_{2}$MoO$_{4}$ [7,8] and TeO$_{2}$ [9] were found to be among the most promising detector materials, containing the $0\nu\beta\beta$ isotopes $^{100}$Mo ($Q_{\beta\beta}$~=~3034 keV $\&$ I.A. = 9.6$\%$) and $^{130}$Te ($Q_{\beta\beta}$ = 2527 keV $\&$ I.A. = 34$\%$) respectively. $^{100}$Mo fulfills the major requirement of a high $Q_{\beta\beta}$ above 2.615 MeV. $^{130}$Te is an attractive candidate for being the $0\nu\beta\beta$ isotope with the highest natural abundance. Li$_{2}$MoO$_{4}$ and TeO$_{2}$ have successfully been studied in LUMINEU and CUORE projects respectively and shown to have a high energy resolution in ROI ($\bigtriangleup_{FWHM}$ = 5 keV) and high internal radiopurity as an effective purification proccesses can be developed for these compounds [6,8,10].

\section{Background rejection technique}
A source of background, as mentioned in  the Introduction, comes from the highest intensity natural $\gamma$ line from $^{208}$Tl (2.615 MeV) in the $^{232}$Th decay chain, which can provide an important $\gamma$ background contribution in the $^{130}$Te ROI. Above 2.615 MeV, the $\gamma$ background contribution can come only from weak lines of $^{214}$Bi and it can be made negligible. Energy-degraded $\alpha$s coming from $^{238}$U and $^{232}$Th radioactive chains due to surface contamination of the crystal and the surrounding materials (copper holder, Polytetrafluoroethylene...) contribute to the background in  ROI in both crystals [8]. In addition, there is a contribution from $\beta$s emitted near the surface from the surrounding materials. Surface radioactive contamination is considered the most challenging limitation to the sensitivity for $0\nu\beta\beta$ bolometeric experiments because of the non-existence of a dead layer for bolometers. The rejection of $\alpha$s can be achieved by exploiting the scintillation and Cherenkov radiation (for Li$_{2}$MoO$_{4}$ and TeO$_{2}$ respectively) emitted by the absorber and detected by a separate light detector (which is also a bolometer made out of Ge wafers) facing the crystal, since $\alpha$s and $\beta$s have different light yield.\\

CROSS proposes another technique of surface background rejection which does not require an additional device. The key idea is to provide the crystal surface with pulse shape modification capability (in terms of the pulse rise-time for example), so that it would be possible to discriminate surface $\alpha$s from the bulk events ($\gamma/\beta$), which include the $0\nu\beta\beta$ signal. This is obtained by depositing a thin layer (few $\mu$m) of superconducting Al film (critical temperature $T_{c}$ = 1.2 K) via evaporation. The superconducting Al film acts as a pulse shape modifier for surface events up to $\sim$ 1 mm depth [18] from the crystal surface coated by the film. The feasibility of this method has been proven in a test done back in 2010 (CSNSM, Orsay) with a TeO$_{2}$ having NbSi thin film (Anderson insulator) [11,12] as a fast thermal sensor ( $\sim$ 1 ms time response).

The tests reported here on small crystals of Li$_{2}$MoO$_{4}$ and TeO$_{2}$ ($20\times20\times10$ mm$^3$) with one surface covered with a few-$\mu$m-thin Al film ($18\times18$ mm$^2$) show that the discrimination capability between surface and bulk events is successfully achieved when using NTD (neutron transmutation doped) Ge thermistor as a sensor (slow sensor) [13]. The NbSi and the NTD sensors showed to work in a different way. The NbSi film is deposited ($13\times13$ mm$^2$) directly on the crystal over a large surface ($14\times14$ mm$^2$, Fig. 1 top right). This makes NbSi sensor sensitive to the prompt athermal component of the phonon population produced by the impinging particle, while NTDs are sensitive rather to the thermal component due to their intrinsic slowness and the glue interface. When a particle releases energy in the absorber close to the Al-film-coated surface it generates high energy phonons that break Cooper pairs in the superconducting film. This effect is less important for bulk events, since the generated phonons reach the surface with a lower average energy due to the quasi-diffusive mode of phonon propagation [14]. The bulk events will be felt directly by the NbSi sensor and later after thermalization by the NTD sensor. For the surface events, a significant fraction of the particle energy will be trapped in the superconductive film in the form of quasiparticles for a few millisecond due to the long recombination time of quasiparticles at a temperature well below the transition energy of the superconducting Al film. The subsequent recombination of quasiparticles to a lower energy $\sim$ 1.2 K phonons will add a delayed component to the signal when read by the NbSi (making bulk events faster than surface events), while it will accelerate the detection in case of the thermal phonon sensitive sensor (NTD) with respect to the bulk events due to the faster thermalization (making surface events faster than bulk events).

\begin{figure}[htbp]
\centering
\includegraphics[width=0.44\linewidth]{presentation1.pdf}
\includegraphics[width=0.365\linewidth]{Presentation2.pdf}
\\
\includegraphics[width=0.3\linewidth]{LMO.jpg}
\includegraphics[width=0.285\linewidth]{big.jpg}
\includegraphics[width=0.3\linewidth]{TeO2.jpg}
\linebreak
\caption{({\it top left}) Photograph of the typical assembly of CROSS detectors, Li$_{2}$MoO$_{4}$ and TeO$_{2}$ (bolo1 and bolo3, see Table 1). ({\it top right}) The  TeO$_{2}$ crystal having two readouts NbSi and NTD Ge (bolo5), and 10 $\mu$m Al film on one 20$\times$5 mm$^{2}$ side. ({\it bottom left}) Small Li$_{2}$MoO$_{4}$ crystal with 10 $\mu$m Al film evaporated on 20$\times$20 mm$^{2}$ side. ({\it bottom middle})  Large Li$_{2}$MoO$_{4}$ crystal (bolo2) with the 10 $\mu$m Al film evaporated on one circular side (on the other side we glue an NTD and a heater). ({\it bottom right}) TeO$_{2}$ crystal with 10 $\mu$m Al film evaporated on 20$\times$20 mm$^{2}$ side. All the faces with Al film are facing a uranium $\alpha$ source (surface events) and for bolo5 another uranium $\alpha$ source deposited directly on the crystal surface with no Al film (the other  20$\times$5 mm$^{2}$ side without Al film).}
\label{fig:one}
\end{figure}

\section{Experimental results}
\subsection{Assembly description}
We have performed several runs on different bolometers (Fig.~\ref{fig:one} and Table~\ref{table:bolo}) with different sensors, film types and thicknesses to achieve the best discrimination capability. A thin film of Al was evaporated on a surface of the bolometers using a dedicated evaporator. An NTD Ge thermistor was used for all the  bolometers. It was glued on a crystal surface with either six 25-$\mu$m-thick glue spots or a thin layer of glue spread all under the NTD (bolo2). In addition, a heater was glued to provide periodically thermal pulses in order to stabilize the response of the bolometers [15]. The crystals were fixed on a copper plate by means of PTFE pieces and brass/copper screws. A uranium $\alpha$ source was placed facing the surface with Al film to generate surface events. The source provides two main $\alpha$ lines at $\sim $ 4.2~MeV and $\sim $~4.7~MeV, and also $\beta$ particles from the isotope $^{234m}$Pa, with a spectrum extending up to $\sim $ 2.27 MeV. The TeO$_{2}$ with NbSi (bolo5) sensor tested in the past [16] was tested again with a modification. An NTD was added as a second sensor to appreciate the different response between two sensors on the same crystal with sensitivities depending on the phonon energy. This crystal has $10~\mu$m Al film faced by a uranium $\alpha$ source (surface events) and another uranium $\alpha$ source deposited directly on the opposite surface to emulate bulk events (there is no pulse shape difference between $\alpha$ and $\beta/\gamma$ for  TeO$_{2}$, while such difference is expected for Li$_2$MoO$_4$ as often happens in bolometers based on crystal scintillators [17,4]).

\begin{table}[t]
\centering
\begin{tabular}{c c c c c c c c}
\hline
 & Bolometer & Sensor (glue) & \thead{Crystal size \\ (mm$^{3}$)} & \thead{Crystal \\ mass (g)} &  \thead{Al film \\ thickness} & DP & \thead{Sensitivity \\ (nV/keV)} \\
[0.5ex]
\hline

bolo1 & Li$_{2}$MoO$_{4}$ & NTD Ge (6 spots)   & 20$\times$20$\times$10    & 12   &  10 $\mu$m  &  12 & 53 \\

bolo2 & Li$_{2}$MoO$_{4}$  & NTD Ge (veil)  & $\varnothing$40$\times$20   & 67 & 10 $\mu$m & 7 & 37 \\

bolo3 & TeO$_{2}$   & NTD Ge (6 spots)   & 20$\times$20$\times$10   & 25 & 1 $\mu$m  & 3 & 43\\

bolo4 & TeO$_{2}$   & NTD Ge (6 spots)   & 20$\times$20$\times$10   & 25 &  10 $\mu$m & 4 & 44 \\

bolo5 &TeO$_{2}$    & \thead{NTD Ge (2 spots)  \\+\\ NbSi (evaporation)}  & 20$\times$20$\times$5  & 12.5 & 10 $\mu$m & & \thead{70 \\ \\ 54}   \\

\hline
\end{tabular}
\linebreak
\caption{Description of the five bolometers that have been tested and their respective DP (Eq. (1)) and sensitivity at 22~mK (bolo5 was operating at 30 mK). A uranium $\alpha$ source irradiated the film in all the detectors. DPs for bolo 5 are missing because the alpha population has a large spread presumably due to position effects. As a results, the distribution of the rise-time is strongly not gaussian and the DP as defined here becomes meaningless.}
\label{table:bolo}
\end{table}

\begin{figure}[t]
\centering
\includegraphics[width=0.47\linewidth]{TeO2_risetime_run46m13_dp.pdf}
\includegraphics[width=0.43\linewidth]{risetime_lmo_big_DP.pdf}
\linebreak
\caption{ The rise-time is plotted as a function of the energy: ({\it left}) TeO$_{2}$ (bolo3) with 1 $\mu$m Al film;  ({\it right}) Li$_{2}$MoO$_{4}$ (bolo2) with 10 $\mu$m Al film.}
\label{fig:two}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.45\linewidth]{NTD_risetime.pdf}
\includegraphics[width=0.446\linewidth]{NbSi_risetime.pdf}
\linebreak
\caption{ The rise-time is plotted as a function of the energy: TeO$_{2}$ with two readouts, NTD Ge ({\it left}) and NbSi ({\it right}).}
\label{fig:three}
\end{figure}

\subsection{Pulse shape discrimination results}
The tests have shown that the surface $\alpha$ events can be safely rejected by their different pulse shape from the bulk events in all of the bolometer listed in Table ~\ref{table:bolo}, while the identification of surface $\beta$ particles requires additional R$\&$D work. The PSD parameter that we discuss here is the rise-time (time from 10$\%$ to 90$\%$ of the pulse maximum amplitude). Fig.~\ref{fig:two}  shows the rise-time parameter plotted as a function of the deposited energy in the crystal for TeO$_{2}$ (bolo3) and Li$_{2}$MoO$_{4}$ (bolo2). The surface $\alpha$s (from the uranium source) are clearly separated from the bulk making their rejection possible. The discrimination capability can be quantified in terms of the discrimination power DP (Eq. (1)), where $\mu$ and $\sigma$ are the mean and standard deviation respectively extracted from the gaussian fit applied to the surface and bulk event distributions (insets of Fig.~\ref{fig:two}):

\begin{equation}
	DP=\frac{ |{\mu_{surface}-\mu_{bulk}}|}{\sqrt{\sigma_{surface}^2+\sigma_{bulk}^2}} ,
\end{equation}

The DP values for all the bolometer are given in Table ~\ref{table:bolo}. Comparing the results of bolo3 and bolo4, which are identical but the film thickness, there is a weak but not conclusive indication that thicker films may lead to a better discrimination. In addition, results showed that depositing Al film does not affect the energy resolution and sensitivity of the bolometer [18].

Fig.~\ref{fig:three} shows the discrimination results for TeO$_{2}$ bolometer instrumented with NTD and NbSi. It is confirmed that surface events are faster (slower) than the bulk read out by the NTD Ge sensor (NbSi sensor) as discussed in section 3. The two alpha lines in the bulk (from the uranium source deposited on the crytsal surface directly) are at their nominal energies (4.2 MeV and 4.7 MeV), whereas the surface alphas are either energy-overestimated (NTD sensor) or energy-underestimated (NbSi sensor). The former observation is because of the better thermal conversion of Al film (a superconductor compared to the crystal which is an insulator), while the later effect is because of lower pulse amplitude caused by a trapped energy in form of quasiparticles. Note that the calibration in both cases was done with the low energy gamma quanta of 352 keV from $^{238}$U decay chain. One can notice that the plots of Fig.~\ref{fig:three} show a linear behaviour (PSD parameters don't depend on deposited energy) for bolo5 unlike bolo2 and bolo3 (Fig.~\ref{fig:two}) and this is because bolo5 was operating at a higher temperature. It has to be noticed that the U source, due to the geometry of the set-up, can produce also a small population of bulk-like (the reason of the double poplution in the bulk) events when $\alpha$ particles hit directly the TeO$_2$ crystal face, which is not fully covered by the Al film [16].

The discrimination power depends on the pulse parameter used. Another pulse shape parameter was tested and was found to provide a higher DP's than the rise-time (a full description of this parameter and the related results can be found in [18]).

\section{Conclusions}
We have described in this paper the technology behind CROSS (Cryogenic Rare-event Observatory with Surface Sensitivity), which showed its capability of rejecting surface $\alpha$ background when a superconducting Al film is deposited on a side of the crystal. This method will allow us to simplify substantially the bolometric structure by getting rid of a light detector. Future goals will be to fully coat the crystals and operate a medium-scale demonstrator with tens of elements in the Canfranc underground laboratory (Spain).

\begin{acknowledgements}
The project CROSS is funded by the European Research Council (ERC) under the European-Union Horizon 2020 program (H2020/2014-2020) with the ERC Advanced Grant no. 742345 (ERC-2016-ADG). The PhD fellowship of H. Khalife has been partially funded by the P2IO LabEx (ANR-10-LABX-0038) managed by the Agence Nationale de la Recherche (France) in the framework of the 2017 P2IO doctoral call.
\end{acknowledgements}


\section{Introduction}\label{intro}}

\IEEEPARstart{C}{linical} data mining is the application of machine learning to medical databases as a tool to assist clinical research and decision. Its impact depends on one hand on the clinical pertinence of the models, but also on their predictive performance, interpretability, and appropriate integration in clinical software~\cite{Patel2009}. While many methodological successes exist, examples that affect patient management are seldom found~\cite{Clifton2015}. Measurement uncertainty limits the application of data mining, as it impairs model performance and generalizability~\cite{Roddick2003,Lavrac1998}. Improving data quality is however resource-intensive and often unfeasible~\cite{Karr2006}. We hypothesize that the acknowledgement of uncertainty, in particular by integrating domain-knowledge about the reliability of each measurement, can improve models and the leverage them as an asset for clinical research and decision making.

The uncertainty of a measurement reflects the lack of knowledge about its exact value, and is often caused by noise in the acquisition~\cite{jcgm2008evaluation}. Clinical measurement uncertainty originates from multiple sources including distinct diagnostic practices~\cite{Cios2002}, inter- and intra-observer variability~\cite{Singh2003}, manufacturer-dependent technology~\cite{Foley2012}, the use of distinct modalities for the same measurement~\cite{Lopez-Minguez2014}, or patient factors such as body habitus or claustrophobia, which affect the choice and the quality of an imaging test~\cite{genders2017quantitative}. E.g. in the determination of device size for left-atrial appendage closure, consistency between \ac{CT}, transesophageal echocardiography and angiography was observed in only 21.6\% of the cases~\cite{Lopez-Minguez2014}. In the estimation of \ac{EF}, the variability between \ac{CMR} and echocardiography resulted in 28\% of the population having opposing device eligibility~\cite{DeHaan2014}. Clinical reasoning involves making guideline-abiding decisions based on such unreliable data. In order to employ scientific evidence in practice, the experienced clinician assesses the reliability of each measurement, and integrates it with his/her training and experience. This endeavor is all the more challenging, considering the emphasis on internal validity of medical research, as opposed to external validity. The contrast between the scrutinous design of populations used for research and the actual populations where evidence is employed has been considered an obstacle to evidence-based medicine~\cite{green2014closing}.

\Acp{DT} are a knowledge-representation structure, where decisions at the internal nodes lead to a prediction at the leaves. They can be efficiently built using the well-known algorithms ID3~\cite{quinlan1979interactive}, C4.5~\cite{Quinlan1993}, CART~\cite{Breiman1984} or CHAID~\cite{kass1980exploratory}. Compared to other methods, \acp{DT} can be interpreted as a sequence of decisions. Although recent methods such as deep neural networks can offer better performance, their output is often a black box~\cite{Weng2017}. Interpretability is all the more necessary as recent European law secures the right to an explanation of all algorithmically-made decisions~\cite{EU}.

Learning \acp{DT} from noisy measurements can overfit to the noise and fail to generalize. Generating predictions for noisy instances can generate incorrect predictions. Each test at a \ac{DT} node compares a measurement with a \textit{hard threshold}, such that small errors can lead to the instance following an opposing path. Moreover, the distance between the measurement and the threshold is disregarded~\cite{Quinlan1987}. Several algorithms explore the idea of \textit{soft thresholds} to make \acp{DT} robust to the uncertainty~\cite{Quinlan1987,DvorakS07,Tsang2011,irsoy2012soft,Yuan1995}. Such approaches weight the contribution of all child branches to the prediction. Some methods focus on cognitive uncertainty, while others focus on statistical uncertainty or noise. The notion of \emph{fractional tuple} was first introduced in C4.5 for handling missing values~\cite{Quinlan1993}.

Fuzzy \acp{DT} use fuzzy logic to handle \textit{cognitive uncertainty}, which deals with inconsistencies in human reasoning~\cite{Yuan1995}. Using fuzzy \acp{DT} to handle statistical noise involves setting the fuzzy membership functions to express uncertainties around the \ac{DT} thresholds, known in advance. Alternatively, discretization algorithms can be used. But those methods are computationally costly~\cite{Wang2000}, and the accuracy of the resulting \acp{DT} strongly depends on the chosen algorithm~\cite{Segatori2017}.

Probabilistic \acp{DT}, on the other hand, focus on \textit{stochastic uncertainty} or random noise, arising from the unpredictable behavior of physical systems and measurement limitations:

\begin{itemize}[leftmargin=5mm]

	\item Quinlan~\cite{quinlan1990probabilistic} proposed a method using \acp{PLT} for evaluating a~\ac{DT}, assuming a two-modal uniform distribution for the noise. The parameters of this distribution are set using a statistical heuristic based on the training data. Dvor{\'{a}}k and Savick{\'{y}}~\cite{DvorakS07} employed a variation of this method, where the parameters were estimated through simulated annealing. Experiments in one dataset led to 2-3\% error rate reductions compared to CART~\cite{Breiman1984}, suggesting the potential of the method and the need for an evaluation on more data. No significant differences were found between the two parameter-estimation methods, but the authors highlight the computational cost of simulated annealing. The approach is applied to a finalized DT, so the uncertainty is not accounted for during training.

	\item The \ac{UDT} algorithm~\cite{Tsang2011} extends the probabilistic splits to the training phase, assuming a Gaussian noise, and achieving accuracy gains in $10$ datasets. The method takes an oversampling strategy, where each measurement is replaced by $s$ points. The authors propose optimized searches to control the therefore increased training time. In \ac{UDT}, the same noise model is used for training and evaluation. In medical research and practice, however, the noise in the data used to obtain evidence can be very different from the noise in the data used in practice~\cite{green2014closing}. A learning algorithm ideally supports independent uncertainty models for the training data and the target data.

	\item Soft approaches have also been proposed for multivariate \acp{DT}. An algorithm~\cite{irsoy2012soft} employs logistic regression, treating the separation of data among the children as a classification problem. The method led to small accuracy gains in 10 classification datasets, with significant reductions in tree size. Hierarchical mixtures of experts are another multivariate tree architecture~\cite{jordan1994hierarchical}, where each test is a softmax linear combination of all variables, with performance gains compared to CART. Multivariate \acp{DT} can compactly model complex phenomena, and generally improve accuracy with fewer nodes. They are however intended to generate predictions that do not need to be understood, lacking the interpretability of univariate~\acp{DT}.
\end{itemize}

The aforementioned approaches show the benefit of probabilistic \acp{DT} in improving \ac{DT} performance. However, the strategies used to model noise are either limited to the evaluation phase, or they do not use independent noise distributions for learning and evaluation. Moreover, the behavior of the methods upon increasing levels of noise was not investigated, and no distinction was made between the impact of noise in the training data and noise in the test data. A practical, flexible and well-understood approach for handling measurement uncertainty has not been established for~\ac{DT} learning.

This manuscript proposes a probabilistic \ac{DT} learning approach to handle uncertainty, modeled as a distribution of noise added to the real measurements. Unlike previous methods, the uncertainty distribution is orthogonally considered:
\begin{enumerate}[leftmargin=5mm]
	\item during training, when searching for the best threshold at each node, denoted \acf{SS};
	\item during training, in the propagation of the training data through the \ac{DT}, denoted \acf{STP};
	\item in the propagation of the test instances when evaluating a finalized tree, or \acf{SE}.
\end{enumerate}
\noindent
To promote interpretability, we consider univariate \acp{DT} where each node contains a decision based on a single variable. The proposed \ac{SS} keeps the computational cost under control. We address the problem of integrating knowledge about the uncertainty coming from clinical experience and from the meta-analysis of clinical studies. As a proof-of-concept, we opted for a normal noise model, and evaluate its impact on the distinct learning phases. We also separately study the effect of corrupting the training or test data.

In the following, we introduce the ID3 and C4.5 algorithms and discuss a probabilistic interpretation. The manuscript then proceeds with a description of the proposed soft \ac{DT} algorithm components, followed by the experiments to illustrate and evaluate them.

	\section{Decision tree learning} \label{background}

Consider the input variable~$\mathbf{X}$ with dimensions $X_j$, $j = 1, ..., M$, and the output~$Y$, related by the unknown distribution~$P_{\mathbf{X}Y}(\mathbf{x}, y)$. Drawing samples from $P_{XY}$ composes the training dataset~$\mathcal{D}$. The supervised learning problem consists in learning a model from $\mathcal{D}$ that predicts $Y$ for an unseen sample~$\mathbf{x}$. \Acp{DT}~algorithms follow an algorithmic approach that does not attempt to learn~$P_{\mathbf{X}Y}(\mathbf{x}, y)$.

Learning an optimal \ac{DT} that maximizes generalization accuracy with a minimum number of decisions is NP-complete~\cite{Hyafil1976}. Although non-greedy methods exist for multivariate~\acp{DT}, locally-optimal approaches offer a good trade-off of accuracy and computational complexity. Notable examples include the \ac{ID3} \cite{Quinlan1986}, CART~\cite{Breiman1984} and C4.5~\cite{Quinlan1993}. The \ac{ID3} selects binary splits of numerical variables using the \textit{information gain}, and the C4.5 extends it to categorical features. For a survey on top-down \ac{DT} induction, refer to Rokach~\textit{et~al.}~\cite{Rokach2005}.

\Acp{DT} are regularized by employing node \textit{pruning} algorithms. C4.5 employs a postpruning approach that cuts branches with high estimates of generalization error, based on the training data

\subsection{Split search} \label{split-search}

In top-down \ac{DT} learning, suppose that a new node~$\mathcal{n}$ sees the training subset $\mathcal{D}^{(\mathcal{n})}\subset\mathcal{D}$. In a binary \ac{DT}, we define the branch function at $\mathcal{n}$ as $B^{(\mathcal{n})}(\mathbf{x})$ that indicates if the instance $\mathbf{x}$ goes to the left or right child of~$\mathcal{n}$, $\mathcal{n}_L$ or $\mathcal{n}_R$. $B^{(\mathcal{n})}(\mathbf{x})$ is parametrized by the attribute index $j^{(\mathcal{n})}$ and the threshold~$\tau^{(\mathcal{n})}$:
\begin{equation}
B^{(\mathcal{n})}(\mathbf{x}) =
\left\{
\begin{array}{ll}
\mathcal{n}_L & \mbox{if  } x_{j^{(\mathcal{n})}} < \tau^{(\mathcal{n})} \\
\mathcal{n}_R & \mbox{otherwise},
\end{array}
\right.
\label{eqn:branch}
\end{equation}

\noindent
where $x_{j^{(\mathcal{n})}}$ is the observation of variable $X_{j^{(\mathcal{n})}}$ for~$\mathbf{x}$. The split search consists in finding the attribute $X_j$ and threshold $\tau$ that split $\mathcal{D}^{(\mathcal{n})}$ with maximum class separation.

The entropy of variable~$Y$ is defined as:
\[
H_\alpha(Y) = - \sum_{y} p(y) {log}_\alpha p(y),
\]
\noindent
We take $\alpha=2$ such that the entropy is measured in \textit{bit}, and omit~$\alpha$. In \ac{ID3}, increasing class purity corresponds to reducing $H(Y)$ in $\mathcal{n}_L$ and $\mathcal{n}_R$, compared to~$\mathcal{n}$. In other words, $j^{(\mathcal{n})}$ and $\tau^{(\mathcal{n})}$ are chosen to maximize the \textit{Mutual Information} between $B^{(\mathcal{n})}$ and~$Y$,~$I(B^{(\mathcal{n})};Y)$:
\begin{equation}
j^{(\mathcal{n})}, \tau^{(\mathcal{n})} = \operatorname*{arg\,max}_{j, \tau} I(B^{(\mathcal{n})};Y).
\label{eqn:infoGain}
\end{equation}
\noindent
In the \ac{DT} literature, $I(B^{(\mathcal{n})};Y)$ is known as \textit{information gain}, and is equal to the difference between the entropy of~$Y$ in~$\mathcal{D}^{(\mathcal{n})}$ and the entropy of~$Y$ in the resulting nodes:
\[
I(B^{(\mathcal{n})};Y) = H(Y) - H(Y|B^{(\mathcal{n})}).
\]
\noindent
The term $H(Y|B^{(\mathcal{n})})$ is the conditional entropy of $Y$ given the split variable~$B^{(\mathcal{n})}$. Equation~\ref{eqn:infoGain} is equivalent to:
\begin{align}
j^{(\mathcal{n})},\tau^{(\mathcal{n})}		& = \operatorname*{arg\,min}_{j, \tau} H(Y|B^{(\mathcal{n})}) \nonumber \\
& = \operatorname*{arg\,min}_{j, \tau} \sum_{b \in \{\mathcal{n}_L ,\mathcal{n}_R\}}{p(b)H(Y|B^{(\mathcal{n})}=b)},
\label{eqn:condEnt}
\end{align}
\noindent where
\begin{equation}
H(Y|B^{(\mathcal{n})}=b)	 = - \sum_{y=1}^{K}{ p(y|b)log_2p(y|b) }.
\label{eqn:condEnt2}
\end{equation}
\noindent
The maximum-likelihood estimates of the probabilities in Equations~\ref{eqn:condEnt} and~\ref{eqn:condEnt2} from the training data are:
\begin{gather}
\hat{p}(b) = \tfrac{N_\mathcal{n}(b)}{N_\mathcal{n}} \text{ and } \hat{p}(y|b) = \tfrac{N_\mathcal{n}(y,b)}{N_\mathcal{n}(b)},
\label{eqn:probs}
\end{gather}

\noindent where $N_\mathcal{n}$ is the number of training instances in~$\mathcal{n}$, $N_\mathcal{n} = \sum_{\mathbf{x} \in \mathcal{D}^{(\mathcal{n})}} 1$, and $N_\mathcal{n}(e)$ is the number of instances in $\mathcal{n}$ for which event $e$ was observed, $N_\mathcal{n}(e) = \sum_{\mathbf{x} \in \mathcal{D}^{(\mathcal{n})}} \mathbf{1}(e)$.

In this manuscript, we consider numerical variables, and employ binary splits and use the information gain as a split criterion. In C4.5 release 8, the best split for each numerical attribute is first selected using information gain. Subsequently, the \textit{gain ratio} criterion is used to compare the best splits found for all variables, categorical and numerical~\cite{Quinlanpersonal}.

\footnote{C4.5 performs $n$-ary splits of categorical variables. Since the information gain is biased towards splits with many outcomes, C4.5 normalizes it by the entropy of the $B^{(\mathcal{n})}$ variable, defining the \textit{gain ratio}. In C4.5 release 8, the best split for each numerical attribute is first selected using information gain. Subsequently, the gain ratio is used to compare the best splits found for all variables, categorical and numerical\cite{Quinlanpersonal}.}

\begin{equation}
\tfrac {H(Y|B_j)} {H(B_j)}.
\end{equation}

\subsection{Probabilistic interpretation of the splits} \label{probabilistic}

Consider $p(y|\mathbf{x})$ the probability of $y$ given instance $\mathbf{x}$, based on which we can make a prediction about $y$. We can estimate $p(y|\mathbf{x})$ from the \ac{DT} rooted in node $\mathcal{n}$, with child nodes $\mathcal{n}_L$ and~$\mathcal{n}_R$ as:
\begin{align}
\begin{split}
\hat{p}(y|\mathbf{x}) & = \textstyle \sum_{b \in \{ \mathcal{n}_L, \mathcal{n}_R \}}{ p(y|b, \mathbf{x}) p(b|\mathbf{x}) } \\
& = \textstyle \sum_{b \in \{ \mathcal{n}_L, \mathcal{n}_R \}}{ p(y|b) p(b|\mathbf{x}) } \\
& = p(y|\mathcal{n}_L) p(\mathcal{n}_L|\mathbf{x}) + p(y|\mathcal{n}_R) p(\mathcal{n}_R|\mathbf{x}),
\label{eqn:splits}
\end{split}
\end{align}

\noindent
where $y$ and $\mathbf{x}$ are conditionally independent given~$B^{(n)}=b$. Equation~\ref{eqn:splits} is based on Bayesian model averaging~\cite{hoeting1999bayesian}, which translates the uncertainty of the distinct models, in this case the two subtrees $\mathcal{n}_L$ and $\mathcal{n}_R$, into uncertainty in the class prediction.

Probabilistic \acp{DT} use this idea to instead express the uncertainty about the observed instance~$\mathbf{x}$. This uncertainty is modeled by the posterior $p(\mathcal{n}_L|\mathbf{x})$, known as the \textit{gating function}, $g_\mathcal{n}(\mathbf{x})$~\cite{irsoy2012soft}. The estimate of $p(y|\mathbf{x})$ becomes:
\begin{equation}
\hat{p}(y|\mathbf{x}) = p(y|\mathcal{n}_L) g_\mathcal{n}(\mathbf{x}) + p(y|\mathcal{n}_R) \big(1 - g_n(\mathbf{x})\big).
\label{eqn:prediction}
\end{equation}

If we assume that the observed value is accurate, node $\mathcal{n}$ performs a \textit{hard split} as in Equation~\ref{eqn:branch}, and so:
\[
g_\mathcal{n}(\mathbf{x}) = \mathbf{1}(x_{j^{(\mathcal{n})}} < \tau^{(\mathcal{n})}),
\]
\noindent In this case, if $x_{j^{(\mathcal{n})}}$ is close to $\tau^{(\mathcal{n})}$, small variations in its value can drastically change the estimate $p(y|\mathbf{x})$, and produce incorrect predictions~\cite{Quinlan1987}. Probabilistic \acp{DT} instead model the uncertainty each variable as a distribution of noise, and $g_\mathcal{n}(\mathbf{x})$ becomes the \ac{CDF} of the chosen distribution. A small variation of $x_{j^{(\mathcal{n})}}$ around the threshold value will then smoothly alter~$p(y|\mathbf{x})$

When searching for the best split for a node $\mathcal{n}$, we use the training subset~$\mathcal{D}^{(\mathcal{n})}$ to obtain the probability estimates of Equation \ref{eqn:probs}. E.g. we can estimate $p(\mathcal{n}_L)$, or equivalently $p(\mathcal{n}_R)$, in terms of~$g_\mathcal{n}(\mathbf{x})$:
\begin{align}
\begin{split}
\hat{p}(\mathcal{n}_L) & = \tfrac{ \sum_{\mathbf{x} \in \mathcal{D}^{(\mathcal{n})}} p(\mathcal{n}_L|\mathbf{x})}{\sum_{\mathbf{x} \in \mathcal{D}^{(\mathcal{n})}}1
= \tfrac{ \sum_{\mathbf{x} \in \mathcal{D}^{(\mathcal{n})}} g_\mathcal{n}(\mathbf{x})}{\sum_{\mathbf{x} \in \mathcal{D}^{(\mathcal{n})}}1
\end{split}
\label{eqn:splits-train}
\end{align}

\section{Proposed approach} \label{proposed-approach}

We propose a probabilistic \ac{DT} approach to handle uncertainty by modeling it as a distribution of noise around the observed value. The noise model should expresses existing knowledge about the uncertainty. This model is independently considered:
\begin{enumerate}
	\item When searching for the split thresholds, during training, or  (Section ~\ref{soft-search}, \acf{SS}),
	\item When propagating the training instances, during training (Section ~\ref{soft-propagation-training}, \acf{STP}),
	\item When propagating test instances through the constructed~\ac{DT}, during evaluation \ac{DT} (Section ~\ref{soft-evaluation}, \acf{SE}).
\end{enumerate}
As a proof-of-concept, we consider that the noise of variable~$X$ is additive and captured by the variable $\mathcal{E}_{X} \sim \mathcal{N}(0, \sigma^2)$. We consider the normal distribution to be a good assumption to study our approach. It can be useful for various types of data, owing to the Central Limit Theorem, and is fully described by its mean and~variance.

\subsection{Soft search} \label{soft-search}

Let us focus on the computation of the information gain for variable~$X$. Suppose that $X$ takes $D$ \textit{distinct} values sorted as $\{ x^{(1)}, ..., x^{(D)} \}$ with $x^{(i)} < x^{(i+1)}$, $i=1, ..., D-1$ in the training data $\mathcal{D}$. Note that $|\mathcal{D}| \ge D$. We focus on finding a split for node~$\mathcal{n}$, and consider that all the concepts in this Section refer to~$\mathcal{n}$, and omit the subscript. $N(x^{(i)})$~is the number of instances with value $x^{(i)}$, and $N(y, x^{(i)})$ the number of such instances that belong to class~$y$. Let $\tau$ denote the candidate threshold for a split based on $X$ at node~$\mathcal{n}$.

In C4.5, the split search is done by computing the information gain in Equation~\ref{eqn:infoGain} for each candidate $\tau = x^{(i)}$, $i=2, ..., D$. Since the method assumes certain measurements, we have $g(\mathbf{x}) = \mathbf{1}(x<\tau)$. The probabilities $p(b)$ and $p(y|b)$ are estimated as in Equation~\ref{eqn:probs}. We now describe how to estimate these probabilities efficiently, considering uncertain measurements.

If we consider the noise model $\mathcal{E}_{X} \sim \mathcal{N}(0, \sigma^2)$, the gating function will be the normal~\ac{CDF}. Let us denote the numerator of Equation~\ref{eqn:splits-train} as $\rho(\tau)$, representing the density of training examples falling on the left child node:
\begin{align}
\begin{split}
\rho(\tau) = \sum_{\mathbf{x} \in \mathcal{D}} p(\mathcal{n}_L|\mathbf{x}) =
\sum_{\mathbf{x} \in \mathcal{D}} g_\mathcal{n}(\mathbf{x})
=  \sum_{\mathbf{x} \in \mathcal{D}} \int_{-\infty}^{\tau} K(\tau - x, \sigma) d\tau,
\end{split}
\end{align}
\noindent where $K(\tau - x, \sigma)$ is the Gaussian kernel function centered on the measurement $x$ with variance~$\sigma^2$:
\begin{equation}
K(\tau - x, \sigma) \propto exp(-\tfrac{(\tau - x)^2}{2\sigma^2}).
\end{equation}
\noindent This corresponds to using Gaussian kernel density estimation for the probabilities used to compute the information gain. Since $\rho(\tau)$ is no longer constant between each $x^{(i)}$ and $x^{(i+1)}$, the candidate thresholds need not be the dataset values. We take:
\begin{equation}
\tau = \tau_{min}, ~ \tau_{min} + \delta \sigma, ~ \tau_{min} + 2 \delta \sigma, ~ ..., ~ \tau_{max}
\label{eqn:discretization}
\end{equation}
\noindent
with $\tau_{min} =  x^{(1)}-\tfrac{\omega}{2}\sigma$ and $\tau_{max} \le x^{(D)}+\tfrac{\omega}{2}\sigma$. The parameter $\delta$ controls the resolution of the search, and $\omega$ is the \textit{window factor}, $\omega>\delta$. For efficiency, we consider the kernel adjusted to be zero for~$|\tau - x| > \omega$. Using Equation~\ref{eqn:discretization}, the number of information gain computations is limited to $\tfrac{1}{\delta} ( \tau_{max} - \tau_{min} )$.

To compute the information gain efficiently for the candidate splits, we keep two running sums of $\rho(\tau)$, on the left and on the right of~$\tau$,  $\rho_L(\tau)$ and  $\rho_R(\tau)$. We initialize the search with  $\rho_L(\tau) = 0$ and  $\rho_R(\tau)=N$. Each time $\tau$ is incremented $\delta\sigma$, the left and right sums are respectively incremented and decremented the quantity $\Delta\rho(\tau)$:
\begin{equation}
\Delta\rho(\tau) =
\begin{cases}
\mathlarger{\sum}_{i=1}^{D}{N(x^{(i)}) \Phi \big(\tfrac{\tau - x^{(i)}}{\sigma} \big) },
\hfill \mbox{ ~~if } \tau=\tau_{min} \\
\mathlarger{\sum}_{i=1}^{D}{N(x^{(i)}) \Big[ \Phi \big(\tfrac{\tau - x^{(i)}}{\sigma} \big) - \Phi \big(\tfrac{\tau - \delta \sigma - x^{(i)} }{\sigma} \big)\Big] }, 	\\
\hfill  \mbox{ ~~if } \tau_{min} < \tau < \tau_{max} \\
\mathlarger{\sum}_{i=1}^{D}{N(x^{(i)}) \Big[ 1 - \Phi \big(\tfrac{\tau - \delta \sigma - x^{(i)} }{\sigma} \big)\Big] },
\hfill  \mbox{ ~~if } \tau=\tau_{max}
\end{cases}
\label{eqn:filter}
\end{equation}
\noindent with $\Phi(\cdot)$ the \ac{CDF} of the standard normal distribution. The contribution of measurement $x^{(i)}$ to $\Delta\rho(\tau)$ is proportional to $\Phi \big(\tfrac{\tau - x^{(i)}}{\sigma} \big) - \Phi \big(\tfrac{\tau - x^{(i)} - \delta \sigma}{\sigma} \big)$. The last point $\tau_{max}$ ensures that the density contributions of $x^{(i)}$ sum to~$N(x^{(i)})$.

Similarly, to estimate $\hat{p}(y|{\mathcal{n}}_L)$, we consider the sum of the densities per class, $\rho(y,\tau)$. The density increments per class $\Delta\rho(y,\tau)$ are computed by replacing the number of instances $N(x^{(i)})$ by $N(y, x^{(i)})$ in Equation~\ref{eqn:filter}. Finally, the estimated probabilities $\hat{p}({\mathcal{n}}_L)$ and $\hat{p}(y|{\mathcal{n}}_L)$ are used to minimize the conditional entropy $H(Y|B)$ in Equation~\ref{eqn:condEnt}. Searching for the threshold using the set of values in Equation~\ref{eqn:discretization} and the density increments $\Delta\rho(\tau)$ and $\Delta\rho(y,\tau)$ acts as a Gaussian filter to the information gain. We refer to this approach as \acf{SS}

The standard deviation of $\mathcal{E}_{X}$ is set to $\sigma=u_{s} \bar x$, where $u_{s}$ is the \ac{SS} \textit{uncertainty factor} and $\bar x$ is the sample mean of~$X$. This choice of $\sigma$ was made, as clinicians often report the uncertainty in terms of the absolute values of the attributes~\cite{d2016two}. In the following experiments, $u_s$ is assumed to be the same for all variables, although it could be specified independently.

\subsection{Soft training propagation} \label{soft-propagation-training}

We can also account for the uncertainty when splitting the training data, which we call \acf{STP}. The two soft training approaches, \ac{SS} and \ac{STP} can be used in combination or independently.

Consider the split variable $X$ and value~$\tau^{(\mathcal{n})}$ at node~$\mathcal{n}$. As before, a soft split is achieved by setting the gating function to the \ac{CDF} of~$\mathcal{E}_{X}$ centered in $x$:
\begin{equation}
g_{\cal{n}}(\mathbf{x}) = F_{\mathcal{E}_{X}}(\tau - x) = \Phi \big(\tfrac{\tau - x}{\sigma} \big).
\label{eqn:soft-splits}
\end{equation}
\noindent with $F_{\mathcal{E}_{X}}(.)$ the \ac{CDF} of~$\mathcal{E}_{X}$. Each training instance is fractionally divided between the child nodes, according to the probability~$g_{\cal{n}}(\mathbf{x})$. The standard deviation is set to $\sigma=u_{p} \bar x$, with $u_{p}$ the \ac{STP} uncertainty factor, and $\bar x$ the variable mean.

\Ac{STP} leads to increased learning times. Using hard splits, each instance is sent down one branch, and the total number of instances remains constant at each level of the tree. The number of information gain computations at any \ac{DT} depth is bounded by the dataset size and the number of attributes, $|\mathcal{D}| M$. With \ac{STP}, each instance is sent down all the branches, with weights determined by the gating function. Therefore, a node sees a \textit{greater} number of instances compared to the hard approach. The number of information gain computations at a given \ac{DT} depth $d$ is raised to~$2^{d} |\mathcal{D}| M$.

\subsection{Soft evaluation} \label{soft-evaluation}

The uncertainty may also be accounted for when classifying test cases with a finalized~\ac{DT}. \Acf{SE} is achieved by setting the gating function as in Equation~\ref{eqn:soft-splits}, but with $\sigma=u_{e} \bar x$. The evaluation uncertainty factor is denoted as~$u_e$ and $\bar x$ is obtained from the \textit{training} data. From Equation~\ref{eqn:prediction}, the class probability estimates for instance~$\mathbf{x}$ become:
\begin{equation}
\hat{p}(y|\mathbf{x}) = p(y| \mathcal{n}_L)\Phi(\tfrac{\tau-x}{\sigma}) + p(y|\mathcal{n}_R)\big(1 - \Phi(\tfrac{\tau-x}{\sigma}) \big).
\label{eqn:comp2}
\end{equation}

\noindent As described in Section~\ref{probabilistic}, this approach adjusts $\hat{p}(y|\mathbf{x})$ to reflect the choice of noise distribution, in this case, Gaussian.

\subsection{Motivation on a toy example} \label{motivation}

\noindent To better motivate the use of the soft algorithm components, we consider a toy example with a input variable $X$ and class~$Y \in \{0, 1\}$.

\subsubsection*{\Acl{SS}}

\begin{figure*}	\centering
	\begin{subfigure}[b]{\textwidth} \centering
		\includegraphics[width=\textwidth]{4point.eps} \caption{Example with $(x_{(1)}, y_{(1)}) = (x_{(2)}, y_{(2)})=(-2,0)$ and $(x_{(3)}, y_{(3)})=(x_{(4)}, y_{(4)})=(2, 1)$.\\} \label{fig:4point-a}
	\end{subfigure}
	\begin{subfigure}[b]{\textwidth} \centering
		\includegraphics[width=\textwidth]{4point-noise-2.eps} \caption{Example with $(x_{(1)}, y_{(1)}) = (-2,0)$ and $(x_{(4)}, y_{(4)})=(2, 1)$, and noisy measurements $(x_{(2)}, y_{(2)})=(0.2,0)$ and $(x_{(3)}, y_{(3)})= (-0.2, 1)$.} \label{fig:4point-b}
	\end{subfigure}
	\caption{Motivating example to show the effect of employing the \acf{SS} on the information gain $I(X;Y)$, with $\sigma_s=0.3$.}
	\label{fig:4point}
\end{figure*}

Suppose we have a training dataset of four examples $\{ (x_{(1)} ,0),  (x_{(2)} ,0), (x_{(3)} ,1), (x_{(4)} ,1)\}$, with $x_{(1)} = x_{(2)} = -2$ and $x_{(3)} = x_{(4)} = 2$. Figure~\ref{fig:4point-a} displays the information gain $I(X;Y)$ computed by C4.5 and the \ac{SS} for the corresponding candidate splits, denoted by~$\tau$. C4.5 would choose $\tau = 2$ as a split, while \ac{SS} would choose a split close to $\tau = 0$, potentially avoiding the misclassification of test examples in the interval $(0,2)$, for which there is no training data. If the measurements are noisy such that e.g. $x_{(2)}=0.2$ and $x_{(3)}=-0.2$, C4.5 would select $\tau=0.2$ or $\tau=2$, as seen in Figure~\ref{fig:4point-b}. The \ac{SS} would select a value close to zero.

\subsubsection*{\Acl{STP}}

Let us now discard two training examples, and keep $\mathbf{x}_{(1)} = (-2,0)$ and $\mathbf{x}_{(4)} = (2,1)$. We model the uncertainty of $X$, such that $x_{(1)} \sim \mathcal{N}(-2, \sigma^2)$ and $x_{(4)} \sim \mathcal{N}(2, \sigma^2)$ are drawn from normal distributions with mean $\mu_{(1)} = -2$ and $\mu_{(2)} = 2$, and common variance $\sigma^2$. Let us also consider the test point $\mathbf{x}_{(t)}=(x_{(t)}, 1)$. We analyze the probability of misclassifying $\mathbf{x}_{(t)}$. We note that the sum and difference of normal distributions are also normally-distributed:
\begin{gather*}
d = x_{(1)} - x_{(4)} \sim \mathcal{N}(-2 -2, \sigma^2 + \sigma^2 ) = \mathcal{N}(-4, 2\sigma^2 ),\\
s = x_{(1)} + x_{(4)} \sim \mathcal{N}(-2 + 2, \sigma^2 + \sigma^2 ) = \mathcal{N}(0, 2\sigma^2 ), \\
m = \tfrac{s}{2} =\tfrac{x_{(1)} + x_{(4)}}{2} \sim \mathcal{N}(\tfrac{1}{2} \cdot 0, \tfrac{1}{2^2}) = \mathcal{N}(0, \tfrac{\sigma^2}{2})
\end{gather*}
\noindent Let us assume that the search selects the midpoint $ m$ as a split value. The prediction $\hat{y}(x_{(t)})$ depends on whether the difference $c=x_{(t)}-m$ is greater than or less than~$0$. Specifically, the probability of misclassifying $\mathbf{x}_{(t)}$ is composed by two terms:

\begin{figure}[h]
	\captionsetup[subfigure]{justification=centering} \centering
	\begin{subfigure}{.49\textwidth} \centering
		\includegraphics[width=.95\textwidth]{perror-1.eps} \vspace{-4mm}
		\caption{ $P\big(\hat{y}(x_{(t)})$ when we assume $x_{(t)}$ to be certain.} \label{fig:perror-ex1a}
	\end{subfigure}
	\begin{subfigure}{.49\textwidth} \centering
		\includegraphics[width=.95\textwidth]{perror-2.eps} \vspace{-4mm}
		\caption{ $P\big(\hat{y}(x_{(t)})$ when $x_{(t)}$ is assumed to have noise.} \label{fig:perror-ex1b}
	\end{subfigure}
	\caption{Probability of misclassifying $(x_{(t)}, 1)$ as function of the standard deviation $\sigma$ of the normal uncertainty model. In \ref{fig:perror-ex1a}, the uncertainty model is considered only for the training instances $x_{(1)}$ and $x_{(4)}$, simulating \acf{STP}, while $x_{(t)}$ is certain. In \ref{fig:perror-ex1b} $x_{(t)}$ has normally-distributed noise, as in \acf{SE}.}
	\label{fig:perror-ex1}
\end{figure}

\begin{gather*}
P\big(\hat{y}(x_{(t)}) = 0\big) = P (d < 0)P (c < 0) + P (d > 0)P (c > 0).
\end{gather*}
\noindent In case $x_{(1)}$ is smaller than~$x_{(4)}$, C4.5 assigns class $0$ to $x_{(t)}$ iff~$x_{(t)} < m$. In case $x_{(1)}$ is greater than~$x_{(4)}$, C4.5 assigns class $0$ to $x_{(t)}$ iff~$m < x_{(t)}$. Each of the above probabilities is given by the \ac{CDF} associated with the corresponding distribution. If we shift each distribution to have $0$ mean, we can express the probabilities as:
\begin{gather*}
P (d < 0) = P (d + 4 < 4) = F_{d+4}(4) = \Phi\big(\tfrac{4}{2 \sigma}\big), \\
P (d > 0) = 1 - P (d < 0) = 1 - \Phi\big(\tfrac{4}{2 \sigma}\big).
\end{gather*}
$F_{d+4}$ is the \ac{CDF} of $d+4$. Figure~\ref{fig:perror-ex1a} displays the probability of misclassifying $x_{(t)}$ as a function of $\sigma$ when $x_{(t)}$ is certain. When $x_{(t)} = 1$, and so $c > 0$, $P\big(\hat{y}(x_{(t)}) = 0\big)$ is nearly zero until the distributions of the training instances start to overlap with increasing~$\sigma$. The inverse occurs for $c < 0$, where the error probability starts to decrease. Figure~\ref{fig:perror-ex1a} shows how the prediction probabilities change as the model expresses less confidence on the training data, when using~\ac{STP}.

\subsubsection*{\Acl{SE}}

Let us now express the uncertainty about the test example, such that $x_{(t)}$ is also normally distributed with variance~$\sigma^2$. We can obtain $P (c < 0)$ and $P (c > 0)$ as:
\begin{align*}
\begin{split}
P (c < 0) = P ( c - 1 < -1 ) = F_{c-1}(-1) & = 1 - F_{c-1}(1) = \\
& = 1 - \Phi\big(\tfrac{\sqrt{2}}{\sqrt{3} \sigma}\big), \\
P (c > 0) = 1 - P ( c < 0 ) = \Phi\big(\tfrac{\sqrt{2}}{\sqrt{3} \sigma}\big).
\end{split}
\end{align*}

\noindent Figure~\ref{fig:perror-ex1b} displays $P\big(\hat{y}(x_{(t)}=0)$ for $x_{(t)} \sim \mathcal{N}(1, \sigma^2)$ and $x_{(t)} \sim \mathcal{N}(-1, \sigma^22)$, which now converges faster to~$0.5$.

\section{Experiment design}

The experiments assess effect of the proposed uncertainty model in the distinct \ac{DT} learning phases, and the impact of corrupting the training versus the test data with noise.

Improving prediction performance consists in maximizing the generalization accuracy, i.e. the fraction of correctly classified test examples, while minimizing model complexity~\cite{Kearns1997}. \Ac{DT} complexity is assessed by measuring its number of leaves.

Each proposed soft component, \ac{SS}, \ac{STP} and \ac{SE} is independently compared to C4.5. The C4.5 pruning and missing-value strategies are equally employed in all experiments. Pruning is extended with the Laplace correction, as recommended to favor the exploration of smaller \acp{DT} for the same accuracy~\cite{Niblett1986}. We also extend the evaluation of the \ac{PLT}~\cite{quinlan1990probabilistic} and \ac{UDT}~\cite{Tsang2011} approaches to more data and noise scenarios.

\subsection{Data description and pruning confidence factor} \label{data}

Table~\ref{tab:datasets} displays the employed datasets. We sought clinical data from open resources, including UCI \ac{ML}~\cite{Lichman2013} and KEEL~\cite{Alcala-Fdez2011} repositories. Non-ordinal features were excluded. We synthesized $5$ additional datasets using an adaptation of the method by Guyon \cite{Guyon2003}, available in Scikit-learn's implementation \textit{make\_classification}~\cite{scikit-learn}. The datasets were generated with varying properties.

Since optimal \ac{DT} size is problem- and dataset-dependent~\cite{jensen1999effects}, the pruning confidence factor is not optimized. Instead, for each dataset the confidence factor is fixed such that the average \ac{DT} size achieved by C4.5 through \ac{CV} is 15 leaves. This corresponds to a binary tree wit nearly $4$ levels, considered a manageable/interpretable number of decisions in a visual clinical guideline. Some of the datasets were too small to reach 15 leaves, so their confidence factor is set to either 10 or 5 leaves. The confidence factor is fixed across all experiments.

\subsection{Experiments} \label{model-evaluation}

For each real dataset, $30$ random train-test permutations were created, containing respectively 70\% and 30\% of the data. Stratified sampling was used, so that class proportions are equal in all samples. For each synthetic dataset, distinct instances were generated and divided into 30 different sets, which were then split into 70\%-30\% train-test samples.

The experiments are performed with varying degrees of noise to the data. The noise added to a variable $X$ in a data subset is sampled according to $\mathcal{N}(0,n \bar x)$, with $n$ the \textit{noise factor} and $\bar x$ the training subset mean of $X$. The same $n$ is used for all its variables. All randomness was generated with fixed seeds for reproducibility.

\begin{table}
	\renewcommand{\arraystretch}{0.97} \centering
	\caption{ Classification datasets used in the experiments. }
	\label{tab:datasets}
	\setlength{\tabcolsep}{3.2pt}
	\begin{tabular}{llccc}
		\toprule
		\multirow{2}{*}{\specialcell{Dataset}} & \multirow{2}{*}{\specialcell{Source}}								& \multirow{2}{*}{\specialcell{No.\\variables}}	& \multirow{2}{*}{\specialcell{No.\\instances}}			& \multirow{2}{*}{\specialcell{No.\\classes}} \\ \\
		\midrule \midrule
		Hepatitis							& KEEL				& 6			& 155				& 2		\\
		Heart disease						& UCI \ac{ML}		& 8			& 303				& 2		\\
		Pima Indians diabetes				& UCI \ac{ML}		& 8			& 768				& 2		\\
		South African heart					& UCI \ac{ML}		& 8			& 462				& 2		\\
		Breast cancer Wisconsin				& UCI \ac{ML}		& 39		& 569				& 2		\\
		Dermatology							& UCI \ac{ML}		& 34		& 366				& 6		\\
		Haberman's breast cancer			& UCI \ac{ML}		& 3			& 306				& 2		\\
		Indian liver patient 				& UCI \ac{ML}		& 9			& 582				& 2		\\
		BUPA liver disorders				& UCI \ac{ML}		& 6			& 345				& 2		\\
		Vertebral column (2 classes)		& UCI \ac{ML}		& 12		& 310				& 2		\\
		Vertebral column (3 classes)		& UCI \ac{ML}		& 12		& 310				& 3		\\
		Thyroid gland						& UCI \ac{ML}		& 5			& 215				& 3		\\
		Oxford Parkinson's disease 			& UCI \ac{ML}		& 22		& 194				& 2		\\
		SPECTF								& UCI \ac{ML}		& 44		& 266				& 2		\\
		Thoracic surgery					& UCI \ac{ML}		& 3			& 470				& 2		\\
		Synthetic 1							& Guyon							& 15			& 30$\times$500		& 2		\\
		Synthetic 2							& Guyon							& 15			& 30$\times$400		& 2		\\
		Synthetic 3							& Guyon							& 20			& 30$\times$300		& 2		\\
		Synthetic 4							& Guyon							& 25			& 30$\times$200		& 3		\\
		Synthetic 5							& Guyon							& 20			& 30$\times$250		& 3		\\
		\bottomrule
	\end{tabular}
\end{table}

\subsubsection{Experiment 1: noise added to the training data}\label{experiment-1}

Experiment~$1$ studies the approaches with noise added to the training data, so we denote the training noise factor as~$n_{train}$. The uncertainty factors $u_s$, $u_t$, $u_e$ of the \ac{SS}, \ac{STP} and \ac{SE} approaches, as well as the \ac{UDT} parameter $w$, control the standard deviation of the uncertainty model. As such, we tune them by \acf{CV} for each dataset and~$n_{train}$. \noindent The \ac{SS} parameters $\omega$ and $\delta$ are respectively set to $6.0$ and $0.1$. Initial experiments showed that they do not significantly impact the results, provided that $\omega$ is large enough to contain most of the density of the uncertainty distribution, and $\delta$ is small enough to ensure a large set of candidate splits. The \ac{PLT} parameters $\tau^-$ and $\tau^+$ are derived using the method proposed by Quinlan~\cite{quinlan1990probabilistic}. We summarize the steps for model tuning and evaluation:
\begin{quote}
	\begin{small}
		For each $70\%-30\%$ train-test permutation, $n_{train}$ and model: \\
		1. Hold out the 30\% test set \\
		2. If model has parameter to set ($u_s$, $u_t$, $u_e$, or $w$), use the 70\% training set to tune it by \ac{CV}: \\
		\text{~~~~~}(a) Compute 10 stratified \ac{CV} folds. \\
		\text{~~~~~}(b) Add noise to the \ac{CV} training folds, distributed \\ \text{~~~~~~~~~} as $\mathcal{N}(0,n_{train} \bar x)$. Do not add noise to the \ac{CV} \\ \text{~~~~~~~~~} validation folds. \\
		\text{~~~~~}(c) Tune the parameter to maximize \ac{CV} accuracy.\\

		3. Add noise to the initial 70\% training set, distributed as $\mathcal{N}(0,n_{train} \bar x)$. \\
		4. Learn a tree using the noisy training set, and the selected parameter value, if applicable. \\
		5. Evaluate the tree on the 30\% test set, to which no noise was added.
	\end{small}
\end{quote}

\noindent C4.5 and \ac{PLT} do not undergo step $2$.

\subsubsection{Experiment 2: noise added to the test data}\label{experiment-2}

Experiment~$2$ evaluates the merit of the models built on data without added noise in predicting the labels of noisy test cases. Noise is added to the data used for evaluation, and we denote $n$ as~$n_{test}$. Model evaluation is done as:
\begin{quote}\begin{small}
		For each $70\%-30\%$ train-test permutation, $n_{test}$ and model: \\
		1. Hold out the 30\% test set \\
		2. If model has parameter to set ($u_s$, $u_t$, $u_e$, or $w$), use the 70\% training set to tune it by \ac{CV}: \\
		\text{~~~~~}(a) Compute 10 stratified \ac{CV} folds. \\
		\text{~~~~~}(b) Do not add noise to the \ac{CV} training folds. \\ \text{~~~~~~~~~} Add noise to the \ac{CV} validation folds, distri-\\ \text{~~~~~~~~~~}buted as $\mathcal{N}(0,n_{test} \bar x)$. \\
		\text{~~~~~}(c) Tune the parameter to maximize \ac{CV} accuracy.\\
		3. Learn a tree using the initial 70\% training set, and the selected parameter value, if applicable \\
		4. Add noise to the 30\% test set, distributed as $\mathcal{N}(0,n_{test} \bar x)$. \\
		5. Evaluate the tree on noisy the 30\% test set.
\end{small}\end{quote}

\noindent As before, C4.5 and \ac{PLT} do not undergo step 2. Note that noise is added to the \ac{CV} validation folds.

\section{Results} \label{results}

Section~\ref{ef-example} illustrates \ac{SS}, \ac{STP} and \ac{SE} on a single variable. In Section~\ref{exp-res}, we show the results of Experiments~$1$ and $2$.

\subsection{Illustration on a single variable} \label{ef-example}

We take the \acf{EF} estimates of the \emph{Data Science Bowl Cardiac Challenge}~\cite{kaggle2015}. \Ac{EF} is a variable of critical importance in cardiology. Implantable device therapy is officially recommended for $EF < 35\%$~\cite{Ponikowski2016}. Therefore, we take cases with $EF < 35\%$ to have positive eligibility, i.e. $Y= \mathbf{1}(EF < 35\%)$, as shown in Figure~\ref{fig:ef1}. Adding random noise to these data results in $7$ \acf{FN} and $5$ \acf{FP}, as seen in Figure~\ref{fig:ef2}.\\

\noindent \textit{\Acl{SS}:} In Section~\ref{motivation}, we motivated the \ac{SS} as way of increasing the set of candidate splits and smoothing the information gain. We now observe this on real measurements. Figure~\ref{fig:ef-disc-a} displays the number of patients for each class and \ac{EF} value, $N(y, x^{(i)})$, like a histogram. Figure~\ref{fig:ef-disc-b} shows the same data with noise. Figure~\ref{fig:ef-disc-c} shows the \ac{SS} density increments $\Delta\rho(y, \tau)$, and \ac{SS} information gain.

\begin{figure}[t]
	\captionsetup[subfigure]{justification=centering}
	\centering
	\begin{subfigure}{.48\textwidth} \centering
		\includegraphics[height=4.3cm]{figure1a.eps}
		\caption{\Ac{EF} measurements.} \label{fig:ef1}
	\end{subfigure}\\
	\begin{subfigure}{.48\textwidth} \centering
		\includegraphics[height=4.3cm]{figure1b.eps}
		\caption{\Ac{EF} measurements with added noise.} \label{fig:ef2}
	\end{subfigure}
	\caption{ \footnotesize (a) \Acf{EF} data of the Data Science Bowl Cardiac Challenge. (b) Same data with noise sampled from $\mathcal{N}(0,0.1\bar{x})$, with $\bar{x}$ the mean. \Ac{EF} $< 35\%$ indicates therapy eligibility leading to \acf{TN}, \acf{TP}, \acf{FP}, \acf{FN}.} \label{fig:ef}
\end{figure}
\begin{figure}[ht]\centering
	\vspace{-3mm}
	\begin{subfigure}{.48\textwidth} \centering
		\includegraphics[trim={0cm 11.5cm 0cm 0cm}, clip, width=.95\textwidth]{figure3-2.eps} \caption{\scriptsize} \label{fig:ef-disc-a}
	\end{subfigure}
	\begin{subfigure}{.48\textwidth} \centering
		\includegraphics[trim={0cm 6.2cm 0cm 5.3cm}, clip, width=.95\textwidth]{figure3-2.eps} \caption{\scriptsize} \label{fig:ef-disc-b}
	\end{subfigure}
	\begin{subfigure}{.48\textwidth} \centering
		\includegraphics[trim={0cm .01cm 0cm 10.7cm}, clip, width=.95\textwidth]{figure3-2.eps} \vspace{-2mm} \caption{\scriptsize} \label{fig:ef-disc-c}
	\end{subfigure}
	\caption{ \footnotesize Information gain computation the \acf{EF} measurements of Figure~\ref{fig:ef}, using the standard search (a,b) or \acf{SS} with $u_s=0.1$ (c). The left (a) and (b) plots show the number of patients for each class and~$x^{(i)}$, $N(y=0, x^{(i)})$ and $N(y=1, x^{(i)})$. The left~(c) plot displays the \ac{SS} density increments, $\Delta\rho(y, \tau)$. The bottom plots show the corresponding information gain.}
	\label{fig:ef-disc}
\end{figure}

\clearpage

Soft search methods such \ac{SS} or \ac{UDT} increase the set of candidate splits. This dataset has $500$ instances. Employing \ac{UDT} with a resampling factor $s=100$ raises the number of information gain computations from $5e2$ to $5e4$~\cite{Tsang2011}. Setting the \ac{SS} parameters e.g. as $\cal{w}=8$, $\delta=0.05$ and $u_s=0.1$ limits this number to a maximum of approximately~$2.3e3$.\\

\noindent \textit{\Acl{STP}:} To illustrate how \ac{STP} alters the class probability estimates, we consider the noisy data of Figure~\ref{fig:ef2}. We learned two single-node two-leaf trees using the standard training propagation and \ac{STP} with $u_{p} = 0.1$. Table~\ref{tab:tab1} shows that the class probability estimates are less extreme when \ac{STP} is employed, as they reflect the choice of noise distribution.\\

\begin{table}[t]
	\centering
	\caption{ Class probability estimates at leaves ${\mathcal{n}}_L$ and ${\mathcal{n}}_R$ of the trees learned on the noisy \acf{EF} dataset of Figure~\ref{fig:ef2}. \Acfp{DT} learned using standard or \acf{STP}.}
	\label{tab:tab1}
	\renewcommand{\arraystretch}{.9}
	\setlength{\tabcolsep}{6pt}
	\begin{tabular}{lcccc}
		\toprule
		& Standard training propagation & \Ac{STP} \\
		\midrule
		$\hat{p}(y=0|\mathcal{n}_L)$ & 0.000 & 0.310	\\
		$\hat{p}(y=1|\mathcal{n}_L)$ & 1.000 & 0.790	\\
		\midrule
		$\hat{p}(y=0|\mathcal{n}_R)$ & 0.981 & 0.979	\\
		$\hat{p}(y=1|\mathcal{n}_R)$ & 0.018 & 0.021	\\
		\bottomrule
	\end{tabular}
\end{table}
\begin{table}[t]
	\centering
	\caption{ Class probability estimates for the misclassified examples of the noisy \acf{EF} data in Figure~\ref{fig:ef2}, estimated by the tree learned with C4.5 on the non-noisy data in Figure~\ref{fig:ef1}.}
	\label{tab:stest}
	\renewcommand{\arraystretch}{.9}
	\setlength{\tabcolsep}{1.5pt}
	\begin{tabular}{llcccc}
		\toprule
		& & \multicolumn{2} {c} {Hard evaluation} & \multicolumn{2} {c}{\Ac{SE}}\\
		\midrule
		& Patient & $\hat{p}(y=0|\mathbf{x})$ & $\hat{p}(y=1|\mathbf{x})$ & $\hat{p}(y=0|\mathbf{x})$ & $\hat{p}(y=1|\mathbf{x})$ \\
		\midrule
		\multirow{6}{*}{False negative}		& 6		& 1.00 & 0.00 & 0.54		& 0.46 \\
		& 33		& 1.00 & 0.00 & 0.83		& 0.17 \\
		& 269	& 1.00 & 0.00 & 0.67		& 0.33 \\
		& 280	& 1.00 & 0.00 & 0.91		& 0.09 \\
		& 295	& 1.00 & 0.00 & 0.51		& 0.49 \\
		& 309	& 1.00 & 0.00 & 0.79		& 0.21 \\
		\midrule
		\multirow{6}{*}{False positive}		& 42		& 0.00 & 1.00 & 0.16		& 0.84 \\
		& 116	& 0.00 & 1.00 & 0.33		& 0.67 \\
		& 167	& 0.00 & 1.00 & 0.39		& 0.61 \\
		& 228	& 0.00 & 1.00 & 0.30		& 0.70 \\
		& 359	& 0.00 & 1.00 & 0.49		& 0.51 \\
		& 478	& 0.00 & 1.00 & 0.30		& 0.70 \\
		\bottomrule
	\end{tabular}
\end{table}

\noindent \textit{\Acl{SE}:} Figure~\ref{fig:ef2} shows that 7 \ac{FN} and 5 \ac{FP} were introduced by the noise added to the \ac{EF} dataset. Table~\ref{tab:stest} shows the probability estimates of those misclassified examples, obtained using hard or soft evaluation, with $u_e = 0.1$. The numbers of \ac{FN} and \ac{FP} are different because the algorithm learned a threshold of $35.37\%$ rather than $35\%$.

\subsection{Experimental results} \label{exp-res}

The average number of leaves, test accuracy and running time are computed over 30 train-test permutations, for each experiment and dataset. The absolute difference to C4.5, averaged over all datasets, is displayed in Tables~\ref{table:res-SS-prop-eval-train-noise} and~\ref{table:res-SS-prop-eval-test-noise}.

However, since absolute results of distinct datasets cannot be directly compared~\cite{Demsar2006}, we focus on standardized metrics. The results of each dataset and method were standardized by the dataset's baseline. The baseline result is obtained by C4.5 without added noise, and estimated with the 30 permutations. The standardization consists in subtracting the baseline mean, and dividing by the baseline standard deviation. E.g. the baseline of the Heart disease dataset has mean $15.0$ leaves and standard deviation $3.0$ leaves. In this case, standardized results $0.0$, $-1.0$ and $1.5$ translate into $15.0$, $12.0$ and $19.5$ leaves, respectively. Tables~\ref{table:res-SS-prop-eval-train-noise} and~\ref{table:res-SS-prop-eval-test-noise} display the standardized metrics, averaged over all datasets, and Figure~\ref{fig:res-box-plots} shows the corresponding boxplots. Computational times are merely indicative, as the experiments were run on a cluster, and the specifications of the machines may vary slightly.

All approaches show an increase in the number of leaves and a reduction in test accuracy as the noise factor grows. The decrease in accuracy was sharper for the predictions made on noisy test data, as seen in Table~\ref{table:res-SS-prop-eval-test-noise}.

\Ac{SS}, \ac{STP} and \ac{UDT} show maintenance or non-statistically significant improvement in accuracy compared to C4.5, in all noise scenarios. In Table~\ref{table:res-SS-prop-eval-train-noise}, we see that \ac{STP} had higher accuracy compared to the other methods for all $n_{train}$ and~$n_{test}$. For the \ac{SS} and \ac{STP} approaches, the maintenance of accuracy was accompanied by statistically significant reductions in the number of leaves compared to C4.5 and \ac{UDT}. The \ac{SS} tree size reduction was statistically significant for noise factors greater than $0.00$. \Ac{STP} had a further reduction in tree size, significant for all $n_{train}$ and~$n_{test}$. The maintenance of accuracy by \ac{UDT} compared to C4.5 was accompanied by an increase in the number of leaves. The method was considerably slower than \ac{SS} and~\ac{STP}.

In Experiment~1, the accuracies obtained by \ac{SE} and \ac{PLT} were equal or smaller than those obtained through hard evaluation, as seen in Table~\ref{table:res-SS-prop-eval-train-noise}. This is particularly evident when~$n_{train} = 0$. The number of leaves remains unchanged as these methods do not affect training.

On the contrary, Table~\ref{table:res-SS-prop-eval-test-noise} shows that \ac{SE} accuracy was superior to that of hard evaluation in Experiment~2 for $n_{test} \ge 0.05$. However, this increase was not statistically significant. It suggests that the uncertainty distribution considered in the \ac{SE} approach better captures the noise added to the data, compared to~\ac{PLT}.

\section{Discussion}

We propose a probabilistic \ac{DT} approach to handle uncertain data, which separates the uncertainty model in three independent algorithm components. Our experiments evaluate these components in their ability to handle varying degrees of noise in the training and test data.

The first observation is that corrupting the data decreases the accuracy of the predictions, specially if the noise is in the test data. Accordingly, learning on data with increasing uncertainty results in \acp{DT} with a larger number of leaves, as the models attempt to learn the particularities of the training set.

The results indicate that \ac{SS}, \ac{STP} and \ac{UDT} are at least as robust to noise as C4.5, with non-significant improvements in accuracy. This was observed both for training or test data noise. \Ac{UDT} results are consistent with previous experiments with accuracy gains in the order of $3\%$, but where \ac{DT} size was not reported~\cite{Tsang2011}.

All soft training methods had longer running times compared to C4.5, the slowest being~\ac{UDT}. When comparing the search approaches, we observe that, by employing the discretization in Equation~\ref{eqn:discretization}, \ac{SS} increases the set of possible thresholds compared to C4.5, while preventing the computation of the information gain for values $x^{(i)}$ that are very close. \Ac{UDT} generates $s$ samples for each measurement. The number of entropy computations per attribute is bounded by~$s |\mathcal{D}|$, and therefore grows with the size of the dataset. Using \ac{SS}, this number is bounded by $\tfrac{1}{\delta} ( \tau_{max} - \tau_{min} )$, and does not grow with $|\mathcal{D}|$.

While maintaining accuracy, \ac{SS} and \ac{STP} led to significantly smaller~\acp{DT}. All approaches built larger trees for increasing $n_{train}$, as they start to overfit to the noise. \Ac{SS} and \ac{STP} were able to cope with this by reducing the number of splits, acting as regularizers. This can be interpreted as a consequence of having class probability estimates that reflect the uncertainty, illustrated in Figure~\ref{fig:perror-ex1a} and Table~\ref{tab:tab1}, on the pruning algorithm. The C4.5 pruning method uses a statistical test on the training data to make a pessimistic estimate the generalization error. Using a soft training approach expresses less confidence in the data, causing the pruning algorithm to remove more nodes. Tree size reductions were also observed for the multivariate sigmoid-split approach~\cite{irsoy2012soft}. However, they were most likely caused by the use of multivariate split functions, which are able to express complex rules more compactly, at the cost of reduced interpretability.

To investigate if the \ac{DT} size reduction could be obtained by changing the C4.5 pruning confidence factor $c$, in Appendix~C we show the result of varying $c$ when using C4.5, \ac{SS} or \ac{STP} for the first $5$ datasets of Table~\ref{tab:datasets}. For lower $c$, \ac{SS} resulted in smaller models with similar accuracy. In the overfitting range, this tendency is inverted. This indicates that the estimated class probabilities are more accurate with \ac{SS}, when the model actually learned representative splits. \Acp{STP} has led to consistently smaller trees than C4.5, except when $c$ is vey close to~$1$.

Given that the uncertainty models in \ac{SS}, \ac{STP} and \ac{UDT} are all Gaussian, an explanation for the disparity in results obtained by the proposed soft training approaches and \ac{UDT} could be the mismatch between the uncertainty model considered in the algorithms and the distribution of the noise added to the data. In our experiments, the same distribution was used in \ac{SS} and \ac{STP} and in the noise model, i.e. the standard deviation was a factor of the variable mean. In~\ac{UDT}, the standard deviation is instead a factor of the range of the data. This is also suggested by the values of the $w$ parameter selected through \ac{CV} displayed in Figures \ref{fig:params-train-noise} and~\ref{fig:params-test-noise}, which are smaller than $u_s$ and~$u_t$. Some degree of correlation between the parameters $u_s$, $u_t$ and the noise factors $n_{train}$, $n_{test}$ was observed. But no significant correlation was observed between $w$ and $n_{train}$ or $n_{test}$.

We hypothesize that the soft learning works more effectively as the uncertainty model approximates the real noise distribution. We therefore recommend the use of uncertainty models specified by domain experts. In our experiments, noise was simulated as an experimental proof-of-concept. To validate our approach on concrete clinical decision problems, the uncertainty distribution and its parameters shall be estimated \textit{a priori} for each variable. Such an estimation may be based, for instance, on the meta-analysis of clinical studies and on empirical clinical knowledge.

\begin{table}[ht]
	\centering
	\setlength{\tabcolsep}{4pt}

	\caption[]{{ \footnotesize Results for Experiment~1. Standardized metrics and absolute difference to C4.5 are presented, averaged over all datasets. The standardized results are tested for statistically significant differences according to a Wilcoxon signed-ranks test~\cite{Wilcoxon1945}, where each result is tested against C4.5 with the same~$n_{train}$. Standardized results with $p<.001$ highlighted in bold.}}
	\scriptsize
	\begin{tabular}{c|l|rrrrrrr}
		\toprule
		\multirow{2}{*}{Metric} 	& \multirow{2}{*}{Method}	& \multicolumn{7}{c}{Training data noise factor ($n_{train}$)} \\
		& 							& 0.00 & 0.05 & 0.10 & 0.20 & 0.30 & 0.40 & 0.50\\
		\midrule \midrule
		\multirow{6}{*}{\specialcell{No. leaves \\ standardized \\ by  C4.5 with \\ $n_{train}=0$}}
		& C4.5		& 0.00	& 0.44	& 0.72	& 1.09	& 1.42	& 1.75	& 2.03	\\
		& \Ac{SS}		& -0.48			 &	\textbf{-0.78} &	\textbf{-0.81} &	\textbf{-0.41} &	\textbf{-0.01} &	\textbf{0.40} &	\textbf{0.75} \\
		& \Ac{STP}		& \textbf{-1.11} &	\textbf{-2.19} &	\textbf{-2.06} &	\textbf{-1.95} &	\textbf{-0.90} &	\textbf{-1.01}	&	\textbf{-0.16} \\
		& \Ac{SE}		& 0.00	& 0.44	& 0.72	& 1.09	& 1.42	& 1.75	& 2.03	\\
		& \Ac{PLT}		& 0.00	& 0.44	& 0.72	& 1.09	& 1.42	& 1.75	& 2.03	\\
		& \acs{UDT}		& 0.46	& 0.64	& 0.83	& 1.25	& 1.81	& 1.97	& 2.17	\\
		\midrule

		\multirow{5}{*}{\specialcell{Absolute diff. \\ in no. leaves \\ to C4. with \\ same $n_{train}$ }}
		& \Ac{SS}		& -2.9	& -5.1 	& -5.6	& -5.4	& -5.4	& -5.3	& -5.0	\\
		& \Ac{STP}		& -5.4 	& -8.4 	& -8.6 	& -9.3	& -8.9	& -9.3	& -8.7	\\
		& \Ac{SE}		& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	\\
		& \Ac{PLT}		& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	\\
		& \acs{UDT}		& 1.5	& 0.43	& 0.6	& 0.7	& 1.1	& 0.4	& 0.2	\\
		\midrule
		\midrule

		\multirow{6}{*}{\specialcell{Accuracy \\ standardized \\ by C4.5 with \\ $n_{train}=0$}}
		& C4.5		& 0.00		& -0.52		& -0.87		& -1.27		& -1.59		& -1.65		& -2.06		\\
		& \Ac{SS}	& 0.26		& -0.14		& -0.39		& -0.79		& -1.05		& -1.23		& -1.40 	\\
		& \Ac{STP}	& 0.32		& -0.03		& -0.30		& -0.70		& -0.86		& -1.06		& -1.24		\\
		& \Ac{SE}	& -0.60		& -0.57		& -0.87		& -1.27		& -1.59		& -1.81		& -1.94		\\
		& \Ac{PLT}	& -0.51		& -0.67		& -1.01		& -1.38		& -1.68		& -1.75		& -2.14 	\\
		& \acs{UDT}	& 0.09		& -0.45		& -0.72		& -1.10		& -1.36		& -1.70		& -1.81 	\\
		\midrule

		\multirow{5}{*}{\specialcell{Absolute diff. \\ in accuracy (\%) \\ to C4. with \\ same $n_{train}$ }}
		& \Ac{SS}	& 0.8		& 1.2		& 1.7		& 1.6		& 1.8	&	1.3		&	2.1		\\
		& \Ac{STP}	& 0.9		& 1.6		& 2.0		& 1.9		& 2.3	&	1.9		&	2.6		\\
		& \Ac{SE}	& -2.0		& -0.3		& 0.0		& -0.1		& 0.0	&	-0.5	&	0.3		\\
		& \Ac{PLT}	& -1.2		& -0.5		& -0.4		& -0.3		& -0.2	&	-0.2	&	-0.2	\\
		& \acs{UDT}	& 0.3		& 0.1		& 0.5		& 0.4		& 0.6	&	-0.3	&	0.7		\\
		\midrule \midrule

		\multirow{5}{*}{\specialcell{Absolute diff. \\ running time (s) \\ to C4. with \\ same $n_{train}$ }}
		& \Ac{SS}	& 2.7	& 5.6	& 4.5	& 4.9	& 5.1	& 4.3	& 5.2	\\
		& \Ac{STP}	& 1.0	& 1.2	& 1.4	& 1.9	& 1.9	& 1.5	& 1.5	\\
		& \Ac{SE}	& 0.0	& 0.1	& -0.1	& 0.1	& 0.0	& 0.1	& 0.2	\\
		& \Ac{PLT}	& 0.1	& 0.2	& 0.2	& 0.2	& 0.2	& 0.2	& 0.2	\\
		& \acs{UDT}	& 40.6	& 45.0	& 37.6	& 49.0	& 53.9	& 37.3	& 44.2	\\

		\bottomrule
	\end{tabular} \label{table:res-SS-prop-eval-train-noise}
\end{table}
\begin{table}[ht]
	\centering
	\setlength{\tabcolsep}{4.2pt}

	\caption[]{{ \footnotesize Results for Experiment~2. Standardized metrics and absolute difference to C4.5 are presented, averaged over all datasets. The standardized results are tested for statistically significant differences according to a Wilcoxon signed-ranks test~\cite{Wilcoxon1945}, where each result is tested against C4.5 with the same~$n_{test}$. Standardized results with $p<.001$ highlighted in bold.}}
	\scriptsize
	\begin{tabular}{l|l|rrrrrrr}
		\toprule
		\multirow{2}{*}{Metric} 	& \multirow{2}{*}{Method}	& \multicolumn{7}{c}{Test data noise factor ($n_{test}$)} \\
		& 							& 0.00 & 0.05 & 0.10 & 0.20 & 0.30 & 0.40 & 0.50\\
		\midrule \midrule
		\multirow{6}{*}{\specialcell{No. leaves \\ standardized \\ by  C4.5 with \\ $n_{test}=0$}}
		& C4.5			& 0.00	& 0.00	& 0.00	& 0.00	& 0.00	& 0.00	& 0.00	\\
		&	\Ac{SS}	& -0.48				&	\textbf{-1.10}	& \textbf{-1.07}	& \textbf{-0.81}	& \textbf{-0.93}	& \textbf{-0.81}	& \textbf{-0.77}	\\
		&	\Ac{STP}	& \textbf{-1.11}	&	\textbf{-1.73}	& \textbf{-1.89}	& \textbf{-2.05}	& \textbf{-1.90}	& \textbf{-2.16}	& \textbf{-1.98}	\\
		&	\Ac{SE}		& 0.00				&	0.00			& 0.00				& 0.00				& 0.00				& 0.00				& 0.00				\\
		&	\Ac{PLT}	& 0.00				&	0.00			& 0.00				& 0.00				& 0.00				& 0.00				& 0.00				\\
		&	\acs{UDT}	& 0.46				&	0.62			& 0.49				& 0.45				& 0.47				& 0.33				& 0.49				\\
		\midrule

		\multirow{5}{*}{\specialcell{Absolute diff. \\ in no. leaves \\ to C4. with \\ same $n_{test}$ }}
		& \Ac{SS}		& -2.9	& -4.3	& -4.3	& -3.4	& -4.0	& -3.5	& -3.4	\\
		& \Ac{STP}		& -5.4	& -7.0	& -7.3	& -7.7	& -7.5	& -7.7	& -7.5	\\
		& \Ac{SE}		& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	\\
		& \Ac{PLT}		& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	\\
		& \acs{UDT}		& 1.5	& 1.7	& 1.3	& 1.2	& 0.9	& 0.9	& 1.5	\\
		\midrule \midrule

		\multirow{6}{*}{\specialcell{Accuracy \\ standardized \\ by C4.5 with \\ $n_{test}=0$}}
		& C4.5								& 0.00				& -0.97				& -1.38				& -2.01		& -2.59	& -3.06	& -3.50	\\
		& \Ac{SS}							& 0.26				& -0.50				& -0.82				& -1.45		& -1.85 & -2.27	& -2.62	\\
		& \Ac{STP}							& 0.32				& -0.40				& -0.72				& -1.18		& -1.54	& -1.91	& -2.15	\\
		& \Ac{SE}							& -0.60				& -0.85				& -1.19				& -1.79		& -2.32	& -2.78	& -3.19	\\
		& \Ac{PLT}							& -0.51				& -1.16				& -1.60				& -2.28		& -2.90	& -3.43	& -3.86 \\
		& \acs{UDT}							& 0.09				& -0.64				& -1.02				& -1.74		& -2.39	& -2.94	& -3.46 \\
		\midrule

		\multirow{5}{*}{\specialcell{Absolute diff. \\ in accuracy (\%) \\ to C4. with \\ same $n_{test}$ }}
		& \Ac{SS}	& 0.8	&	1.5		&	2.0		&	1.8		&	2.2		&	2.4		&	2.7		\\
		& \Ac{STP}	& 0.9	&	1.7		&	2.1		&	2.5		&	3.2		&	3.4		&	4.1		\\
		& \Ac{SE}	& -2.0	&	-0.2	&	0.0	&	-0.1	&	-0.1		&	-0.2		&	-0.1		\\
		& \Ac{PLT}	& -1.2	&	-0.5	&	-0.5	&	-0.6	&	-0.6	&	-0.7	&	-0.7	\\
		& \acs{UDT}	& 0.3	&	0.6		&	0.8		&	0.5		&	0.3		&	0.1		&	-0.1	\\

		\midrule \midrule

		\multirow{5}{*}{\specialcell{Absolute diff. \\ running time (s) \\ to C4. with \\ same $n_{train}$ }}
		& \Ac{SS}	& 2.05	& 4.8	& 6.0	& 4.8	& 4.7	& 3.3	& 5.29	\\
		& \Ac{STP}	& 1.7	& 2.8	& 2.8	& 3.2	& 3.7	& 3.4	& 3.40	\\
		& \Ac{SE}	& 0.1	& 0.1	& 0.1	& 0.1	& 0.0	& 0.0	& 0.0	\\
		& \Ac{PLT}	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	& 0.0	\\
		& \acs{UDT}	& 40.9	& 37.7	& 33.9	& 38.1	& 39.0	& 40.3	& 39.9	\\

		\bottomrule

	\end{tabular} \label{table:res-SS-prop-eval-test-noise}
\end{table}

\clearpage

\begin{figure*}
	\captionsetup[subfigure]{justification=centering} \centering
	\begin{subfigure}{\textwidth} \centering
		\includegraphics[width=.8\textwidth]{figure-res-train-noise-leaves.eps} \vspace{-3mm}
		\caption{Standardized number of leaves of Experiment 1 - training data noise.\\~\\} \label{fig:res-train-noise-leaves}
	\end{subfigure}
	\begin{subfigure}{\textwidth} \centering
		\includegraphics[width=.8\textwidth]{figure-res-train-noise-acc.eps} \vspace{-3mm}
		\vspace{-1mm}
		\caption{Standardized test accuracy (\%) of of Experiment 1 - training data noise.\\~\\} \label{fig:res-train-noise-acc}
	\end{subfigure}

	\begin{subfigure}{\textwidth} \centering
		\includegraphics[width=.8\textwidth]{figure-res-test-noise-leaves.eps} \vspace{-3mm}

		\caption{Standardized number of leaves of Experiment 2 - test data noise.\\~\\} \label{fig:res-test-noise-leaves}
	\end{subfigure}
	\begin{subfigure}{\textwidth} \centering
		\includegraphics[width=.8\textwidth]{figure-res-test-noise-acc.eps} \vspace{-3mm}
		\caption{Standardized test accuracy (\%) of Experiment 2 - test data noise.\\~\\} \label{fig:res-test-noise-acc}
	\end{subfigure}

	\caption{Results of the experiments displayed as boxplots of the standardized metrics for all datasets. (a,b)~-~Experiment 1: models trained on data with increasing levels of  noise~$\sim~\mathcal{N}(0, n_{train}\bar{x})$, and evaluated on data without added noise. (c,d)~-~Experiment 2: models trained on data without added noise, and evaluated on data with noise~$\sim~\mathcal{N}(0, n_{test}\bar{x})$. The baseline was obtained with~C4.5.}
	\label{fig:res-box-plots}
\end{figure*}

\clearpage

\Ac{SE} and \ac{PLT} did not improve accuracy given noisy training data, compared to the standard hard split approach. When noise was added to the test data, \ac{SE} led to non-significant increases in accuracy. This suggests that modeling uncertainty to target training data noise is only effective when this model is incorporated in the training phase, and not during evaluation. As such, we do not recommend the use of soft evaluation to handle training noise.

The disparity between the \ac{SE} and the \ac{PLT} results may also be explained by the consistency between the uncertainty model considered by the algorithms, and the model used to corrupt the data. \Ac{PLT} has demonstrated 2-3\% error rate reductions on a previous experiment using a single dataset, where the shape of the uncertainty model had been optimized.

\section{Conclusions}

This paper presents a probabilistic \ac{DT} learning approach to handle the uncertainty in the data, where the uncertainty model is separated in three independent algorithm components. The context is providing interpretable models for clinical decision support, with the motivation that the acknowledgement of uncertainty will facilitate the adoption of automated learning approaches in practice.

Previous \ac{DT} algorithms suggest the potential of probabilistic approaches to improve prediction robustness, and the need for an evaluation on more datasets and levels of noise. The impact of considering an uncertainty model in the learning phase or during evaluation was not however reported, as well as the impact of having noise is the training examples or the test examples.

In our approach, the uncertainty representation is incorporated: in the learning phase when searching for the optimal thresholds (\ac{SS}), when propagating the training data through the tree (\ac{STP}), and in the evaluation phase when obtaining predictions for unseen data (\ac{SE}). Any model can be chosen to capture the uncertainty. Our purpose is to incorporate clinical knowledge about the reliability of measurements. But as a proof-of-concept, we model the uncertainty as normally-distributed noise.

In our experiments, corrupting the test data seems to have a more severe impact on accuracy than adding noise to the training data. Upon increased noise, the soft training components, \ac{SS}, \ac{STP} and \ac{UDT}, show maintained or improved accuracy compared to C4.5. \Ac{STP}~and \ac{SS} act as regularizers, showing significant reductions in tree size, with \ac{STP} outperforming the latter. This was not the case of \ac{UDT}, possibly given the disparity between the noise model in the data and the uncertainty model in the algorithm. The running times of \ac{SS} and \ac{STP} were lower than those of \ac{UDT}. None of the soft evaluation approaches shows significant benefit compared to hard evaluation. Overall, we recommend using \ac{SS} and \ac{STP} with an uncertainty model that approximates as much as possible the real noise in the data. Finally, our study shows the importance of the acknowledgement of data uncertainty when learning decision models. Ideally, when designing clinical studies, an assessment of the reliability of each measurement should be considered part of the database.

Future work directions include evaluating the approach with domain-specific noise distributions, and studying the conditions under which the soft training provides benefit, regarding the complexity of the data. For highly separable data, the benefit of any soft approach is expected to be limited.

\ifCLASSOPTIONcompsoc

  \section*{Acknowledgments}
 This work was supported by the European Union Horizon 2020 research and innovation programme under grant agreement 642676 (Cardiofunxion), by the Spanish Ministry of Economy and Competitiveness (grant TIN2014-52923-R; Maria de Maeztu Units of Excellence Programme - MDM-2015-0502), by the European Union FP7 for research, technological development and demonstration under grant agreement VP2HF (611823), and FEDER.\\
\else

  \section*{Acknowledgment}
\fi

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\IEEEtriggeratref{32}

\bibliographystyle{IEEEtran}