Keywords: Prediction step, classification procedure on mixed data

Continuous data

Radar emitter patterns are partly described by continuous features that can be partially observed and approximately measured due to a noisy electromagnetic environment and sensor deficiencies. This chapter focuses on the development of a model that handles outliers and missing values to enable classification and clustering of radar emitters. First, continuous features of a radar emitter pattern are presented in Section 2.1 before introducing an experimental protocol developed to acquire realistic data since real military data are often classified. Then, the proposed model is explained in Section 2.2 where latent variables are introduced to model outliers and missing values. Inference procedure is processed through a Variational Bayesian Approximation in Section 2.3. Finally, evaluation of the model is proposed through two experiments and performance of the method are detailed in Section 2.4.

List of Abbreviations

In this thesis, several abbreviations are used and they are summarized here in order to facilitate reading and examination.

The abbreviations related to operational context are :

• EW : Electronic Warfare

• ESM : Electronic Support Measures

• ELINT : Electronic Intelligence

• SNR : Signal-to-Noise Ratio

The abbreviations related to state-of-the-art algorithms are :

• DA : Discriminant Analysis The abbreviations related to radar signal pattern are :

• PDW : Pulse Description Word

• TOA : Time-Of-Arrival

• RF : Radio Frequency

Introduction

In Electronic Warfare (EW) [START_REF] Schleher | Introduction to Electronic Warfare[END_REF], radar signal identification is a crucial component of Electronic Support Measures (ESM) systems [START_REF] Rogers | ESM processor system for high pulse density radar environments[END_REF]. ESM functions enable surveillance of enemy forces such as movements of enemy planes and warning of imminent attack such as launches of rockets. By providing information about the presence of threats, classification of radar signal has a self protection role ensuring that countermeasures against enemies are well-chosen by ESM systems [Wil82]. Furthermore, Electronic Intelligence (ELINT) functions focus on the interception and the analysis of unknown radar signals to update and improve EW databases. Then clustering of radar signals can play a significant role by detecting unknown signal waveforms and supporting ESM functions. Through its classification and clustering aspects, identification of radar signals is a supreme asset for decision making in military tactical situations. Depending on the information available in databases, the identification process can be distinguished into Source Emission Identification, also known as Radar Emitter Classification (REC), which concerns the classification of types of emission sources and Specific Emitter Identification (SEI) which focuses on recognition of copies of electromagnetic emission sources which are of the same type [START_REF] Dudczyk | Radar emission sources identification based on hierarchical agglomerative clustering for large data sets[END_REF]. REC practically relies on statistical analysis of radar signal patterns from distinct emitters [SL02, HZWT09, PJR13, YWY + 13, LJLC16, ZWCZ16, Che17, Sun18] whereas SEI aims to extract distinctive features in the process of signal processing to identify even a single copy of an emission source [GBB + 03, KO04, DK13, CLH14, Shi14, DK15b, Dud16]. Generally, the main di culties for SEI result from the lack of a precise and detailed description of a source emission model in databases [DK15a]. Indeed, information required for SEI can be exhaustive and are not always available in databases provided by operational entities. Therefore, this work only focuses on REC to meet operational constraints.

A radar signal [START_REF] Mark | Fundamentals of radar signal processing[END_REF] is conceived as a pulse-to-pulse modulation pattern in order to perform a specific role such as surveillance, missile guidance or short range tracking. The ability of a radar to perform such a role relies on its capacity to measure range and velocity of its targets. As a pulse-to-pulse modulation pattern, a radar signal pattern is decomposed into a relevant arrangement of sequences of pulses where each pulse is defined by continuous features and each sequence is characterized by categorical features. The continuous features of a pulse mainly refer to its time of interception, its radio frequency, its duration and its amplitude whereas the categorical features of a sequence refer to modulations of the continuous features. Then, a radar signal pattern is chosen as the combination of continuous and categorical features that minimizes ambiguities related to range measurements and velocity measurements. Depending on its expected function and the military context, radar signal patterns can be either simple or extremely complex. As an example, Multi-Function Radars (MFR) emitters [START_REF] Mackay | Tracking and control in multi-function radar[END_REF]BWW15], widely used in surveillance and tracking, are able to adapt their emitted patterns to a specific tactical situation they are operating in. Indeed, their emitted waveforms are designed to fit with characteristics of intercepted targets. Therefore in presence of multiple targets, their patterns can become extremely complex and provide to ESM systems a real challenge in terms of identification. As for the military context, it deeply has an e↵ect on the choice of patterns for radar emitters. Indeed, the continuous and categorical features related to radar emitter patterns mostly remain unchanged in a peace context and radar emitters can be mainly identified through their continuous features such as their pulse frequencies and pulse durations which refer to their spectral signature. On the contrary in a war context, radar emitters likely change their spectral signature to avoid being identified by enemy ESM systems which listed their continuous features during the previous peace context. Nonetheless, some categorical and continuous features remain identical since they characterize the way radar emitters operate in the electromagnetic environment. As an example, temporal evolution of radar emitter amplitudes completely characterizes scanning behaviours of radar emitters regardless of the military context.

Most of the time, ESM systems receive mixtures of signals from di↵erent radar emitters in the electromagnetic environment. Before identifying radar emitters, ESM systems have to isolate each radar signal from the received mixtures of radar signals. To this end, deinterleaving methods [START_REF] Mardia | New techniques for the deinterleaving of repetitive sequences[END_REF][START_REF] Milojevic | Improved algorithm for the deinterleaving of radar pulses[END_REF][START_REF] Moore | Deinterleaving pulse trains using discretetime stochastic dynamic-linear models[END_REF] are deployed as source separation algorithms to transform the homogeneous signal into a set of heterogeneous signals. Nonetheless, deinterleaving techniques cannot always manage to group all the pulses that belong to a radar emitter which results in a partial observation of its pattern. Furthermore, EW sensor deficiency and low Signal-to-Noise Ratio (SNR) values in sensors can also cause measurement errors [KP16] that disable detection of modulations related to radar signal patterns. When measurements are known to be erroneous, considering them as missing measurements can also be a more reliable approach than using them or discarding them. These material constraints introduce outliers in continuous radar data and missing components in both continuous and categorical radar data. At last, military databases are filled by human beings and may also be imperfect by gathering outliers and missing data due to human errors.

In statistical words, a radar signal pattern is described by continuous and categorical data which can be partially missing and erroneous. Depending on the complexity of radar signal patterns, the classification and clustering procedure should take into consideration any type of data and model a dependence structure to handle outliers and missing data. Classification and clustering problems are closely connected with pattern recognition [Bis06] where many general algorithms [HW79, EKS + 96, Bre01] have been developed and used in various fields [START_REF] Sander | Density-based clustering in spatial databases: The algorithm GDBSCAN and its applications[END_REF][START_REF] Anil | Data clustering: 50 years beyond K-means[END_REF]. However, most algorithms cannot handle missing data and imputation methods [TCS + 01] are required to generate data to use them. Hence, the main objective of this work is to define a classification and clustering framework that handles both outliers and missing values. Here, an approach based on mixture models is preferred since mixture models provide a mathematically based, flexible and meaningful framework for the wide variety of classification and clustering requirements [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. More precisely, a scale mixture of Normal distributions [AM74] is updated to handle outliers and missing data issues for any types of data. Exact inference in that Bayesian approach is unfortunately intractable, therefore a Variational Bayesian (VB) inference [WMR + 96] is used to find approximate posterior distributions of parameters and to provide a lower bound on the model log evidence used as a criterion for selecting the number of clusters.

Outline of the thesis is as follows. In Chapter 1, classification and clustering methods from state of the art are first presented according to their degree of supervision. After detailing supervised and unsupervised learning methods dedicated to classification and clustering, the selected approach is introduced through its theoretical aspects by defining mixture models as a flexible probabilistic framework that can handle both classification and clustering applications. However, estimation of parameters can turn out to be a cumbersome task and an approximation method is proposed to overcome this issue. Once theoretical aspects of mixture models have been presented, mixture models for continuous data are studied in Chapter 2 where generalizations of standard Gaussian mixture models are developed to handle outliers and missing data issues. The resulting model is performed on realistic simulated data obtained through an experimental protocol reproducing faults from real acquisition systems. Chapter 3 is dedicated to the extension of the model for mixed data composed of continuous and categorical features. After presenting categorical features of radar emitters by defining di↵erent types of modulations, a dependence structure between mixed data is investigated in order to develop a mixture model that handles both continuous and categorical data even in presence of outliers and missing values. The resulting model is performed on simulated data issued from a real-world database gathering various radar emitter patterns. Then Chapter 4 focuses on integration of temporal evolution radar data into the mixture model framework to take into consideration temporal evolution of radar emitter amplitudes that significantly reveal radar emitter scanning behaviours. To this end, the temporal evolution of radar emitter amplitudes is assumed to be either parabolic or piecewise parabolic and resulting models are evaluated on real operational cases to exhibit their performance. At last, an overall conclusion and work perspectives are given to conclude this thesis.

Chapter 1

State of the art and selected approach

Classification is used mostly as a supervised learning method to achieve a predictive goal [START_REF] Veyssieres | Identification of vegetation state and transition domains in California's hardwood rangelands[END_REF] by extrinsically adding unlabelled groups of data to reference classes. As for clustering, it is used for unsupervised learning to achieve a descriptive goal by discovering new groups of interest in data via an intrinsic assessment [RM05]. The main goal of this work is to develop a framework that handles both classification and clustering for mixed-type data gathering outliers and missing values. Then in this chapter, general state-of-the-art algorithms are first detailed in Section 1.1 before introducing the selected approach based on mixture models in Section 1.2.

State of the art

According to classification and clustering methods, three families of methods can be distinguished. Partitioning methods focus on relocating observations by moving them from one cluster to another conditionally to an initial partitioning. These partitioned-based methods generally require an a priori number of clusters set by the user. Then, density-based methods assume that each cluster is distributed according to a specific probability distribution [START_REF] Je↵rey | Model-based Gaussian and non-Gaussian clustering[END_REF] involving that the resulting marginal distribution of the data follows a mixture of distributions related to clusters. These density-based methods are designed for discovering clusters of arbitrary shape and for identifying their distribution parameters. At last, model-based methods aim to optimize the fit between data and chosen mathematical models by finding characteristic descriptions for each class or cluster [RM05]. Classification and clustering methods from these three families are successively presented in subsections 1.1.1 and 1.1.2.

Supervised learning

General state-of-the-art classification algorithms are introduced in this subsection.

Discriminant Analysis

Discriminant Analysis (DA) is a generalization of Fisher's linear discriminant [Fis36], a method used to find a linear combination of features that characterizes or separates two or more classes of objects or events. The resulting combination may be used as a linear classifier, or, more commonly, for dimensionality reduction before later classification. DA works when the measurements made on independent variables for each observation are continuous quantities. DA is used on labelled data to learn model parameters and then DA can perform classification. Moreover, the analysis is quite sensitive to outliers and the size of the smallest group must be larger than the number of predictor variables. In the case where there are more than two classes, the analysis used in the derivation of the Fisher discriminant can be extended to find a subspace which appears to contain all of the class variability.

DA approaches the problem by assuming that the conditional probability density functions p(x|y = 0) and p(x|y = 1) are both normally distributed with mean and covariance parameters (µ 0 , 0 ) and (µ 1 , 1 ). Under this assumption, a classification rule is built by computing the log of likelihoods ratio and testing if it is higher than some threshold '. The obtained classification rule is given by

(x ≠ µ 0 ) T ≠1 0 (x ≠ µ 0 ) + ln 0 ≠ (x ≠ µ 1 ) T ≠1 1 (x ≠ µ 1 ) ≠ ln 1 > ' . (1.1)
The resulting classifier is referred to Quadratic Discriminant Analysis (QDA) since the classification rule (1.1) is quadratic according to data x. The classification rule (1.1) can be linearly relaxed by assuming homoscedasticity ( 0 = 1 ) and becomes

a • x + b > 0 (1.2) with a = (µ 1 ≠ µ 0 ) T ≠1 0 , b = 1 2 1 µ T 0 ≠1 0 µ 0 ≠ µ T 1 ≠1 0 µ 1 ≠ ' 2 .
In that case, the resulting classifier is referred to Linear Discriminant Analysis (LDA). Performance of LDA and QDA classifiers are illustrated on Figure 1.1.

Logistic regression

Logistic regression [START_REF] Cox | The regression analysis of binary sequences[END_REF] is a powerful statistical way of modeling a binomial outcome with one or more explanatory variables. It measures the relationship between the categorical dependent variable and one or more independent variables by estimating probabilities using a logistic function, which is the cumulative logistic distribution. The predictive variables can be from any type ranges from continuous to categorical. Logistic regression can be extended to multinomial logistic regression [J.88] which is used to predict the probabilities of the di↵erent possible outcomes of a categorically distributed dependent variable, given a set of independent variables.

Mathematically, the outcome variable y is modeled conditionally to the explanatory variable x as follows

y = I 1 if -• x + ' > 0 0 otherwise . (1.3)
where ' follows a standard logistic distribution. The classification rule (1.3) can be interpreted in probabilistic way such that

y ≥ Y _ _ ] _ _ [ p(y = 0|x) = 1 1 + e ≠-•x p(y = 1|x) = e ≠-•x 1 + e ≠-•x

Decision trees and random forests

Decision trees [START_REF] Breiman | Classification and regression trees[END_REF] are decision support tools that use tree-like graphs or models of decisions and their possible consequences, including chance-event outcomes, resource costs, and utility. Each internal node represents a test on an attribute, each branch represents the outcome of the test, and each leaf node represents a class label. The paths from root to leaf represent classification rules. Decision trees are simple to understand and interpret since they can be generated from experts' rules based on mixed data. However, they are not robust to outliers since a small change in the data can lead to a large change in the structure of the optimal decision tree. This can be remedied by replacing a single decision tree with a random forest of decision trees, but a random forest is not as easy to interpret as a single decision tree.

Random Forests [Bre01] are an ensemble learning method [Rok10] for classification and regression that operates by constructing a multitude of decision trees at training time and outputting the mode of the classes (classification) or mean prediction (regression) of the individual trees. Random forests use boosting and bagging techniques [Sch90, Fre95, FS97] in order to correct for decision trees' habit of overfitting to their training set and having high variances [FHT + 00]. Indeed as an aggregation of multiple decision trees randomly trained on di↵erent feature sets of the training dataset [START_REF] Tin | The random subspace method for constructing decision forests[END_REF], random forests reduce the correlation between trees by avoiding over-focusing on features that appear highly significant in the training set but reveal less relevant in the test set. The ability of random forests to learn smoother decision rules than classical decision trees is visible on Figure 1.3. As predictive tools, random forests can not provide a description of features' relationships in datasets leading to infer on missing data.

Support vector machines

Support vector machines (SVM) [CV95] are supervised learning models with associated learning algorithms that analyze data used for classification and regression analysis. Given a set of binary examples, the SVM method maps that set into a transformed feature space [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF] where data are linearly separable and generates an hyperplane to separate those points into two distinct groups. SVM are scalable algorithms and can perform classification on large and sparse datasets. However, SVM methods focus on maximizing the margin between the two groups, they can not provide information about clusters' structure used to infer on missing data. An example of a decision rule learned by a SVM algorithm is presented on Figure 1.3.

Neural networks

Neural Networks (NN) [START_REF] Mcclulloch | A logical calculus of the ideas immanent in neurons activity[END_REF] were popularized by [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF] with networks called perceptrons making predictions based on a linear predictor function combining a set of weights with feature vectors. Perceptrons' limitations had been stated by [START_REF] Minsky | Perceptrons[END_REF] before [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF] introduced internal representation that enables a non linear mapping of data and ensures a better representation of the problem by adding hidden units in neural networks architecture. Thanks to the back-propagation procedure which distributes pattern recognition errors throughout the network, [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF] generalized the learning rule for multilayer networks [HSW89] demonstrated that standard multilayer feed-forward networks are capable of approximating any measurable function f to any desired degree of accuracy. Indeed, feed-forward networks define a mapping y = g(x; ◊) and learn values of parameter ◊ in order to find the best function approximation. Their architectures are straightforward in a sense that information flows through the function being evaluated from x, through the intermediate computations used to define g, and finally to the output y.

Mathematically, a feed-forward neural network with K layers is a function from a subset X 0 oe R n to a subset X K oe R p recursively defined by :

X k = g k (W k X k≠1 + b k ), k oe {1, . . . , K} (1.5)
where W k and b k are the k th layer weights and bias. X 0 and X K are the input and output layers whereas (X k ) 2:(K≠1) are hidden layers. The mapping g results in the composition of the (g k ) 1:K called activation functions.

Therefore, NN show strong results in classification [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] since they can extract features and learn classification rules (Figure 1.3) for a given architecture. However, NN are not descriptive models and can not provide explainable and relevant information about structure of classes.

Unsupervised learning

General state-of-the-art clustering algorithms are introduced in this subsection.

k-means

k-means algorithm [START_REF] John | Algorithm AS 136: A k-means clustering algorithm[END_REF] partitions observations into K clusters in which each observation belongs to the cluster with the nearest mean defined as the cluster whose mean has the least squared Euclidean distance (Figure 1.4). The k-means algorithm is described in procedure 1.2. The k-means algorithm is easily scalable and can be applied to large datasets without extra computational costs. Nonetheless, a key limitation is its cluster model which is based on isotropic clusters that are separable so that the mean converges towards the cluster center. The clusters are expected to be of similar size, so that the assignment to the nearest cluster center is the correct assignment. In addition, this algorithm is sensitive to noisy data and outliers and is limited to numeric attributes since it uses euclidean distance as metric. At last, it requires a value for the number of clusters k which is not trivial when no prior knowledge is available The k-prototypes algorithm was presented by [Hua98] to handle categorical data by defining the k-modes algorithm which uses a simple matching dissimilarity measure to deal with categorical objects, replaces the means of clusters with modes, and uses a frequency-based method to update modes in the clustering process to minimise the clustering cost function. Then, the k-prototypes algorithm results from integrating the k-means and k-modes algorithms to enable clustering of mixed-type data.

As for outliers and noise handling, [KR87] proposed the k-medoids method which di↵ers from the k-means mainly in its representation of the di↵erent clusters. Each cluster is represented by the most centric object in the cluster, rather than by the implicit mean that may not belong to the cluster. Hence, the k-medoids method is more robust than the k-means algorithm in the presence of noise and outliers since a medoid is less influenced by outliers or other extreme values than a mean. At last, the kernel k-means method was introduced as an extension of the k-means method by mapping the input data points non-linearly into a higher-dimensional feature space via a kernel function [DGK04]. The kernel k-means method enables discovering clusters with no arbitrary shape by relaxing the assumption on isotropic clusters.

Both presented updates of the k-means algorithm are more complex in nature and have a larger time complexity than the standard k-means algorithm. Moreover both methods still require the user to specify the a priori number of clusters K. Update each cluster center by computing the mean of all points assigned to it :

'k oe {1, . . . , K}, µ k = q N n=1 I zn=k x n q N n=1 I zn=k
until convergence return labels z and cluster centroids µ 1 , µ 2 , . . . , µ K

Density-based spatial clustering of applications with noise

Density-based spatial clustering of applications with noise (DBSCAN) [EKS + 96] is a densitybased clustering algorithm which groups together points with many nearby neighbors and marks as outliers points that lie alone in low-density regions. The mean idea behind DBSCAN is to continue growing a cluster as long as the density in the neighborhood exceeds a given threshold ' under the constraint that the neighborhood has to contain at least a minimum number of data points minP ts. An abstract algorithm for the DBSCAN algorithm is proposed in procedure 1.3. Advantages of the DBSCAN algorithm lie in its ability to discover clusters of arbitrary shapes even for large spatial databases. Indeed, DBSCAN is one of the most common clustering algorithms and also most cited in scientific literature since it can better learn the underlying structure of data (Figure 1.5). In addition to its robustness to outliers, DBSCAN does not Find the connected components of core points on the neighbor graph, ignoring all non-core points

Assign each non-core point to a nearby cluster if the cluster is an ' neighbor, otherwise assign it to noise return the estimated partition of the dataset

Spectral clustering

The main idea behind spectral clustering techniques [START_REF] Von | A tutorial on spectral clustering[END_REF] lies in the use of the spectrum of a given similarity matrix of the data to perform dimensionality reduction before clustering in fewer dimensions (Figure 1.5). The similarity matrix is provided as a symmetric matrix whose each element represents a measure of the similarity between data points. The general approach to spectral clustering is to use a standard clustering method on relevant eigenvectors of a Laplacian matrix of the similarity matrix. For computational e ciency, these eigenvectors are often computed as the eigenvectors corresponding to the largest several eigenvalues of a function of the Laplacian. When the relevant eigenvectors are processed through a k-means algorithm, the spectral clustering can be reformulated as a weighted kernel k-means problem [DGK04] where the kernel function is assimilated to the dimensionality reduction step leading to creation of the relevant eigenvectors. Moreover, spectral clustering can also be related to DBSCAN clustering [START_REF] Hess | The SpectACl of nonconvex clustering: a spectral approach to density-based clustering[END_REF] since optimal spectral clusters can correspond to density-connected components obtained by an asymmetric neighbor graph with edges removed when source points are not dense. Limitations of the spectral clustering lie in its computational cost and the choice of the similarity when data do not have a trivial structure.

Selected approach : mixture models

Here, an approach based on mixture models is preferred since mixture models provide a mathematicalbased, flexible and meaningful framework for the wide variety of classification and clustering requirements [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. Indeed mixture models, as generative models, enable the creation of a latent space where each latent variable can model a constraint of the problem of interest. Moreover, mixture models incorporate every degree of supervision since they handle both unsupervised or supervised classification problems. Finally, number of classes can be selected using criteria built on the model likelihood.

This section can be considered as a general theoretical framework used as a building block in the following chapters.

Definition

Mixture modeling [START_REF] Michael | Hierarchical mixtures of experts and the EM algorithm[END_REF] is a natural framework for classification and clustering. It can be formalized as :

p(x j | , K) = ÿ koeK a k  k (x j |◊ k ) , (1.6)
where x j oe X is an observation variable on the observation space X , K = {1, . . . , K} is the set of clusters and = (a, ◊ 1 , . . . , ◊ K ), with a = [a 1 , . . . , a K ] Õ , stands for parameters. Each probability distribution  k stands for the k th component distribution with a weight a k where a k Ø 0 and

q k a k = 1.
Assuming a dataset x oe X J of i.i.d observations (x 1 , . . . , x J ), the log likelihood function is given by

log p(x| , K) = ÿ joeJ log ÿ koeK a k  k (x j |◊ k ) ,
(1.7)

where J = {1, . . . , J}.

According to the degree of supervision, three problems can be distinguished : supervised classification, semi-supervised classification and unsupervised classification known as clustering. Supervised classification consists in parameters estimation of K known classes through a set of training data. Semi-supervised requires estimation of parameters of K unknown clusters whereas clustering proceeds to estimation of both parameters and number of clusters K.

Hence, parameters estimation is required to proceed to these three techniques. Unfortunately, classical Maximum Likelihood estimation turns out to be a complex problem since maximizing the log likelihood function (1.7) requires to deal with the summation over k that appears inside the logarithm and leads to a non closed form solution [Bis06].

One way of solving that estimation problem is to consider it as an incomplete data problem where only x is observed and where the complete data are composed of x and latent variables h such as labels of observations. The likelihood function for the complete dataset simply takes the form log p(x, h| ) and maximization of this complete-data log likelihood function should be straightforward. However since only the incomplete data x are given in practice, the complete data likelihood cannot be used and its expected value under the posterior distribution of the latent variable p(h|x, ) is considered. An elegant and powerful method for solving that issue is called the expectation-maximization algorithm (EM) [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] and consists in performing an expectation (E) step, which creates a function for the expectation of the complete log-likelihood evaluated using the current estimate for the parameters, and a maximization (M) step, which computes parameters maximizing the expected log-likelihood found on the E step. These parameter-estimates are then used to determine the distribution of the latent variables in the next E step.

The EM algorithm iterative scheme can be formalized as follows :

• E step : Calculate the expected value of the log likelihood function, with respect to the conditional distribution of h given x under the current estimate of the parameters t

Q( | t ) = E h|x, t [log p(x, h| )] ,
• M step: Find the parameters that maximize this quantity

t+1 = arg max Q( | t ) .
Nonetheless for some models, it can be infeasible to evaluate the posterior distribution or to compute expectations with respect to this distribution since the dimensionality of the latent space is too high to work with directly or because the posterior distribution has a highly complex form for which expectations are not analytically tractable [Bis06]. Hence approximation schemes are needed and rely on stochastic or deterministic approximations. Stochastic techniques such as Markov chain Monte Carlo have enabled the widespread use of Bayesian methods across many domains. They generally have the property that given infinite computational resource, they can generate exact results, and the approximation arises from the use of a finite amount of processor time. In practice, sampling methods can be computationally demanding, often limiting their use to small-scale problems. Also, it can be di cult to know whether a sampling scheme is generating independent samples from the required distribution. Deterministic approximation schemes, some of which scale well to large applications, are based on analytical approximations to the posterior distribution by assuming that it factorizes in a particular way or that it has a specific parametric form. As such, they can never generate exact results, and so their strengths and weaknesses are complementary to those of sampling methods [Bis06]. In this study, a deterministic approximation method known as Variational Bayes is developed for parameter estimation.

Latent variables

A mixture can be formalized as a latent model since the component label associated to each data point is unobserved. To this end, a categorical variable z j oe K can be considered to describe the index of the component distribution generating the observation variable x j . Then, the mixture distribution (1.6) is expressed as

p(x j | , K) = ÿ z j oeK p(x j |z j , , K)p(z j | , K) (1.8) where p(x j |z j , , K) = Ÿ koeK Â k (x j |◊ k ) " k z j , p(z j | , K) = Ÿ koeK a " k z j k .
and " k z j denotes the Kronecker symbol which is 1 if z j = k and 0 otherwise. The latent representation (1.8) can also be viewed in a hierarchical way as

x j |z j = k ≥ Â k (x j |◊ k )
(1.9)

z j ≥ Categorical(a) (1.10)
Then the joint distribution is

p(x j , z j | , K) = Ÿ koeK [a k  k (x j |◊ k )] " k z j .
Depending on the target problem, other latent variables can be introduced to model the data. h will refer to latent variables in the following parts.

Bayesian framework

Assuming a dataset x oe X J of i.i.d observations (x 1 , . . . , x J ) and independent latent data h = {h j } J j=1 , likelihood functions can be expressed as

p(x| , K) = Ÿ joeJ p(x j | ) , (1.11) p(x, h| , K) = Ÿ joeJ p(x j , h j | ) , (1.12)
where p(x, h| , K) is called the complete likelihood since it represents the joint distribution of the observed and latent data and J = {1, . . . , J}. That Bayesian framework imposes to specify a prior distribution for the parameters

p( |K) = p(a) Ÿ koeK p(◊ k ) .
Eventually, the posterior distribution of interest is obtained as

p(h, |x, K) = p(x, h| , K)p( |K) p(x|K) (1.13)
with the marginal distribution of data given by

p(x|K) = ⁄ p(x, h| , K)p( |K)ˆhˆ .

Inference

The Variational Bayesian inference was introduced in [WMR + 96] as an ensemble learning method for the mixtures of experts in order to avoid over-fitting and noise level under-estimation problems of traditional maximum likelihood inference. In [START_REF] Attias | Inferring parameters and structure of latent variable models by variational Bayes[END_REF], the Variational Bayesian inference was generalized for di↵erent types of mixture distributions and took the name Variational Bayes (VB). VB can be viewed as a Bayesian generalization of the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] combined with a Mean Field Approach [START_REF] Opper | Advanced mean field methods: Theory and practice[END_REF]. It consists in approximating the intractable posterior distribution p(h, |x, K) by a tractable one q(h, ) whose parameters are chosen via a variational principle to minimize the Kullback-Leibler (KL) divergence

KL [q||p] = ⁄ q(h, ) log 3 q(h, ) p(h, |x, K) 4 dhd .
Noting that p(h, |x, K) = p(x,h, |K) p(x|K) , the KL divergence can be written as

KL [q||p] = log p(x|K) ≠ L(q|K) .
L(q|K) is considered as a lower bound for the log evidence log p(x|K) and can be expressed as

L(q|K) = E h, [log p(x, h, |K)] ≠ E h, [log q(h, )] , (1.14)
where E h, [•] denotes the expectation with respect to q(h, ). Then, minimizing the KL divergence is equivalent to maximizing L(q|K). Assuming that q(h, ) can be factorized over the latent variables h and the parameters , a free-form maximization with respect to q(h) and q( ) leads to the following update rules :

VBE-step : q(h) Ã exp (E [log p(x, h| , K)]) , VBM-step : q( ) Ã exp (E h [log p( , x, h|K)]) .
The expectations E h [•] and E [•] are respectively taken with respect to the variational posteriors q(h) and q( ). Thereafter, the algorithm iteratively updates the variational posteriors by increasing the bound L(q|K).

Classification and clustering

According to the degree of supervision, three problems can be distinguished : supervised classification, semi-supervised classification and unsupervised classification known as clustering.

The supervised classification problem is decomposed into a training step and a prediction step. The training step consists in estimating parameters given the number of classes K and a set of training data x with known labels z. Then, the prediction step results in associating label z ú of a new sample x ú to its class k ú chosen as the Maximum A Posteriori (MAP) solution

k ú = arg max koeK p(z ú = k|x ú , , K)
given the previous estimated parameters .

In the semi-supervised classification, only the number of classes K is known and both labels z of the dataset x and parameters have to be determined. As for the prediction step, the MAP criterion is retained for a↵ecting observations to classes such that

k ú = arg max koeK p(z = k|x, , K) .
Given a set of data x, the clustering problem aims to determine the number of clusters K, labels z of data and parameters . Selecting the appropriate K seems like a model selection issue and is usually based on a maximized likelihood criterion given by

K = arg max K log p(x|K) .
( Then according to an a priori range of numbers of clusters {K min , . . . , K max }, the semisupervised classification is performed for each K oe {K min , . . . , K max } and both z K and K are estimated. Finally, the number of classes K in (1.15) is chosen as the maximizer of the lower bound L(q|K) :

K = arg max K L(q|K) .
(1.17)

After determining K, only z K and K are kept as estimated labels and parameters.

Conclusion

In this chapter, state-of-the-art classification and clustering algorithms have been presented. Some of them are dedicated to create boundaries to separate data into heterogeneous clusters such as LDA, SVM or NN whereas others focus on learning underlying structure of data to build them such as k-NN or k-means algorithms. However both types of algorithms do not provide an internal framework that infers on missing data and copes with any degree of supervision. Therefore, an approach based on mixture models is proposed and developed through its theoretical aspects. As hierarchical graphical models, mixture models provide a flexible framework to handle classification and clustering issues by introducing a latent space where each latent variable focuses on a specific constraint. However, the resulting model is not tractable and model learning is processed through an approximation method known as Variational Bayes Approximation. Eventually whatever degree of supervision is required, the number of classes K and parameters can be estimated to perform classification and clustering tasks. Next chapters deal with implementations of such models with di↵erent types of data.

Data

In this section, typical continuous radar features are first presented before introducing an acquisition system designed to generate realistic radar data which naturally embed outliers.

Continuous radar features

Continuous features of a radar emitter (Figure 2.1) are traditionally extracted from its Pulse Description Words (PDW). Each PDW gathers information related to a given pulse in the radar signal pattern such as

• its Time of Arrival (TOA) which is the time in µs at which the pulse is detected,

• its Amplitude (A) which the average measured amplitude of the pulse,

• its Radio Frequency (RF) which is the average measured frequency of the pulse in GHz,

• its Pulse Width (PW) which is the pulse duration in µs,

• its Pulse Repetition Interval (PRI) which is the di↵erence in µs between its TOA and the TOA of the previous pulse in the radar signal pattern.

The TOA of a pulse can be taken as the instant that a threshold is crossed. In presence of low signal-to-noise ratio (SNR), this measurement may be not precise and retaining the TOA of the first 3 dB is preferable [START_REF] Davies | Automatic processing for ESM[END_REF]. Moreover, the TOA is not an invariant feature since TOA sequence depends on the first observed TOA. Therefore, the PRI is retained as an invariant feature since it is the di↵erence between times of arrival of two successive pulses. Then, the RF of a pulse can be either fixed or modulated pulse-to-pulse. The RF is said to be frequency agile if it is randomly modulated pulse-to-pulse within fixed bounds and frequency hopping if it has systematic variations. Finally, the PW is the pulse duration chosen to ensure that a radar emits su cient energy such that reflected pulses are always detectable by its receiver. The amount of energy that can be delivered to a distant target is the product of two things; the output power of the transmitter, and the duration of the transmission. Therefore, pulse width constrains the maximum detection range of a target. Depending on sensor sensitivity, the PW may not be reliable and considering PW as missing data can be preferable.

Realistic data acquisition

In this subsection, an experimental protocol is introduced to acquire unclassified realistic data from di↵erent radar emitters. This protocol consists in an acquisition step followed by a feature extraction step.

The acquisition system is composed of two Software Defined Radio (SDR) platforms based on Ettus USRP E312 (Emitter) and B200 (Receiver) boards, linked to a laptop to record the data. As in [SEG + 16], this setup was chosen because it allows quick development and experimentation tasks on radiofrequencies from 70 MHz to 6 GHz, it is quite cheap and is available o↵-the-shelf. Radar waveforms, emitted by the URSP E312 board, are generated from bin files coded from a database gathering more than 40 typical radar waveforms with agile and hopping frequencies and jittered and staggered PRI. The developed system is presented in Figure 2.2. In order to meet hardware constraints, RF range was mapped to a 4MHz bandwidth and patterns of TOA and PW were slightly modified but their dynamic was preserved. Then, a threshold algorithm, provided by [START_REF] Davies | Automatic processing for ESM[END_REF], is used to detect pulses in the recorded signal s r (t) of duration T . Each pulse at TOA t is characterised by a triplet (PRI t , RF t , PW t ) t where PRI t is the di↵erence between TOA t and TOA t≠1 and the RF t feature is estimated with a Fast Fourier Transform (FFT) algorithm. Figure 2.3 shows parameters measurement on real data. For a given recorded signal s gathering n s pulses, the following PDW matrix is obtained Finally for each signal s j gathering n j pulses of the J recorded signals (s j ) J j=1 , a matrix P DW j is created from (2.1) and an observation vector x j is defined according to (2.2) such that

P DW = Q c c c c c c c
(d) (f) (b) (d) (f) (a) (c) (e) (a) (c) (e) (a)
x j = ( RF j , P W j , P RI j ) (2.2)
where RF j is the average value of RF, P W j is the average value of PW and P RI j is the average value of PRI defined in (2.3), (2.4) and (2.5).

RF j = 1 n j n j ÿ m=1 RF m , ( 2.3)
P W j = 1 n j n j ÿ m=1 P W m , ( 2.4)
P RI j = 1 n j n j ÿ m=1 P RI m . (2.5)
Once all observation vectors (x j ) J j=1 have been constructed, they are normalized to meet constraints of machine learning algorithms. Figure 2.5 shows the distribution of 150 normalized observation vectors of a radar emitter, where three outliers are visible.

Model

In this section, K emitters defined by continuous data are considered. Therefore, the main objective is to develop a mixture model which can build K distinct clusters even in presence of outliers and missing values. First, state-of-the-art approaches are reviewed. Then, the standard Gaussian mixture model is presented and two varieties of this model are introduced to handle outliers and missing values. At last, the proposed mixture model is developed into a Bayesian framework.

State of the art

Radar emitter classification relies on statistical analysis of Pulse Description Words (PDW) of a radar signal that gather its basic measurable parameters such as Radio Frequency (RF), Amplitude, Pulse Width (PW) or Pulse Repetition Interval (PRI). In terms of classification and clustering of emission sources from di↵erent types, many approaches based on data fusion and machine learning have been developed and traditionally proceed to feature extraction, dimensionality reduction and classification or clustering. For example, [SL02, [START_REF] Petrov | Radar emitter signals recognition and classification with feedforward networks[END_REF][START_REF] Li | Work mode identification of airborne phased array radar based on the combination of multi-level modeling and deep learning[END_REF][START_REF] Sun | Radar emitter classification based on unidimensional convolutional neural network[END_REF] propose various neural classification approaches based on the PDW structure of observed signals whereas [YWY + 13] introduce a hybrid radar emitter recognition method based on rough k-means and relevance vector machine and [Che17] develop an e cient classification method using weighted-xgboost model for complex radar signals in large datasets. As regards the clustering problem, [START_REF] He | Multi-parameter signal sorting algorithm based on dynamic distance clustering[END_REF] develop a dynamic clustering algorithm that uses designed distances and dynamic cluster centers and does not require fixing the number of classes which depends on the input data, [ZWCZ16] also introduce a clustering framework composed of local processing and multi-sensor fusion processing and use a Minimum Description Length criterion to update dynamically the number of clusters rather than setup in advance. These practical approaches mostly result from more general algorithms such as Random Forests [Bre01], Neural Networks [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF], Density-Based Spatial Clustering of Applications with Noise algorithm (DBSCAN) [EKS + 96] and k-means algorithm [START_REF] John | Algorithm AS 136: A k-means clustering algorithm[END_REF] which are also considered as state-of-the-art algorithms since they are used in various fields [START_REF] Sander | Density-based clustering in spatial databases: The algorithm GDBSCAN and its applications[END_REF][START_REF] Anil | Data clustering: 50 years beyond K-means[END_REF]. However, these practical and general algorithms can not handle missing data and imputation methods [TCS + 01] are required to generate data to use them. Hence, an approach based on mixture models is preferred since mixture models provide a mathematically based, flexible and meaningful framework for the wide variety of classification and clustering requirements [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. More precisely, a scale mixture of Normal distributions [AM74] is updated to handle outliers and missing data issues. On the one hand, this model is robust to outliers by accounting for the uncertainties of variances and covariances since the associated marginal distributions are heavy-tailed [START_REF] Archambeau | Robust Bayesian clustering[END_REF]. On the other hand, dependencies between features can easily be modelled through a multivariate Gaussian distribution in order to infer on missing values by benefiting from attractive Gaussian properties.

Standard Gaussian mixture models

Gaussian mixture models [QR78, JJ94] (GMM) are the most well-known mixture models for continuous data and have been widely used for decades. As a natural framework for classification and clustering, a GMM can be formalized as :

'j oe J , p(x j | , K) = ÿ koeK a k N (x j |µ k , k ) , ( 2.6)
where x j oe R d is an observation vector, K = {1, . . . , K} is a finite and known set of clusters and = (a, (µ k , k ) koeK ), with a = [a 1 , . . . , a K ] Õ , stands for parameters. Moreover, µ k and k are respectively the mean and the covariance matrix of the k th component distribution with a weight a k where a k Ø 0 and

q koeK a k = 1.
The GMM can be formalized as a latent model since the component label associated to each data point is unobserved. To this end, a categorical variable z j oe K can be considered to describe the index of the component distribution generating the observation variable x j . Then, the mixture distribution (1.6) is expressed as

p(x j | , K) = ÿ z j oeK p(x j |z j , , K)p(z j | , K) , ( 2.7)
where

p(x j |z j , , K) = Ÿ koeK N (x j |µ k , k ) " k z j , (2.8) p(z j | , K) = Cat(z j |a) = Ÿ koeK a " k z j k (2.9)
and " k z j denotes the Kronecker symbol which is 1 if z j = k and 0 otherwise.

Assuming a dataset x = (x j ) joeJ of i.i.d observations and independent labels z = (z j ) joeJ , the complete likelihood is obtained as follows

p(x, z| , K) = Ÿ joeJ Ÿ koeK [a k N (x j |µ k , k )] " k z j .
(2.10)

At last, the Bayesian framework imposes to specify priors for the parameters . The resulting conjugate priors are

Y _ _ _ _ _ _ ] _ _ _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) p(µ| , K) = Ÿ koeK N (µ k |µ 0 , ÷ ≠1 0 k ) p( |K) = Ÿ koeK IW( k |" 0 , 0 ) . (2.11)
where the Dirichlet and the Inverse Wishart distributions are defined as follows :

D(a|Ÿ) = c D (Ÿ) Ÿ koeK a Ÿ k ≠1 k , IW( |", S) = c IW (", S)| | ≠ "+d+1 2 exp 3 ≠ 1 2 tr(S ≠1 ) 4 ,
where c D (Ÿ) and c IW (", S) are normalizing constants such that

c D (Ÿ) = ( q koeK Ÿ k ) r koeK (Ÿ k ) , c IW (", S) = |S| " 2 2 d" 2 d ( " 2 )
.

The standard Gaussian model is shown on Figure 2.6.

CHAPTER 2. CONTINUOUS DATA K J " 0 0 µ 0 ÷ 0 k µ k x j z j a Ÿ 0 Figure 2
.6: Graphical representation of the standard Gaussian mixture model. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate represents the independent identically distributed observations x j and the indicator variables z j . Known quantities, respectively unknown quantities, are in blue, respectively in red.

Gaussian mixture models with missing data

As weighted sums of Gaussian distributions, GMMs benefit from attractive Gaussian properties that enable modeling dependencies between features to infer on missing data. Indeed, missing values can be handled by decomposing the features vector x j oe R d into observed features x obs Reminding that conditionally to its index cluster the features vector x j is Gaussian distributed as

x j = A x miss j x obs j B |z j = k ≥ N A µ k = A µ miss k µ obs k B , k = A miss k cov k cov Õ k obs k BB ,
the latent variable x miss j can be expressed as a Gaussian distributed variable such that

x miss j |x obs j , z j = k ≥ N 1 x miss j |µ x miss jk , x miss k 2
(2.12)

where

µ x miss jk = µ miss k + cov k obs≠1 k (x obs j ≠ µ obs k ) , (2.13
)

x miss k = miss k ≠ cov k obs≠1 k cov Õ k . (2.14)
Then, a marginal distribution for x obs

j is obtained such that x obs j ≥ N 1 x obs j |µ obs k , x obs k 2 with x obs k = 3 obs≠1 k + 2 ◊ obs≠1 k cov Õ k 1 x miss k 2 ≠1 cov k obs≠1 k 4 ≠1
.

Eventually, a Gaussian mixture model handling missing data is obtained by integrating missing data distributions (2.12) into the complete likelihood (2.10) such that

p(x, z| , K) = Ÿ joeJ Ÿ koeK [a k N (x j |µ k , k )] " k z j = Ÿ joeJ Ÿ koeK Ë a k N 1 x miss j |µ x miss jk , x miss k 2 N 1 x obs j |µ obs k , x obs k 2È " k z j . ( 2

Gaussian mixture models with outliers

A major limitation of GMMs is their lack of robustness to outliers that can lead to over-estimate the number of clusters since they use additional components to capture the tails of the distributions [START_REF] Svensén | Robust Bayesian mixture modelling[END_REF]. Nonetheless, outlier values in x j can be handled by introducing a latent variable u j to scale each mixture component covariance matrix k . That family of mixture models is known as scale mixtures of Normal distributions [AM74] and benefits from heavy-tailed marginal distributions accounting for the uncertainties of variances and covariances [START_REF] Archambeau | Robust Bayesian clustering[END_REF]. Introducing the latent positive variable u j into (2.8), the following scale component distribution is obtained

p(x j |u j , z j , , K) = Ÿ koeK N (x j |µ k , u ≠1 j k ) " k z j , (2.16)
and the joint distribution of (x j , u j ) is derived from (2.16) such that

p(x j , u j |z j , , K) = Ÿ koeK Ë N (x j |µ k , u ≠1 j k )p k (u j ) È " k z j
(2.17)

where p k (u j ) is the prior distribution of u j conditionally to z j = k.

Conditionally to the choice of a prior distribution for u j , the marginal distribution p(x j |z j , , K) = s OE 0 p(x j , u j |z j , , K)ˆu j of x j over u j can take di↵erent forms [START_REF] Martin | Scale mixtures of Gaussians and the statistics of natural images[END_REF]. [SB05, AV07, SZKL17, NW14] mainly propose using a Gamma distribution parametrized by a deterministic parameter ‹ k such that the joint distribution of (x j , u j ) from (2.17) becomes

p(x j , u j |z j , , K) = Ÿ koeK 5 N (x j |µ k , u ≠1 j k )G 3 u j | ‹ k 2 , ‹ k 2 46 " k z j (2.18)
Then the resulting marginal distribution p(x j |z j , , K) follows a Student-t distribution obtained

K J " 0 0 µ 0 ÷ 0 k µ k ‹ k u j x j z j a Ÿ 0 Figure 2
.8: Graphical representation of the Gaussian mixture model handling outliers. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations x j , the scale variables u j and the indicator variables z j . Known quantities, respectively unknown quantities, are in blue, respectively in red.

by

p(x j |z j , , K) = ⁄ +OE 0 p(x j , u j |z j , , K)ˆu j = ⁄ +OE 0 N 1 x j |µ k , u ≠1 j k 2 G 3 u j | ‹ k 2 , ‹ k 2 4 ˆuj = ( d+‹ k 2 ) ( ‹ k 2 )(‹ k fi) d 2 ◊ | ≠1 k | 1 2 ◊ 5 1 + 1 ‹ k (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) 6 ≠ d+‹ k 2 = T (x j |µ k , k , ‹ k ) (2.19)
where d is the dimension of the feature space and ‹ k is the degree of freedom of the Student-t distribution. Eventually, a Gaussian mixture model handling outliers is obtained by integrating (2.18) into the complete likelihood (2.10) such that

p(x, u, z| , K) = Ÿ joeJ Ÿ koeK 5 a k N 1 x j |µ k , u ≠1 j k 2 G 3 u j | ‹ k 2 , ‹ k 2 46 " k z j
.

where u = (u j ) joeJ are the scale latent variables related to continuous data x = (x j ) joeJ . As for prior distributions, ‹ k does not require a prior distribution since it is deterministic and other parameters are a priori distributed according to (2.11). Then, a graphical representation of the model is shown on Figure 2.8.

The degree of freedom variable ‹ k has been considered as a deterministic variable updated via an optimization argument during the maximization step of the VB inference [START_REF] Peel | Robust mixture modelling using the t distribution[END_REF] and [SB05, AV07, SZKL17, NW14] did not assume any prior distribution for ‹ k since there do not exist any known conjugate priors for ‹. For the sake of keeping conjugacy between prior and posterior distributions and adopting a full Bayesian treatment, a Gamma distribution G(u j |k ,k ) with shape and rate parameters

(-k , -k ) is chosen for p k (u j ) = p(u j |z = k) = G(u j |-k , -k ) such that the joint distribution of (x j , u j ) from (2.17) becomes p(x j , u j |z j , , K) = Ÿ koeK Ë N (x j |µ k , u ≠1 j k )G (u j |-k , -k ) È " k z j (2.20)
As in (2.19), the resulting marginal distribution p(x j |z j , , K) is also a Student-t distribution [DLGMD17] which is obtained as follows

p(x j |z j , , K) = ⁄ +OE 0 p(x j , u j |z j , , K)ˆu j = ⁄ +OE 0 N 1 x j |µ k , u ≠1 j k 2 G (u j |-k , -k ) du = (-k + d 2 ) (-k )(2-k fi) d 2 ◊ | ≠1 k | 1 2 ◊ 5 1 + 1 2-k (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) 6 ≠(-k + d 2 ) = T (x j |µ k , k , -k , -k )
Then, a Gaussian mixture model handling outliers is obtained by integrating (2.20) into the complete likelihood (2.10) such that

p(x, u, z| , K) = Ÿ joeJ Ÿ koeK Ë a k N 1 x j |µ k , u ≠1 j k 2 G (u j |-k , -k ) È " k z j . (2.21)
In a full Bayesian treatment, bothk andk require a prior distribution. Hence, a new conjugate prior distribution is introduced to avoid a non closed-form posterior distribution for (k ,k ). This prior distribution is defined below :

p(-k , -k |p 0 , q 0 , s 0 , r 0 ) Ã p -k ≠1 0 e ≠q 0 -k -s 0 -k k (-k ) r 0 I {-k >0} I {-k >0} (2.22)
where p 0 , q 0 , s 0 , r 0 > 0. The previous expression can be reformulated as :

p(-k , -k |p 0 , q 0 , s 0 , r 0 ) = p(-k |-k , s 0 , q 0 )p(-k |p 0 , q 0 , s 0 , r 0 ) with p(-k |-k , s 0 , q 0 ) = G(-k |s 0 -k + 1, q 0 ) , p(-k |p 0 , q 0 , s 0 , r 0 ) = 1 M 0 p -k ≠1 0 (s 0 -k + 1) q s 0 -k +1 0 (-k ) r 0 I {-k >0}
where

M 0 = ⁄ p -k ≠1 0 (s 0 -k + 1)
q s 0 -k +1 0 (-k ) r 0 I {-k >0} ˆ-k .
The normalization constant M 0 is intractable and a Laplace approximation method is derived to estimate it. As for other parameters, they are a priori distributed according to (2.11). At last, the resulting mixture model is shown on Figure 2.9.

Proposed mixture model

Varieties of the standard GMM have been introduced in (2.15) and (2.21) to handle outliers or missing values. The proposed model results from combining these two varieties to enable the handling of both outliers and missing values. Assuming a dataset x = (x j ) joeJ of i.i.d observations decomposed into observed features x obs = (x obs j ) joeJ and missing features x miss = (x miss j ) joeJ and independent latent variables u = (u j ) joeJ and z = (z j ) joeJ , the proposed model

K J " 0 0 µ 0 ÷ 0 k µ k -k -k u j x j z j a Ÿ 0 p 0 r 0 q 0 s 0 Figure 2
.9: Graphical representation of the Gaussian mixture model handling outliers with hyperparameters (-, -). The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations x j , the scale variables u j and the indicator variables z j . Known quantities, respectively unknown quantities, are in blue, respectively in red.

results in

p(x obs , h| , K) = Ÿ joeJ Ÿ koeK C a k N AA x miss j x obs j B - -µ k , u ≠1 j k B G(u j |-k , -k ) D " k z j = Ÿ joeJ Ÿ koeK Ë a k N 1 x miss j |µ x miss jk , u ≠1 j x miss k 2 N 1 x obs j |µ obs k , u ≠1 j x obs k 2 G(u j |-k , -k ) È " k z j (2.23)
where h = (x miss , u, z) is the set of latent variables and

= ! a = (a k ) koeK , µ = (µ k ) koeK , = ( k ) koeK , -= (-k ) koeK , -= (-k ) koeK "
is the set of parameters. Finally, the required prior distribution for is decomposed as follows

p( |K) = p(a)p(-, -|K)p(µ| , K)p( |K)
where prior distributions for a, µ and are given in (2.11) and the prior distribution of (-, -) is a product of distributions defined (2.22) such that

p(-, -|K) = Ÿ koeK p(-k , -k |p 0 , q 0 , s 0 , r 0 ) .
(2.24)

A graphical representation of the proposed model is shown on Figure 2.10.

K J " 0 0 µ 0 ÷ 0 k µ k -k -k u j x miss j x obs j z j a Ÿ 0 p 0 r 0 q 0 s 0 Figure 2
.10: Graphical representation of the proposed model. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations x j , the scale variables u j and the indicator variables z j . Known quantities, respectively unknown quantities, are in blue, respectively in red.

Inference

Direct inference on the proposed model is not trivial since the posterior distribution of latent missing data and parameters is intractable. Therefore, the Variational Bayes (VB) procedure is processed to estimate parameters of the mixture model defined in (2.23). Variational posterior distributions are obtained from the VB Expectation (VBE) and VB Maximization (VBM) steps. These variational posterior distributions are similarly obtained from classical posterior related in [SB05, AV07] and [MP04]. In addition to standard results, missing values are incorporated as latent variables in posterior calculations and a posterior distribution for missing data is proposed. At last, a lower bound on the log evidence is defined to master the convergence of the VB procedure.

Variational posterior distributions

Recalling that the VB procedure consists in approximating the intractable posterior distribution p(h, |x obs , K) by a tractable factorized distribution q(h, ) = q(h)q( ) that maximizes

L(q|K) = E h, Ë log p(x obs , h, |K) È ≠ E h, [log q(h, )] , ( 2.25)
where E h, [•] denotes the expectation with respect to q(h, ), variational posterior distributions are obtained by performing a free-form maximization through the following update rules :

VBE-step : q(h) Ã exp 1 E Ë log p(x obs , h| , K) È2 , VBM-step : q( ) Ã exp 1 E h Ë log p(x obs , h, |K) È2 .
According to a conditional factorization of q(h) and q( ) given by q(h) = q(x miss |u, z)q(u|z)q(z) , q( ) = q(a)q(µ, )q(-, -) , the following conjugate variational posterior distributions are obtained from the VB procedure

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ q(x miss |u, z) = Ÿ joeJ Ÿ koeK N 3 x miss j | μx miss jk , u ≠1 j ˜ x miss k 4 " k z j , q(u|z) = Ÿ joeJ Ÿ koeK G(-jk , -jk ) " k z j , q(z) = Ÿ joeJ Cat(z j |r j ) , q(a) = D(a| k) , q(-, -) = Ÿ koeK p(-k , -k |p k , qk , sk , rk ) , q(µ, ) = Ÿ koeK N (µ k | μk , ÷≠1 k k )IW( k |" k , ˜ k ) . ( 2

VBE-step

The VBE-step consists in deriving the following expectation

E Ë log p(x obs , h| , K) È = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] ≠ 1 2 3 d(log 2fi ≠ log u j ) + E [log | k |] + u j E Ë (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) È 4 + E [-k ] E [log -k ] + (E [-k ] ≠ 1) log u j ≠ E [log (-k )] ≠ u j E [-k ] 4 (2.27)
where '(j, k) oe J ◊ K :

E Ë (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) È = (x j ≠ μk ) T "k ˜ ≠1 k (x j ≠ μk ) + d ÷k (2.28)
is obtained from properties of the variational distribution q(µ, |K) in (2.26). Hence continuous data x are distributed a posteriori according to a product of normal distributions conditionally to latent variables u and labels z such that

x|u, z ≥ Ÿ joeJ Ÿ koeK N 1 x j | μk , u ≠1 j "≠1 k ˜ k 2 " k z j
Mean parameters ( μk ) koeK and variance parameters (" ≠1 k ˜ k ) koeK of these normal distributions are obtained from (2.28). By decomposing x into (x miss , x obs ) and by exploiting properties of the multivariate normal distribution, the following variational posterior distribution is obtained for missing values x miss :

q(x miss |u, x miss c , z, K) = Ÿ joeJ Ÿ koeK N 3 x miss j | μx miss jk , u ≠1 j ˜ x miss k 4 " k z j with '(j, k) oe J ◊ K : μx miss jk = μmiss k + ˜ cov k ˜ obs ≠1 k (x obs j ≠ μobs k ) , ˜ x miss k = ˜ miss k ≠ ˜ cov k ˜ obs ≠1 k ˜ cov Õ k "k .
Then by marginalising over x miss in (2.27), the expectation (2.27) becomes

⁄ E [log p(x, u, z| , K)] ˆxmiss = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + E [log | k |] ≠ 1 2 3 d obs j (log 2fi ≠ log u j ) ≠ log | ˜ x miss k | + u j 3 1 x obs j ≠ μobs k 2 T ˜ x obs j ≠1 k 1 x obs j ≠ μobs k 2 + d ÷k 4 4 + E [-k ] E [log -k ] + (E [-k ] ≠ 1) log u j ≠ E [log (-k )] ≠ u j E [-k ] 44 (2.29) with 'k oe K, ˜ x obs k = A ˜ obs ≠1 k + 2 ◊ ˜ obs ≠1 k ˜ cov Õ k 3 ˜ x miss j k 4 ≠1 ˜ cov k ˜ obs ≠1 k B ≠1 "k .
Conditionally to z, the scale latent variables u are distributed according to a product of Gamma distribution whose parameters are obtained by aggregating terms related to u such that

q(u|z, K) = Ÿ joeJ Ÿ koeK G 1 u j |-jk , -jk 2 " k z j with '(j, k) oe J ◊ K : -jk = E [-k ] + d obs j 2 , -jk = E [-k ] + 1 2 3 1 x obs j ≠ μobs k 2 T ˜ x obs j ≠1 k 1 x obs ≠ μobs k 2 + d ÷k 4 .
Finally, variational posterior categorical distributions are obtained for labels z by marginalising over u in (2.29) such that

⁄ E [log p(x, u, z| , K)] ˆxmiss ˆu = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + log (-jk ) ≠ -jk log -jk + E [-k ] E [log -k ] ≠ E [log (-k )] ≠ 1 2 3 d obs j log 2fi + E [log | k |] ≠ log | ˜ x miss k | 4 4 = ÿ joeJ ÿ koeK " k z j log fl jk (2.30)
where 'j oe J , k oe K,

log fl jk = E [log a k ] ≠ 1 2 3 d obs j log 2fi + E [log | k |] ≠ log | ˜ x miss k | 4 + E [-k ] E [log -k ] ≠ E [log (-k )] + log (-jk ) ≠ -jk log -jk .
(2.31)

Hence the variational categorical distributions are deduced from (2.30) and are given by

q(z|K) = Ÿ joeJ Cat(z j |r j )
where probabilities (r j ) joeJ are obtained from (2.31) such that 'j oe J , k oe K, rjk = fl jk q koeK fl jk .

VBM-step

The VBM-step consists in deriving the following expectation

E h Ë log p(x obs , h, |K) È = E h Ë log p(x obs , h| , K) È + p( |K) = ÿ joeJ ÿ koeK E h Ë " k z j È 3 log a k + ≠ 1 2 3 log | k | + d(log 2fi ≠ E h [log u j ]) + E h Ë u j (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) È 4 + -k log -k + (-k ≠ 1)E h [log u j ] ≠ log (-k ) ≠ E h [u j ] -k 4 + ÿ koeK (Ÿ 0 k ≠ 1) log a k + log c D (Ÿ 0 ) ≠ 1 2 3 (" 0 + d + 1) log | k | + Trace 1 0 ≠1 k 2 4 + c IW (" 0 , 0 ) + 1 2 3 d(log ÷ 0 k ≠ log 2fi) ≠ log | k | ≠ ÷ 0 k 1 µ k ≠ µ 0 k 2 T ≠1 k 1 µ k ≠ µ 0 k 2 4 ≠ log M 0 + (-k ≠ 1) log p 0 ≠ r 0 log (-k ) + s 0 -k log -k ≠ q 0 -k , (2.32)
where '(j, k) oe J ◊ K :

E h Ë u j (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) È = E h [u j ] (E h [x j ] ≠ µ k ) T ≠1 k (E h [x j ] ≠ µ k ) + Trace 1 V h [x j ] ≠1 k 2 (2.33)
is obtained from properties of the variational distribution q(h|K) with

E h [u j ] = -jk -jk , E h [x j ] = A μx miss j k x obs j B , V h [x j ] = Q a ˜ x miss k 0 d miss j ◊d obs j 0 d obs j ◊d miss j 0 d obs j ◊d obs j R b .
By factorizing terms related to a in (2.32), the following Dirichlet distribution is obtained

q(a|K) = D(a|Ÿ)
where

'k oe K, Ÿk = Ÿ 0 k + ÿ joeJ E h Ë " k z j È .
Then, (-, -) are a posteriori distributed according to the distribution defined in (2.22) such that

q(-, -|K) = Ÿ koeK p(-k , -k |p k , qk , sk , rk ) ,
where

pk = p 0 exp Q a ÿ joeJ E h Ë " k z j È E h [log u j ] R b , qk = q 0 + ÿ joeJ E h Ë " k z j È E h [u j ] , sk = s 0 + ÿ joeJ E h Ë " k z j È , rk = r 0 + ÿ joeJ E h Ë " k z j È .
By aggregating and factorizing terms related to each µ k in (2.32), a Normal distribution is obtained for each µ k such that

q(µ| , K) = Ÿ koeK N 1 µ k | μk , ÷≠1 k k 2 where 'k oe K, ÷k = ÿ joeJ E h Ë " k z j È E h [u j ] + ÷ 0 k , μk = q joeJ E h Ë " k z j È E h [u j ] E h [x j ] + ÷ 0 k µ 0 k ÷k .
Eventually, variance parameters are a posteriori distributed according to Inverse Wishart distributions given by

q( |K) = Ÿ koeK IW( k |" k , ˜ k )
where

"k = " 0 + ÿ joeJ E h Ë " k z j È , ˜ k = 0 + ÿ joeJ E h Ë " k z j È 1 E h [u j ] E h [x j ] E h [x j ] T + V h [x j ] 2 + ÷ 0 k µ 0 k µ T 0 k ≠ ÷k μk μT k .

Lower bound

The lower bound L(q|K) on the log evidence (2.25) is decomposed into the free energy E h, Ë log p(x obs , h, |K)

È and the entropy of the approximate posterior q(h, |K) given by E h, [log q(h, |K)]. The free energy can be developed as

E h, Ë log p(x obs , h, |K) È = E h, Ë log p(x obs , h| , K) È + E [log p( |K)]
where

E h, Ë log p(x obs , h| , K) È = ÿ joeJ ÿ koeK E h Ë " k z j È 3 E [log a k ] ≠ 1 2 3 d(log 2fi ≠ E h [log u j ]) + E [log | k |] + E h, Ë u j (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) È 4 + E [-k ] E [log -k ] + (E [-k ] ≠ 1)E h [log u j ] ≠ E [log (-k )] ≠ E h [u j ] E [-k ] 4 with E h, Ë u j (x j ≠ µ k ) T ≠1 k (x j ≠ µ k ) È = E h [u j ] 3 (E h [x j ] ≠ μk ) T "k ˜ ≠1 k (E h [x j ] ≠ μk ) d ÷k 4 + Trace 1 V h [x j ] "k ˜ ≠1 k 2 and E [log p( |K)] = ÿ koeK (Ÿ 0 k ≠ 1)E [log a k ] + log c D (Ÿ 0 ) ≠ 1 2 3 (" 0 + d + 1)E [log | k |] + Trace 1 0 E Ë ≠1 k È2 4 + c IW (" 0 , 0 ) + 1 2 3 d(log ÷ 0 k ≠ log 2fi) ≠ E [log | k |] ≠ ÷ 0 k E 5 1 µ k ≠ µ 0 k 2 T ≠1 k 1 µ k ≠ µ 0 k 2 6 4 ≠ log M 0 + (E [-k ] ≠ 1) log p 0 ≠ r 0 E [log (-k )] + s 0 E [-k ] E [log -k ] ≠ q 0 E [-k ] with E 5 1 µ k ≠ µ 0 k 2 T ≠1 k 1 µ k ≠ µ 0 k 2 6 = 1 μk ≠ µ 0 k 2 T "k ˜ ≠1 k 1 μk ≠ µ 0 k 2 + d ÷k .
As for the entropy term, the following decompositon is obtained

E h, [log q(h, |K)] = E h Ë log q(x miss q , u, z|K) È + E [log q( |K)] = E h Ë log q(x miss q |u, z, K) È + E h [log q(u|z, K)] + E h [log q(z|K)] + E [log q( |K)]
where

E h Ë log q(x miss q |u, z, K) È = ÿ joeJ ÿ koeK E h Ë " k z j È 1 2 3 d miss j (E h [log u j ] ≠ log 2fi ≠ 1) ≠ log | ˜ x j k | 4 , E h [log q(u|z, K)] = ÿ joeJ ÿ koeK E h Ë " k z j È 3 -jk log -jk ≠ log (-jk ) + (-jk ≠ 1)E h [log u j ] ≠ -jk E h [u j ] 4 , E h [log q(z|K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log rjk
and

E [log q( |K)] = ÿ koeK (Ÿ k ≠ 1)E [log a k ] + log c D (Ÿ k ) ≠ 1 2 3 (" k + d + 1)E [log | k |] + Trace 1 ˜ k E Ë ≠1 k È2 4 + c IW (" k , ˜ k ) + 1 2 3 d(log ÷k ≠ log 2fi) ≠ E [log | k |] ≠ ÷k E Ë (µ k ≠ μk ) T ≠1 k (µ k ≠ μk ) È 4 ≠ log M k + (E [-k ] ≠ 1) log pk ≠ rk E [log (-k )] + sk E [-k ] E [log -k ] ≠ qk E [-k ] with E Ë (µ k ≠ μk ) T ≠1 k (µ k ≠ μk ) È = d ÷k .

Expectations from variational distributions

Expectations developed in variational calculations are derived from properties of variational posterior distributions and are obtained as follows. Categorical distribution properties lead to 'j oe J , 'k oe K :

E h Ë " k z j È = rjk .
Dirichlet distribution properties lead to

'k oe K : E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b ,
where Â(•) is the digamma function. Gamma distribution properties lead to 'j oe J , 'k oe K :

E h [u j ] = -jk -jk , E h [log u j ] = Â(-jk ) ≠ log -jk .
Normal distribution properties lead to

'k oe K : E [µ k ] = μk , E Ë µ k µ T k È = V [µ k ] + E [µ k ] E [µ k ] T = ÷≠1 k k + μk μT k ,
Inverse Wishart distribution properties lead to

E Ë ≠1 k È = "k ˜ ≠1 k , E [log | k |] = log | ˜ k | ≠ d ÿ i=1 Â 3 "k + 1 ≠ i 2 4 ≠ d log 2 .
Posterior expectations ofk are derived from the posterior Gamma distribution (2.26) properties and can easily be computed conditionally to

-k E [-k ] = sk E [-k ] + 1 qk , E [log -k ] = E [Â (s k -k + 1)] ≠ log qk .
However, expectations depending onk are intractable

E [Â (s k -k + 1)] = ⁄ Â (s k -k + 1) p(-k |p k , rk )d-k , (2.34) E [-k ] = ⁄ -k p(-k |p k , rk )d-k , (2.35) E [log (-k )] = ⁄ log (-k )p(-k |p k , rk )d-k . (2.36)
Since lower bound calculation is required as a stop criterion, expectations (2.34), (2.35) and (2.36) have to be approximated. A deterministic method [START_REF] Tierney | Accurate approximations for posterior moments and marginal densities[END_REF] based on Laplace approximation is then applied. This method consists in approximating integrals of a smooth function times the posterior h(-)p(-|p, q, s, r) with an approximation proportional to a normal density in

◊ such that E[h(-] = h(-0 )p(-0 |p, q, s, r)(2fi) d-/2 | ≠ u ÕÕ (-0 )| 1/2 ,
where d -is the dimension of -, u(-) = log (h(-)p(-|p, q, r, s)) and -0 is the point at which u(-) is maximized.

In the case of unnormalized density q(-|p, q, r, s), Laplace's method can be applied separately to hq and q to evaluate the numerator and denominator here :

E[h(-] = s h(-)q(-|p, q, s, r)d- s q(-|p, q, s, r)d- .

Experiments

In this section, the proposed method is performed on the set of acquired data. For comparison, a standard neural network (NN), the k-nearest neighbours (k-NN) algorithm, Random Forests (RdF) the k-means algorithm and the DBSCAN are also evaluated. Two experiments are carried out to evaluate classification and clustering performance with respect to a range of percentages of missing values. First, characteristics for realistic data acquisition and imputation methods for missing data are detailed. Then, both experiments are described with their error measure and their performance are shown to exhibit the e↵ectiveness of the proposed model.

Data

Real data are acquired from the system detailed in Section 2.1. For each recording, the sampling frequency and the observation time T are respectively chosen as 4.17 MHz and 20 ms. The database exactly gathers 42 di↵erent radars waveforms and 150 observations are recorded for each waveform. Outliers and missing values are naturally embedded in observations due to material defects and real conditions detailed in Section 2.1. The dataset is shown in Figure 2.11. However, extra missing values are added to evaluate limits of the proposed approach. Missing information are introduced by randomly deleting coordinates of (x j ) 150 j=1 for each of the 42 radar emitters. Percentages of deletion range from 5% to 40%. Nevertheless, comparison algorithms do not handle datasets including missing values. Discarding observations that contain missing values can be a restrictive solution, therefore imputation methods have been developed [GLSGFV10]. In this chapter, two classical imputation methods, based on statistical analysis and machine learning, are performed. First, the mean imputation consists in filling a missing component of an observation by the average of observed values of that component. This method has the obvious disadvantage that it under represents the variability and also ignores correlations between observations [Sch97]. Then, imputation can be processed through a K-nearest neighbours method [HTS + 01] in order to replace missing values of an observation with a weighted mean of the k nearest completed observations where the weights are inversely proportional to the distances from the neighbours. Since replacements are influenced only by the most similar cases, the KNN method is more robust with respect to the amount and type of missing data [TCS + 01]. These imputation methods are compared with the proposed approach in terms of classification, clustering and reconstruction performance. For the comparison of reconstruction performance, mean-squared errors between original data and previous imputation methods are compared with Ÿ 0 ÷ 0 " 0 p 0 r 0 q 0 s 0 0.5 10 ≠4 4 0.9 1 1 1 the mean-squared error between original data and the variational posterior marginal mean of missing data given by

'j oe J , E q [x miss j ] = E q ( ⁄ ⁄ q(x miss j , u j , z j )du j dz j ) = ÿ koeK rjk μx miss jk .
(2.37)

Classification experiment

The classification experiment evaluates the ability of each algorithm to assign unlabeled data to one of the K classes trained by a set of labeled data. The classification task is decomposed into a training step and a prediction step defined in procedures 2.1 and 2.2. The training step consists in estimating variational parameters of q( ) defined in (2.26) given a set of training data with known labels. As for the prediction step, it results in associating new data to the class that maximizes their posterior probabilities. Since comparison algorithms do not handle datasets including missing values, a complete dataset is used to enable their training. During the prediction step, incomplete observations are either discarded and gathered in a reject class or completed thanks to the mean and KNN imputation methods. Standard configurations provided by Matlab are chosen for the RnF, the NN and the KNN algorithm. The proposed model and comparisons algorithms are trained on 70% of the initial database without extra missing values and tested on the remaining 30% of the database whose elements are randomly deleted according to di↵erent proportions of missing values. The RnF gathers 50 trees. The NN is composed of one hidden layer of 70 neurons and a softmax output layer and is trained with a cross-entropy loss. An accuracy metric is chosen for the classification experiment and observations belonging to the reject class are considered as misclassification errors. For each experiment, hyper-parameters are initialised as in Table 2.1 and 100 simulations are performed to take into account randomness of data deletion. Initialise Ÿ 0 , " 0 , ÷ 0 , µ 0 , 0 , p 0 , r 0 , s 0 and q 0 For the classification experiment, results are shown in Figure 2.12. Without missing data, both algorithms perfectly classify the 42 radar emitters. When the proportion of missing values increases, the proposed model outperforms comparisons algorithms and achieves an accuracy of 85% for 40% of deleted values whereas the accuracy of NN and KNN is lower than 50% with or without missing data imputation. As for the RnF, it outperforms both NN and KNN by achieving accuracies of 67% and 72% with standard imputation methods for 40% of deleted values. This higher performance of the proposed model reveals that the proposed method embeds a more e cient inference method than other imputation methods. That result is confirmed on Figure 2.12 where comparison algorithms are applied on data reconstructed by the proposed model. Indeed when the proposed inference is chosen, performance of NN and KNN increase up to 80% for 40% of deleted values and the RnF has almost the same performance than the proposed model. The Figure 2.12 also reveals that the proposed approach is more robust to missing data since it has a lower variance than other algorithms and imputation methods. Finally, this e ciency is shown on Figure 2.13 where the proposed model exhibits a lower mean-squared error for missing data imputation than the mean and KNN imputation methods. E↵ectiveness of the proposed model can be explained by the fact that missing data imputation methods can create outliers that deteriorate performance of classification algorithms whereas the inference on missing data and labels prediction are jointly estimated in the proposed model. Indeed, embedding the inference procedure into the model framework allows properties of the model, such as outliers handling, to counterbalance drawbacks of imputation methods such as outlier creation.

for iter = 1 to itermax do Update -jk , -jk , μx miss jk , ˜ x miss k Update Ÿk , ÷k , "k , pk , rr , sk , qk , μk , ˜ k Calculate the lower bound L if L iter ≠ L iter≠1 AE tol ◊ L iter≠1 then return ˜ train = 1 Ÿk , ÷k , "k , pk
Concerning the computational burden of the proposed approach, Figure 2.14 shows the evolution of computing times taken by model learning of the proposed model and comparison algorithms according to di↵erent numbers of observations. Considering that the learning of the proposed model is done o✏ine and that its code can be drastically optimized since it is only developed under Matlab, the computational burden of the proposed approach is acceptable. Indeed the proposed model is ten times slower than the RnF but shares similar computing times with the NN when the number of observations increases. Moreover, once the model learning has been performed o✏ine, predictions can be done online in real time.

Clustering experiment

The clustering experiment is composed of two experiments that aim to exhibit the clustering ability of each algorithm according to an a priori number of clusters K oe {K min , . . . , K max }. As developed in the previous chapter, the clustering algorithm is decomposed into two parts. First, a semi-supervised classification is performed for each K ranges from K min to K max to estimate variational parameters of q( , h) in (2.26) and labels of data in a mixture of K components. Then, the value of K that maximizes the lower bound (2.25) is retained as the posterior number of clusters as well as its associated parameters. According to the dataset visualised in Procedure 2.3 Semi-supervised classification procedure on continuous data Input: Unlabelled dataset x and number of classes K Output: Labels z and parameters ˜ Initialise Ÿ 0 , " 0 , ÷ 0 , µ 0 , 0 , p 0 , r 0 , s 0 and q 0 for iter = 1 to itermax do Update -jk , -jk , μx miss jk ,

˜ x miss k , rjk Update Ÿk , ÷k , "k , pk , rr , sk , qk , μk , ˜ k Calculate the lower bound L if L iter ≠ L iter≠1 AE tol ◊ L iter≠1 then return ˜ = 1 Ÿk , ÷k , "k , pk , rr , sk , qk , μk , ˜ k 2 koeK
and z such that each zj = arg max koeK rjk end if end for Procedure 2.4 Clustering procedure on continuous data Input: Unlabelled dataset x and a priori range of numbers of clusters K oe {K min , . . . , K max } Output: Labels z, parameters ˜ and optimal number of clusters K

for K = K min to K max do Perform semi-supervised classification with K classes Stock labels zK , parameters ˜ K and L K end for return z K and ˜ K such that K = arg max K L K Table 2.

2: Initialisation of hyper-parameters values for clustering on continuous data

Ÿ 0 ÷ 0 " 0 p 0 r 0 q 0 s 0 0.5 100 4 1 1 1 1 by using an heuristic proposed in the original paper [EKS + 96]. A supervised initialisation is retained for the proposed model due to its sensitivity to initialisation. It consists in initialising prior component means µ 0 from results of a k-means algorithm and prior component covariance matrices 0 from diagonal matrices whose diagonal elements are variances of observed features. Since comparison algorithms do not handle observations with missing values and do not provide a clustering result for them, missing data are either discarded and gathered in a reject class or completed thanks to the mean and KNN imputation methods before running these algorithms.

For each experiment, hyper-parameters are initialised as in Table 2.2 and 100 simulations are performed to take into account randomness of data deletion.

The first clustering experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K oe {K min , . . . , K max }. Performance are evaluated through the Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF] that compares estimated partitions of data with the ground-truth. Results of the first experiment on realistic data are shown in Figures 2.15 and 2.16. Without the presence of missing values, performance of DBSCAN and the proposed model are similar with an ARI of 97% (Figure 2.15) whereas the k-means algorithm ARI reaches 95% (Figure 2.16). When the proportion of missing values increases, the proposed model outperforms both DBSCAN and k-means and achieves an ARI of 87% for 40% of deleted values whereas the ARI of comparison algorithms is lower than 30% with standard missing data imputation. This higher performance reveals that the proposed method embeds a more e cient inference method than other imputation methods. That result is confirmed on both Figure 2.15 and Figure 2.16 where DBSCAN and k-means are applied on data reconstructed by the proposed model. Indeed, performance of both algorithms increase up to 77% and 69% for 40% of deleted values when the proposed inference is chosen.

The second experiment tests the ability of each algorithm to find the true number of clusters K among {K min , . . . , K max }. The lower bound (2.25) and the average Silhouette score [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF] are criteria used to select the optimal number of clusters for the proposed model and the k-means algorithm. Indeed, the ARI can not be used since it requires the ground-truth and DBSCAN automatically selects a number of clusters for a given dataset. Results of the second experiment on realistic data are visible on Figures 2.18 and 2.17. Figure 2.18 shows the evolution of the number of clusters estimated by DBSCAN according to di↵erent proportions of missing values and imputation methods. Since DBSCAN automatically estimates the number of clusters and manages outliers by creating new clusters, results on Figure 2.18 can be used to evaluate performance of imputations methods. For mean and k-NN imputation methods, DBSCAN estimates a number of clusters greater than 140 as proportion of missing values is equal or greater than 5%. When DSBCAN is run on the posterior reconstruction (2.37), the estimated number of clusters stays under 50 until 20% of missing values and reaches 70% for 40% of missing values. These performance indicate that the proposed approach creates less outliers than other imputation methods by providing a more robust inference on missing data since DBSCAN localizes less outliers in the posterior reconstruction (2.37) than in standard imputation methods. Figure 2.17 presents numbers of clusters selected by the lower bound and average Silhouette scores for the proposed model and k-means algorithm according to di↵erent proportions of missing values and imputation methods. Without missing data, the correct number of clusters (K=42) is selected by the two criteria for the k-means algorithm and the proposed model. In presence of missing values, the average Silhouette score always selects K = 72 when the k-means algorithm is run on data completed by standard imputation methods. When, the k-means algorithm performs clustering on the posterior reconstruction (2.37), the average Silhouette score correctly selects K = 42 until 15% of missing values and chooses K oe {62, 72} when the proportion of missing values is greater than 20%. Eventually when the proposed model does clustering, the two criteria select the correct number of clusters K = 42 for every proportion of missing values. These results show two main advantages of the proposed model. As previously, the proposed model provides a more robust inference on missing data since the average Silhouette score chooses more representative number of clusters when the k-means algorithm is run on the posterior reconstruction (2.37) than on data completed by standard imputation methods. Furthermore, since the lower bound criterion also selects the correct number of clusters as the average Silhouette score, it can be used as a valid criterion for selecting the optimal number of clusters and does not require extra computational costs as the Silhouette score since it is computed during the model parameters estimation. Finally, the proposed approach provides a more robust inference on missing data and a criterion for selecting the optimal number of clusters without extra computations.

Figure 2.19 shows the evolution of computing times taken by model learning of the proposed model and comparison algorithms according to di↵erent numbers of clusters and observations. As the model learning in Subsection 2.4.2, clustering is only performed o✏ine to extract information from radar signals recorded during operational missions. Even if the proposed method is ten times slower than the k-means algorithm, the computational burden of the proposed approach is still acceptable and meets operational requirements.

Conclusion

In this section, we propose a mixture model to classify and cluster radar emitters. Radar signals are often partially observed due to imperfect conditions of acquisition and deficient hardwares. Therefore to account for missing data and outliers, a scale mixture of Normal distributions, known for its robustness to outliers and its flexible framework for classification and clustering, is chosen. Moreover, thanks to the introduction of latent variables, the proposed model can infer on missing data. Since the posterior distribution is intractable, model learning is processed through a Variational Bayes inference where a variational posterior distribution is proposed for missing values. Experiments on various real data showed that the proposed approach handles both outliers and missing values and can outperform standard algorithms in classification and clustering tasks. Indeed the main advantage of our approach is that it allows properties of the model, such as outliers handling, to counterbalance drawbacks of imputation methods by embedding the inference procedure into the model framework.

Chapter 3

Mixed data

Continuous data describing radar emitter waveforms such as the Carrier Frequency, the Pulse Width and the Pulse Repetition Interval have been previously taken into account in order to cluster radar emitters. Nonetheless, these continuous features are frequently modulated to enhance functions of the radar emitters. Therefore, these modulations can be exploited as categorical features to cluster radar emitters. According to types of modulations, a dependence structure can be established to model conditional relations between continuous and categorical features. This dependence structure is then included into the previous mixture model to take advantage of specific patterns related to each radar emitter. This chapter contains four sections which focus on the integration of mixed data to enhance classification and clustering performance. Section 3.1 presents assumptions on continuous and categorical features of radar emitters. Section 3.2 introduces the dependence structure of mixed data and the proposed model. Section 3.3 details the inference procedure for the estimation of parameters related to the proposed model. At last in Section 3.4, various experiments are carried out to exhibit performance of the proposed approach.

Data

In this chapter, data consist of J pulses gathering J continuous features x q = (x qj ) joeJ and J categorical features x c = (x cj ) joeJ from K distinct emitters. Let x j = (x qj , x cj ) the j th observation vector of mixed variables where

• x qj oe R d is a vector of d continuous radar features such as the Radio Frequency, the Pulse Width and the Pulse Repetition Interval,

• x cj = 1 x 0 cj , . . . , x q≠1 cj 2
oe C q is a vector of q categorical radar modulations such as intrapulse modulations, pulse-to-pulse modulations or scanning types.

Radar features and distributions related to continuous and categorical data are presented in the following subsections.

Assumptions on continuous data

In this subsection, continuous radar features are first recalled from the previous chapter. Then, the distribution of continuous data is presented.

Radar Features

As in the previous chapter, continuous features of a radar emitter are extracted from its Pulse Description Words (PDW). Each PDW gathers the radio frequency (RF), the pulse width (PW) and the pulse repetition interval (PRI) of a given pulse in the radar signal pattern. These continuous features are exhibited on Figure 3.1.

Distribution of continuous features

The continuous features of the j th observation are modeled through x qj oe R d which is a vector of d continuous variable distributed according to a multivariate normal distribution with mean and variance parameters (µ, ). Indeed through its properties, the multivariate distribution can enhance the dependence structure between continuous variables in order to handle outliers and missing values.

Outliers Outliers for continuous data x q = (x qj ) joeJ are handled as the previous chapter by introducing the scale latent variables u = (u j ) joeJ such that

'j oe J , x qj |u j ≥ N (x qj |µ, u ≠1 j ) , u j ≥ G(u j |-, -) ,
where each u j follows a Gamma distribution with shape and rate parameters (-, -) oe R ú+ ◊ R ú+ .

Missing values Since continuous data x q = (x qj ) joeJ can be partially observed, they are decomposed into observed features x obs q = (x obs qj ) joeJ and missing features x miss q = (x miss qj ) joeJ such that = µ obs ,

'j oe J , x qj = A x miss qj x obs qj B with (x miss qj , x obs qj ) oe R d miss j ◊ R d obs
x miss q = miss ≠ cov obs ≠1 cov' ,

x obs q = 3 obs ≠1 + 2 ◊ obs ≠1 cov' 1 x miss q 2 ≠1 cov obs ≠1 4 ≠1
, and parameters (µ, ) are decomposed as follows

µ = A µ miss µ obs B , = A miss cov cov' obs B .

Assumptions on categorical data

In this subsection, categorical radar features are first introduced. Then, the distribution of categorical data is presented.

Radar Features

Categorical radar features are mainly related to continuous radar features since they describe modulations of a radar emitter pattern. Depending on the nature of continuous features, di↵erent types of categorical features can be taken into consideration. Considering the continuous features RF, PRI and PW of a radar emitter, categorical features are modulations related to these quantities. Then as regards a pattern of pulses, modulations of RF, PRI and PW are called pulse-to-pulse modulations when they are applied on a group of pulses whereas modulations of RF are called intrapulse modulations when they are applied on single pulses. These two families of modulations are used by radar emitters to achieve a common goal. Indeed, both families of modulations aim to obtain a higher resolution of targets by reducing ambiguities related to target range and target velocity. Furthermore, intrapulse modulations also focus on avoiding identification of radar emitters since they enable the minimization of these ambiguities even in presence of noise. Considering amplitudes of a radar emitter, the related categorical feature is the scanning behaviour of the radar emitter. These three types of categorical features are presented below.

Pulse-to-pulse modulations

= a ◊ i + v
where a and v are the slope and the intercept of the linear function.

Dwell and Switch modulation When parameter values (p i ) 1AEiAEn are piecewise constant, the emission of the n pulses is known to be dwelled or switched. Hence, for J disjoint subsets V j forming a partition of {1, . . . , n}, parameter values p i are dwelled or switched if they are constant on each V j such that

'i oe {1, . . . , n}, p i = J ÿ j=1 v j I ioeV j
where v j is the value of the piecewise V j .

Stagger modulation

The emission of the n pulses is staggered when parameter values (p i ) 1AEiAEn are distributed according a sequence of q moments (v 1 , . . . , v q ) such that

'i oe {1, . . . , n}, p i = v m(i)
where m : {1, . . . , n} ae {1, . . . , q} is a surjective application associating moments to parameters values. If m produces a periodic sequence of the q moments, the stagger is regular.

Wobble modulation

The emission of the n pulses is wobulated when parameter values (p i ) 1AEiAEn are repetitively modulated through a periodic pattern such that

'i oe {1, . . . , n}, p i = f (p i )
where f is a periodic pattern usually chosen as a sinus wave or a triangular wave. Jitter modulation The emission of n pulses is jittered when parameter values (p i ) 1AEiAEn are randomly generated around a nominal value v. Jittered emissions are commonly Gaussian such that 'i oe {1, . . . , n},

p i = v + '
where v is the nominal value and ' ≥ N (0, ‡ 2 ) is a Gaussian noise.

Intrapulse modulations To avoid identification of operating radars by ESM systems, radar designers have developed Low Probability of Intercept (LPI) waveforms. Theses waveforms are either frequency-modulated or phased-modulated in order to improve resolution for the radar emitter at the expense of a suboptimal signal-to-noise ratio (SNR) [START_REF] Levanon | Radar signals[END_REF]. In other words, theses pulse modulations enable the maximization of the target range and the range resolution of radars. On the contrary, ESM resolution is less accurate since LPI signals are embedded in much noise and the identification task can be compromised. Intrapulse modulations are presented below.

Frequency Modulation Signal By spreading energy over a modulation bandwidth, Frequency Modulation (FM) signals provide a better range resolution than constant frequency signals. A Linear FM (LFM) signal, also known as a chirp, is obtained by swepting linearly the frequency band during the pulse duration. However, the LFM can involve relatively high autocorrelation sidelobes. Therefore to counter that drawback, Nonlinear FM (NLFM) [Cos84] signals are used to provide a more accurate spectrum which is shaped by deviating the constant Phase Modulation Signal Phase coding is one of the first methods for pulse compression. The concept rests on dividing a pulse of duration T into M bits of identical duration t b = T M and assigning a di↵erent phase value to each bit. The main advantage of phase coding over frequency modulation is low peak side lobe level [START_REF] Levanon | Radar signals[END_REF]. Barker codes [START_REF] Barker | Group synchronization of binary digital systems[END_REF] are the most famous family of phase codes gathering 13 known binary sequences which were reported by [START_REF] Barker | Group synchronization of binary digital systems[END_REF] and [START_REF] Turyn | On Barker codes of even length[END_REF] and are given in Table 3.1. Other polyphase codes such as the Frank code [FZH62] and the Zado↵ code [Z + 63] are widely used for pulse compression. Scanning types A radar emitter can truly di↵er from another one through its scanning pattern. While searching for targets across the environment, its beam steering can behave di↵erently depending on its antenna shape, its composition and its mission. The most common scanning types are presented below and illustrated on Figure 3.5.

Circular Scan A circular scanning radar is a constant rotational scanning radar that provides accurate target range and azimuth information. It uses an antenna system that continuously scans 360°in azimuth making it ideal for the roles of early warning and initial target acquisition.

Sector Scan A sector scanning radar is a radar which scans unidirectionally or bidirectionally in a delimited sector. For example, the helical scan is a unidirectional scan pattern that enables a "pencil" beam to search a 360°pattern. As a bidirectional scan, the raster scan uses a thin beam to cover a rectangular area by scanning in azimuth and elevation.

Track-While-Scan (TWS) A track-while-scan (TWS) system generates two or more distinct radar beams that enable a radar to track multiple targets while scanning for others.

Electronic Scan An electronic scanning radar provides a computer-controlled scanning in which radar beams are electronically steered to point in di↵erent directions without moving the antenna.

Distribution of categorical features

The categorical features of the j th observation are modeled through

x cj = 1 x 0 cj , . . . , x q≠1 cj 2
oe C q which is a vector of q categorical variables where

C q = C 0 ◊ . . . ◊ C q≠1 is the tensor gathering each space C i = Ó m i 1 , . . . , m i |C i | Ô
of events that x i cj can take 'i oe {0, . . . , q ≠ 1}.

As continuous data x q , categorical data x c = (x cj ) joeJ can be partially observed. Hence x c are decomposed into observed features x obs c = (x obs cj ) joeJ and missing features

x miss c = (x miss cj ) joeJ such that 'j oe J , x cj = A x miss cj x obs cj B
with (x miss cj , x obs cj ) oe C q miss j ◊ C q obs j and q miss j + q obs j = q .

where C q miss j and C q obs j , are disjoint subsets of C q embedding missing features x miss cj and observed features x obs cj .

If each categorical feature x i cj of the j th observation x cj is assumed to be independent from other features and distributed according to a categorical distribution, a dependence structure between features cannot be modeled and inference on missing categorical features cannot be handled. Hence, a multivariate categorical distribution integrating a dependence structure for x cj is proposed to tackle that issue. As detailed below, x cj follows a multivariate categorical distribution MC(x cj |fi) if

p(x cj ) = Ÿ coeCq fi " c x cj c (3.1)
where 'c = (c 0 , . . . , c q≠1 ) oe C q = C 0 ◊ . . . ◊ C q≠1 :

ÿ coeCq fi c = 1 , " c x cj = I 1 if x 0 cj = c 0 , . . . , x q≠1 cj = c q≠1 0 otherwise .
Noting that the dependence structure between categorical features is modeled through Kronecker symbols (" c x cj ) coeCq , this dependence structure can be exploited to handle missing features such that the multivariate categorical distribution in (3.1) can be written as with fi c miss ,c obs the joint probability fi c defined in (3.1) for c = (c miss , c obs ). Outliers are not considered for categorical data since in our case only reliable categorical variables are filled in databases and unreliable ones are processed as missing data.

p(x miss cj , x obs cj ) = Ÿ c miss ,c obs oeC q miss j ◊C q obs j fi " c
p(x obs cj ) = ÿ x miss cj oeC q miss j p(x miss cj , x obs cj ) = Ÿ c obs oeC q obs j Q c c a ÿ c miss oeC q miss j fi c miss ,c obs R d d b " c obs x obs cj , p(x miss cj |x obs cj ) = p(x miss cj , x obs cj ) p(x obs cj ) = Ÿ c miss oeC q miss j Q c c c c c c a Ÿ c obs oeC q obs j Q c c c c a fi c miss ,c obs q c miss oeC q miss j fi c miss ,c obs R d d d d b " c obs x obs cj R d d d d d d b " c miss

Model

In this section, K emitters presenting mixed data are considered. Therefore, the main objective is to develop a mixture model which can build K distinct clusters even in presence of outliers and missing values. Before introducing a mixture model that handles mixed data, state-of-the-art mixture models for mixed data are reviewed. Then, component distributions for mixed data are defined and the mixture model is developed into a Bayesian framework.

State of the art

Two families of models emerge from finite mixture models fitting mixed-type data :

• The location mixture model [START_REF] Lawrence | Mixture separation for mixed-mode data[END_REF] that assumes that continuous variables follow a multivariate Gaussian distribution conditionally on both component and categorical variables.

• The underlying variables mixture model [START_REF] Everitt | A finite mixture model for the clustering of mixed-mode data[END_REF] that analyzes data sets with continuous and ordinal variables. It assumes that each discrete variable arises from a latent continuous variable and that all continuous variables (observed and unobserved) follow a Gaussian mixture model.

These two families are first detailed before introducing the retained approach.

Location Mixture Model

[LK96] introduced a location mixture model by assuming that the continuous variables are distributed as a finite mixture of Gaussians conditionally on the categorical variables. In other words, a Gaussian mixture exists for the continuous variables and its component mean vectors depend on the specific combination of categories modeled by the categorical variables. As pointed out by [WB99] the mixture of location models is not identifiable without imposing some constraints on the mean parameters of the Gaussian distributions. This is due to the indeterminacy of class memberships at each location. Even if all within component dependences are taken into account, we note that each combination of categories identifies a set of clusters. It follows that the total number of clusters can be unnecessarily large. A more parsimonious model is given by [START_REF] Hunt | Theory & methods: Mixture model clustering using the MULTIMIX program[END_REF], according to which the variables are decomposed into conditionally independent blocks containing a set of continuous variables or one categorical variable. This generally works as well as the within-independence assumption is realistic, and we know that there are cases where it is not. However, the local independence assumption represents a strong limitation, since it could lead to a solution with too many clusters, as shown by [START_REF] Jeroen | Latent class cluster analysis[END_REF]. By relaxing this assumption, a simpler solution with a lower number of groups can be obtained yielding a better classification. An even better classification can be reached by assuming di↵erent dependences in each group.

Underlying Variables mixture model [START_REF] Everitt | A finite mixture model for the clustering of mixed-mode data[END_REF] and [START_REF] Brian | The clustering of mixed-mode data: a comparison of possible approaches[END_REF] proposed a model according to which both the continuous and the categorical ordinal variables follow a homoscedastic Gaussian mixture model. However, as regards the ordinal variables, the mixture variables are only partially observed through their ordinal counterparts.

In other words, the ordinal variables are modeled following the Underlying Response Variable (UVR) approach. This satisfies the two main requirements: dealing with ordinal data properly and modeling dependences between ordinal and continuous variables. It is interesting to note that this model can be rewritten in terms of copulas [MBV17]. The main drawback of this model is that, in practice, it cannot be estimated through a full maximum likelihood approach, due to the presence of multidimensional integrals making the estimation time consuming. In sight of this, [Mor12] proposed a model-based clustering for mixed binary and continuous variables: each binary attribute is generated by dichotomizing a latent continuous variable, while the scores of the latent variables are estimated from the binary data. The estimated scores of the latent variables and the observed continuous data follow a multivariate Gaussian mixture model. Thus the estimation is carried out in two steps where the scores for binary data are firstly estimated before estimating the parameters of the mixture model. Eventually, [START_REF] Ranalli | Mixture models for mixed-type data through a composite likelihood approach[END_REF] proposed a model with no local independence or conditionally independent blocks assumption where the dependences between variables can be easily measured by adopting the URV approach for the ordinal variables and assuming that each component of the mixture follows a multivariate normal distribution such that the corresponding covariance matrices capture all the dependences regardless the nature of variables.

Retained approach

In this work, the location mixture model approach is retained since it better models relations between continuous and categorical radar features. Indeed, a radar pattern is mostly designed by first choosing a pattern of modulation features (categorical variables) to achieve a specific goal and then choosing continuous features (continuous variables) that meet constraints related to the chosen pattern and the tactical environment. Hence, continuous radar features are mainly chosen conditionally to categorical radar features and the location mixture model naturally responds to that dependence structure by assuming that continuous variables are normally distributed conditionally to categorical variables. Moreover, the local independence assumption proposed by [START_REF] Hunt | Theory & methods: Mixture model clustering using the MULTIMIX program[END_REF] is not retained in order to take advantage of the dependence structure between continuous and categorical data to infer on missing data.

Assumptions on mixed data

In this subsection, a joint distribution for mixed data is introduced to model the dependence structure between continuous and categorical data. Then, outliers and missing values are tackled by taking advantage of the joint distribution.

Distribution of mixed data

Considering that the retained approach focuses on conditioning continuous data x q = (x qj ) joeJ according to categorical data x c = (x cj ) joeJ , the following joint distribution is introduced

'j oe J , p(x qj , x cj ) = Ÿ coeCq (fi c N (x qj |µ c , )) " c x cj (3.5)
where continuous variables x qj are normally distributed according to categorical variables x cj with means (µ c ) coeCq and variance . As for categorical variables x cj , they are jointly distributed according to the multivariate categorical distribution defined in (3.1) and parametrized by fi = (fi c ) coeCq . Indeed, these conditional and marginal distributions can be obtained from (3.5) as follows :

'j oe J , p(x cj ) = ⁄ p(x qj , x cj )ˆx qj = Ÿ coeCq fi " c x cj c = MC(x cj |fi) , p(x qj |x cj ) = p(x qj , x cj ) p(x cj ) = Ÿ coeCq N (x qj |µ c , ) " c
x cj .

(3.6)

Outliers

As developed in Section 3.1, outliers are only considered for continuous data x q = (x qj ) joeJ and they are handled as in subsection 3.1.1 by introducing scale latent variables u = (u j ) joeJ . Nonetheless, the latent variables u are introduced conditionally to categorical data x c due to the dependence structure established in (3.5) and (3.6) such that 'j oe J ,

x qj |u j , x cj ≥ Ÿ coeCq N 1 x qj |µ c , u ≠1 j 2 " c x cj , u j |x cj ≥ Ÿ coeCq G (u j |-c , -c ) " c x cj ,
where each u j follows conditionally to categorical data x cj a Gamma distribution with rate and shape parameters

(-c , -c ) oe R ú+ ◊ R ú+ .

Missing Data

Both quantitative and categorical data (x qj , x cj ) joeJ can be partially observed. Hence (x qj , x cj ) joeJ are decomposed into observed features (x obs qj , x obs cj ) joeJ and missing features (x miss qj , x miss cj ) joeJ such that 'j oe J , with (x miss cj , x obs cj ) oe C q miss j ◊ C q obs j and q miss j + q obs j = q .

x
where (R d miss j , C q miss j ) and (R d obs j , C q obs j ), are disjoint subsets of (R d , C q ) embedding missing features (x miss qj , x miss cj ) and observed features (x obs qj , x obs cj ).

Missing continuous data x miss q = (x miss qj ) joeJ are nearly handled as in subsection 3.1.1 by taking advantage of properties of the multivariate normal distribution to obtain a distribution for missing values. The only di↵erence with the subsection 3.1.1 lies in the fact that continuous data are also distributed conditionally to categorical data x c due to the dependence structure established in (3.5) and (3.6). Hence, the following distributions are obtained 'j oe J ,

x miss qj |x obs qj , x cj ≥ Ÿ coeC N 3 x miss qj |µ x miss q jc , x miss q 4 " c x cj , x obs qj |x obs qj , x cj ≥ Ÿ coeC N 3 x obs qj |µ x obs q jc , x obs q 4 " c x cj
, where 'j oe J , 'c oe C q :

µ x miss q jc = µ miss c + cov obs ≠1 1 x obs qj ≠ µ obs c 2 , µ x obs q jc = µ obs c ,
x miss q = miss ≠ cov obs ≠1 cov' ,

x obs q = 3 obs ≠1 + 2 ◊ obs ≠1 cov' 1 x miss q 2 ≠1 cov obs ≠1 4

Proposed model

In this subsection, the retained approach is developed into a Bayesian framework where the proposed mixture model handles mixed-type data. Component distributions of clusters are first introduced before detailing the proposed model and its Bayesian framework.

Component Distribution

Assuming independent labels z = (z j ) joeJ for continuous and categorical observations x = (x qj , x cj ) joeJ and according to assumptions on mixed data defined in subsection 3.2.2, a component distribution for each cluster k oe K is obtained as follows

'j oe J , 'k oe K, p(x j |u j , z j = k) = p(x qj |u j , x cj , z j = k)p(x cj |z j = k) ,
where

p(x qj |u j , x cj , z j = k) = Ÿ coeCq N 1 x qj |µ kc , u ≠1 j k 2 " c x cj , p(x cj |z j = k) = Ÿ coeCq fi " c x cj kc .
Finally, the complete component distribution for each cluster k oe K results in

'j oe J , 'k oe K, p(x j , u j |z j = k) = Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c
x cj , (3.7) with

• u = (u j ) joeJ the scale latent variables handling outliers for quantitative data x q and distributed according to a Gamma distribution with shape and rate parameters (-, -) = (kc ,kc ) (k,c)oeK◊Cq ,

• (µ, ) = ((µ kc ) coeCq , k ) koeK the mean and the variance parameters of quantitative data

x q for each cluster, • fi = (fi k ) koeK the weights of the multivariate Categorical distribution of categorical data x c for each cluster.

K J " 0 0 µ 0 ÷ 0 k µ k fi 0 fi k -k -k u j x miss

Mixture model

Recalling that p(z j = k) = a k where a = (a k ) koeK are the weights related to component distributions, the mixture model is obtained from (3.7) such that 'j oe J ,

p(x j , u j | ) = ÿ koeK a k Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c x cj (3.8)
where = (a, fi, -, -, µ, ) is the set of parameters.

Bayesian framework

As in previous chapters, a Bayesian framework is used to estimate parameters . Assuming a dataset x = (x q , x c ) of i.i.d observations (x j = (x qj , x cj ) joeJ , independent labels z = (z j ) joeJ and scale latent variables u = (u j ) joeJ , the complete likelihood associated to (3.8) is defined by

p(x, z, u| , K) = Ÿ joeJ Ÿ koeK Q a a k Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c x cj R b " k z j
Eventually, the prior distribution required for is chosen as

p( |K) = p(

a|K)p(fi|K)p(-, -|K)p(µ, |K)

where

Y _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) , p(fi|K) = Ÿ koeK D(fi k |fi 0 ) , p(µ, |K) = Ÿ koeK Ÿ coeCq N 1 µ kc |µ 0 kc , ÷ ≠1 0 kc k 2 IW( k |" 0 , 0 ) , p(-, -|K) = Ÿ koeK Ÿ coeCq p(-kc , -kc |p 0 , q 0 , s 0 , r 0 ) .
and p(•, •|p, q, s, r) is the prior distribution defined in the previous chapter such that

'(-, -) oe R ú+ ◊ R ú+ , p(-, -|p, q, s, r) = 1 M p -≠1 e ≠q--s- (-) r (3.9)
where p, q, s, r > 0 and M =

Inference

Direct inference on the proposed model is not trivial since distributions of latent missing data and parameters may not be defined when both continuous and categorical features are missing. To overcome that issue, latent data and parameters are assumed to be independent a posteriori and their posterior distributions can be defined while keeping dependencies between parameters of these distributions. Therefore, the Variational Bayes (VB) procedure is processed to estimate parameters of the mixture model defined in (3.8). Variational posterior distributions are obtained from the VB Expectation (VBE) and VB Maximization (VBM) steps and a Lower Bound on the log evidence is defined to master the convergence of the VB procedure.

Variational posterior distributions

As previously, a factorized posterior distribution q(x miss q , u, x miss c , z, |K) = q(x miss q , u, x miss c , z|K)q( |K) is chosen as an approximation of the intractable posterior joint distribution p(x miss q , u, x miss c , z, |x obs q , x obs c , K) such that latent variables h = (x miss q , u, x miss c , z) and parameters are a posteriori independent and q(h|K) = q(x miss q |u, x miss c , z, K)q(u|x miss c , z, K)q(x miss c |z, K)q(z|K) , q( |K) = q(a|K)q(fi|K)q(-, -|K)q(µ, |K) .

According to VB assumptions, the following conjugate variational posterior distributions are obtained from the VB procedure

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ q(x miss q |u, x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC N 3 x miss qj | μx miss q jkc , u ≠1 j ˜ x miss q k 4 " c x cj " k z j , q(u|x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC G 1 u j |-jkc , -jkc 2 " c x cj " k z j , q(x miss c |z, K) = Ÿ joeJ Ÿ koeK MC 1 x miss cj |r x miss c jk 2 " k z j q(z|K) = Ÿ joeJ Cat(z j |r j ) , q(a|K) = D(a|Ÿ) , q(fi|K) = Ÿ koeK D(fi| fk ) , q(µ, |K) = Ÿ koeK Ÿ coeCq N 1 µ kc | μkc , ÷≠1 kc k 2 IW( k |" k , ˜ k ) , q(-, -|K) = Ÿ koeK Ÿ coeCq p(-kc , -kc |p k , qk , sk , rk ) , (3.10)
where the variational posterior distributions of (-, -) are defined in (3.9). Their respective parameters are estimated during the VBE and VBM-steps by developing expectations E [log p(x, u, z| , K)] and E h [log p(x, u, z, |K)].

VBE-step

The VBE-step consists in deriving the following expectation

E [log p(x, u, z| , K)] = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + ÿ coeCq " c x cj 3 E [log fi kc ] ≠ 1 2 3 d(log 2fi ≠ log u j ) + E [log | k |] + u j E Ë (x qj ≠ µ kc ) T ≠1 k (x qj ≠ µ kc ) È 4 + E [-kc ] E [log -kc ] + (E [-kc ] ≠ 1) log u j ≠ E [log (-kc )] ≠ u j E [-kc ] 44 (3.11)
where

'(j, k, c) oe J ◊ K ◊ C q : E Ë (x qj ≠ µ kc ) T ≠1 k (x qj ≠ µ kc ) È = (x qj ≠ μkc ) T "k ˜ ≠1 k (x qj ≠ μkc ) + d ÷kc (3.12)
is obtained from properties of the variational distribution q(µ, |K) in (3.10). Hence quantitative data x q are distributed a posteriori according to a product of normal distributions conditionally to categorical data x c , latent variables u and labels z such that

x q |u, x c , z ≥ Ÿ joeJ Ÿ koeK Ÿ coeC N 1 x qj | μkc , u ≠1 j "≠1 k ˜ k 2 " c x cj " k z j
Mean parameters ( μkc ) (k,c)oeK◊Cq and variance parameters (" ≠1 k ˜ k ) koeK of these normal distributions are obtained from (3.12). By decomposing x q into (x miss q , x obs q ) and by exploiting properties of the multivariate normal distribution, the following variational posterior distribution is obtained for missing values x miss q : q(x miss

q |u, x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC N 3 x miss qj | μx miss q jkc , u ≠1 j ˜ x miss q k 4 " c x cj " k z j with '(j, k, c) oe J ◊ K ◊ C q : μx miss q jkc = μmiss kc + ˜ cov k ˜ obs ≠1 k (x obs qj ≠ μobs kc ) , ˜ x miss q k = ˜ miss k ≠ ˜ cov k ˜ obs ≠1 k ˜ cov Õ k "k .
Then by marginalising over x miss q in (3.11), the expectation (3.11) becomes

⁄ E [log p(x, u, z| , K)] ˆxmiss q = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + ÿ coeCq " c x cj 3 E [log fi kc ] ≠ 1 2 3 d obs j (log 2fi ≠ log u j ) + E [log | k |] ≠ log | ˜ x miss q k | + u j 3 1 x obs qj ≠ μobs kc 2 T ˜ x obs qj ≠1 k 1 x obs qj ≠ μobs kc 2 + d ÷kc 4 4 + E [-kc ] E [log -kc ] + (E [-kc ] ≠ 1) log u j ≠ E [log (-kc )] ≠ u j E [-kc ] 44 (3.13)
with 'k oe K,

˜ x obs qj k = A ˜ obs ≠1 k + 2 ◊ ˜ obs ≠1 k ˜ cov Õ k 3 ˜ x miss q k 4 ≠1 ˜ cov k ˜ obs ≠1 k B ≠1
"k .

Conditionally to x c and z, the scale latent variables u are distributed according to a product of Gamma distribution whose parameters are obtained by aggregating terms related to u such that

q(u|x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC G 1 u j |-jkc , -jkc 2 " c x cj " k z j with '(j, k, c) oe J ◊ K ◊ C q : -jkc = E [-kc ] + d obs j 2 , -jkc = E [-kc ] + 1 2 3 1 x obs qj ≠ μobs kc 2 T ˜ x obs qj ≠1 k 1 x obs qj ≠ μobs kc 2 + d ÷kc 4 .
Then by marginalising over u in (3.13), the expectation (3.13) becomes

⁄ E [log p(x, u, z| , K)] ˆxmiss q ˆu = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + ÿ coeCq " c x cj 3 E [log fi kc ] ≠ 1 2 3 d obs j log 2fi + E [log | k |] ≠ log | ˜ x miss q k | 4 + E [-kc ] E [log -kc ] ≠ E [log (-kc )] + log (-jkc ) ≠ -jkc log -jkc 44 = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + ÿ coeCq " c x cj log fl x cj kc 4 (3.14)
where '(j, k, c) oe J ◊ K ◊ C q , log fl

x cj kc = E [log fi kc ] ≠ 1 2 3 d obs j log 2fi + E [log | k |] ≠ log | ˜ x miss q k | 4 + E [-kc ] E [log -kc ] ≠ E [log (-kc )] + log (-jkc ) ≠ -jkc log -jkc . (3.15)
By decomposing each x cj oe C q into (x miss cj , x obs cj ) oe C q miss j ◊ C q obs j and c oe C q into (c miss , c obs ) oe C q miss j ◊ C q obs j , the marginalised expectation (3.14) can be developed as Hence the variational categorical distributions are deduced from (3.18) and are given by

⁄ E [log p(x, u, z| , K)] ˆxmiss q ˆu = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + ÿ c miss oeC q miss j " c miss
q(z|K) = Ÿ joeJ Cat(z j |r j )
where probabilities (r j ) joeJ are obtained from (3.19) such that 'j oe J , k oe K, rjk = fl jk q koeK fl jk .

VBM-step

The VBM-step consists in deriving the following expectation

E h [log p(x, u, z, |K)] = E h [log p(x, u, z| , K)] + p( |K) = ÿ joeJ ÿ koeK E h Ë " k z j È 3 log a k + ÿ coeCq E h Ë " c x cj È 3 log fi kc ≠ 1 2 3 log | k | + d(log 2fi ≠ E h [log u j ]) + E h Ë u j (x qj ≠ µ kc ) T ≠1 k (x qj ≠ µ kc ) È 4 + -kc log -kc + (-kc ≠ 1)E h [log u j ] ≠ log (-kc ) ≠ E h [u j ] -kc 44 + ÿ koeK (Ÿ 0 k ≠ 1) log a k + log c D (Ÿ 0 ) ≠ 1 2 3 (" 0 + d + 1) log | k | + Trace 1 0 ≠1 k 2 4 + c IW (" 0 , 0 ) + ÿ coeCq 1 2 3 d(log ÷ 0 kc ≠ log 2fi) ≠ log | k | ≠ ÷ 0 kc 1 µ kc ≠ µ 0 kc 2 T ≠1 k 1 µ kc ≠ µ 0 kc 2 4 ≠ log M 0 + (-kc ≠ 1) log p 0 ≠ r 0 log (-kc ) + s 0 -kc log -kc ≠ q 0 -kc + (fi 0 kc ≠ 1) log fi kc + log c D (fi 0 k ) , (3.20)
where

'(j, k, c) oe J ◊ K ◊ C q : E h Ë u j (x qj ≠ µ kc ) T ≠1 k (x qj ≠ µ kc ) È = E h [u j ] (E h [x qj ] ≠ µ kc ) T ≠1 k (E h [x qj ] ≠ µ kc ) + Trace 1 V h [x qj ] ≠1 k 2 (3.21)
is obtained from properties of the variational distribution q(h|K) with

E h [u j ] = -jkc -jkc , E h [x qj ] = A μx miss q jkc x obs qj B , V h [x qj ] = Q a ˜ x miss q k 0 d miss j ◊d obs j 0 d obs j ◊d miss j 0 d obs j ◊d obs j R b .
By factorizing terms related to a in (3.20), the following Dirichlet distribution is obtained

q(a|K) = D(a|Ÿ)
where

'k oe K, Ÿk = Ÿ 0 k + ÿ joeJ E h Ë " k z j È .
Like the variational distribution of a, variational posterior distributions of fi are obtained by factorizing terms related to fi in (3.20) and are given by

q(fi|K) = Ÿ koeK D(fi k | fk ) where '(k, c) oe K ◊ C q , fkc = fi 0 kc + ÿ joeJ E h Ë " k z j È E h Ë " c x cj È . with E h Ë " c x cj È
is obtained by decomposing categorical features into observed and missing features such that

E h Ë " c x cj È = " c obs x obs cj ◊ E h 5 " c miss x miss cj 6
.

Then, (-, -) are a posteriori distributed according to the distribution defined in (3.9) such that

q(-, -|K) = Ÿ koeK Ÿ coeCq p(-kc , -kc |p kc , qkc , skc , rkc ) ,
where

pkc = p 0 exp Q a ÿ joeJ E h Ë " k z j È E h Ë " c x cj È E h [log u j ] R b , qkc = q 0 + ÿ joeJ E h Ë " k z j È E h Ë " c x cj È E h [u j ] , skc = s 0 + ÿ joeJ E h Ë " k z j È E h Ë " c x cj È , rkc = r 0 + ÿ joeJ E h Ë " k z j È E h Ë " c x cj È .
By aggregating and factorizing terms related to each µ kc in (3.20), a Normal distribution is obtained for each µ kc such that

q(µ| , K) = Ÿ koeK Ÿ coeCq N 1 µ kc | μkc , ÷≠1 kc k 2
where 'k oe K and 'c oe

C q ÷kc = ÿ joeJ E h Ë " k z j È E h Ë " c x cj È E h [u j ] + ÷ 0 kc , μkc = q joeJ E h Ë " k z j È E h Ë " c x cj È E h [u j ] E h [x qj ] + ÷ 0 kc µ 0 kc ÷kc .
Eventually, variance parameters are a posteriori distributed according to Inverse Wishart distributions given by q

( |K) = Ÿ koeK IW( k |" k , ˜ k )
where

"k = " 0 + ÿ joeJ E h Ë " k z j È , ˜ k = 0 + ÿ joeJ ÿ coeCq E h Ë " k z j È E h Ë " c x cj È E h [u j ] E h [x qj ] E h [x qj ] T + ÿ joeJ E h Ë " k z j È V h [x qj ] + ÿ coeCq ÷ 0 kc µ 0 kc µ T 0 kc ≠ ÷kc μkc μT kc .

Lower Bound

Recalling that the Lower Bound on the log evidence is given by

L(q|K) = E h, [log p(x, h, |K)] ≠ E h, [log q(h, |K)] (3.22)
where E h, [log p(x, h, |K)] is the free energy and E h, [log q(h, |K)] is the entropy of the approximate posterior q(h, |K). The free energy can be developed as

E h, [log p(x, h, |K)] = E h, [log p(x, h| , K)] + E [log p( |K)]
where

E h, [log p(x, h| , K)] = ÿ joeJ ÿ koeK E h Ë " k z j È 3 E [log a k ] + ÿ coeCq E h Ë " c x cj È 3 E [log fi kc ] ≠ 1 2 3 d(log 2fi ≠ E h [log u j ]) + E [log | k |] + E h, Ë u j (x qj ≠ µ kc ) T ≠1 k (x qj ≠ µ kc ) È 4 + E [-kc ] E [log -kc ] + (E [-kc ] ≠ 1)E h [log u j ] ≠ E [log (-kc )] ≠ E h [u j ] E [-kc ] 44 with E h, Ë u j (x qj ≠ µ kc ) T ≠1 k (x qj ≠ µ kc ) È = E h [u j ] 3 (E h [x qj ] ≠ μkc ) T "k ˜ ≠1 k (E h [x qj ] ≠ μkc ) d ÷kc 4 + Trace 1 V h [x qj ] "k ˜ ≠1 k 2
and

E [log p( |K)] = ÿ koeK (Ÿ 0 k ≠ 1)E [log a k ] + log c D (Ÿ 0 ) ≠ 1 2 3 (" 0 + d + 1)E [log | k |] + Trace 1 0 E Ë ≠1 k È2 4 + c IW (" 0 , 0 ) + ÿ coeCq 1 2 3 d(log ÷ 0 kc ≠ log 2fi) ≠ E [log | k |] ≠ ÷ 0 kc E 5 1 µ kc ≠ µ 0 kc 2 T ≠1 k 1 µ kc ≠ µ 0 kc 2 6 4 ≠ log M 0 + (E [-kc ] ≠ 1) log p 0 ≠ r 0 E [log (-kc )] + s 0 E [-kc ] E [log -kc ] ≠ q 0 E [-kc ] + (fi 0 kc ≠ 1)E [log fi kc ] + log c D (fi 0 ) . with E 5 1 µ kc ≠ µ 0 kc 2 T ≠1 k 1 µ kc ≠ µ 0 kc 2 6 = 1 μkc ≠ µ 0 kc 2 T "k ˜ ≠1 k 1 μkc ≠ µ 0 kc 2 + d ÷kc .
As for the entropy term, the following decompositon is obtained

E h, [log q(h, |K)] = E h Ë log q(x miss q , u, x miss c , z|K) È + E [log q( |K)] = E h Ë log q(x miss q |u, x miss c , z, K) È + E h Ë log q(u|x miss c , z, K) È + E h Ë log q(x miss c |z, K) È + E h [log q(z|K)] + E [log q( |K)]
where

E h Ë log q(x miss q |u, x miss c , z, K) È = ÿ joeJ ÿ koeK ÿ coeCq E h Ë " k z j È E h Ë " c x cj È 1 2 3 d miss j (E h [log u j ] ≠ log 2fi ≠ 1) ≠ log | ˜ x qj k | 4 , E h Ë log q(u|x miss c , z, K) È = ÿ joeJ ÿ koeK ÿ coeCq E h Ë " k z j È E h Ë " c x cj È 3 -jkc log -jkc ≠ log (-jkc ) + (-jkc ≠ 1)E h [log u j ] ≠ -jkc E h [u j ] 4 , E h Ë log q(x miss c |z, K) È = ÿ joeJ ÿ koeK ÿ c miss oeC q miss j E h Ë " k z j È E h 5 " c miss x miss cj 6 log rx miss c jkc miss E h [log q(z|K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log rjk
and

E [log q( |K)] = ÿ koeK (Ÿ k ≠ 1)E [log a k ] + log c D (Ÿ k ) ≠ 1 2 3 (" k + d + 1)E [log | k |] + Trace 1 ˜ k E Ë ≠1 k È2 4 + c IW (" k , ˜ k ) + ÿ coeCq 1 2 3 d(log ÷kc ≠ log 2fi) ≠ E [log | k |] ≠ ÷kc E Ë (µ kc ≠ μkc ) T ≠1 k (µ kc ≠ μkc ) È 4 ≠ log M k + (E [-kc ] ≠ 1) log pkc ≠ rkc E [log (-kc )] + skc E [-kc ] E [log -kc ] ≠ qkc E [-kc ] + (f kc ≠ 1)E [log fi kc ] + log c D ( fk ) , with E Ë (µ kc ≠ μkc ) T ≠1 k (µ kc ≠ μkc ) È = d ÷kc .

Expectations from variational distributions

Expectations developed in variational calculations are derived from properties of variational posterior distributions and are obtained as follows. Categorical distribution properties lead to

'j oe J , 'k oe K, 'c miss oe C q miss j : E h Ë " k z j È = rjk , E h 5 " c miss x miss cj 6 = rx miss c jkc miss .
Dirichlet distribution properties lead to

'k oe K, 'c oe C q : E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b , E [log fi kc ] = Â(f kc ) ≠ Â Q a ÿ coeCq fkc R b ,
where Â(•) is the digamma function. Gamma distribution properties lead to

'j oe J , 'k oe K, 'c oe C q : E h [u j ] = -jkc -jkc , E h [log u j ] = Â(-jkc ) ≠ log -jkc .
Normal distribution properties lead to

'k oe K, 'c oe C q : E [µ kc ] = μkc , E Ë µ kc µ T kc È = V [µ kc ] + E [µ kc ] E [µ kc ] T = ÷≠1 kc k + μkc μT kc ,
Inverse Wishart distribution properties lead to

E Ë ≠1 k È = "k ˜ ≠1 k , E [log | k |] = log | ˜ k | ≠ d ÿ i=1 Â 3 "k + 1 ≠ i 2 4 ≠ d log 2 .
Posterior expectations ofcan easily be computed conditionally tosuch that '(k, c) oe K ◊ C q :

E [-kc ] = skc E [-kc ] + 1 qkc , E [log -kc ] = E [Â (s kc -kc + 1)] ≠ log qkc .
However, expectations depending onkc are intractable

E [Â (s kc -kc + 1)] = ⁄ Â (s kc -kc + 1) p(-kc |p kc , rkc )d-kc , E [-kc ] = ⁄ -kc p(-kc |p kc , rkc )d-kc , E [log (-kc )] = ⁄ log (-kc )p(-kc |p kc , rkc )d-kc .
As in previous chapter, the deterministic method introduced by [TK86] is applied to estimate these expectations.

Experiments

In this section, the proposed method is performed on 3 sets of realistic simulated data which are composed of continuous, categorical and mixed data. For comparison, a standard neural network (NN), the k-nearest neighbours (KNN) algorithm, Random Forests (RdF) the k-means algorithm and the DBSCAN are also evaluated. Two experiments are carried out to evaluate classification and clustering performance with respect to a range of percentages of missing values. First, characteristics for realistic data acquisition and imputation methods for missing data are detailed. Then, both experiments are described with their error measure and their performance are shown to exhibit the e↵ectiveness of the proposed model.

Data

Realistic data are generated from an operational database gathering 55 radar emitters presenting various patterns. Each pattern consists of a sequence of pulses which are defined by a triplet of continuous features x q = (RF,PW,PRI) and a fourtet of categorical features x c referring to pulse-to-pulse modulations of RF and PRI, intrapulse modulations of RF and scanning types. 42 combinations of the categorical features are observed among the 55 emitters and they are composed of modulations listed in Table 3.2. For each radar emitter, 100 observations (x j ) 100 j=1 are simulated from its pattern of pulses such that an observation x j = (x qj , x cj ) is made up of continuous features x qj and categorical features x cj related to one of the pulses. Then, continuous observations are noised by applying a multivariate Gaussian noise with a diagonal covariance matrix whose diagonal elements are [ ‡ 2 RF , ‡ 2 P W , ‡ 2 P RI ] = [1MHz, 50ns, 1µs]. Negative features issued from the generated noise are thresholded to zero. Hence, outliers are embedded in observations due to the thresholding step. The dataset is shown in Figure 3.7. Moreover, extra missing values are added to evaluate limits of the proposed approach by randomly deleting coordinates of (x qj ) 100 j=1 and (x cj ) 100 j=1 for each of the 55 radar emitters. Percentages of deletion range from 5% to 90%. At last the continuous dataset, categorical dataset and the mixed dataset are composed of (x qj ) 5500 j=1 , (x cj ) 5500 j=1 and (x j ) 5500 j=1 = (x qj , x cj ) 5500 j=1 . As in the previous chapter, imputation methods [GLSGFV10] are used to handle missing data for comparison algorithms. Mean and k-nearest neighbours imputation methods are still implemented for continuous data. Regarding missing categorical data, they are handled through the k-nearest neighbours and mode imputation methods. The mode imputation consists in filling a missing component of an observation by the mode of observed values of that component. This method has the obvious disadvantage that it under represents the variability and also ignores correlations between observations [Sch97]. These imputation methods are compared with the proposed approach in terms of classification, clustering and reconstruction performance. For the comparison of reconstruction performance on continuous data, mean-squared errors between original continuous data and previous imputation methods are compared with the mean-squared error between original continuous data and the variational posterior marginal mean of missing continuous data given by 'j oe J , xmiss (3.24)

Classification experiment

The classification experiment evaluates the ability of each algorithm to assign unlabeled data to one of the K classes trained by a set of labeled data. As in the previous chapter, the classification task is decomposed into a training step and a prediction step defined in procedures 3.1 and 3.2.

The training step consists in estimating variational parameters of q( ) given a set of training data with known labels. As for the prediction step, it results in associating new data to the class that maximizes their posterior probabilities. Ÿ 0 fi 0 ÷ 0 " 0 p 0 r 0 q 0 s 0 µ 0c 0 0.5 0.5 10 ≠4 4 0.9 1 1 1 [0, 0, 0] I 3 Procedure 3.1 Classification procedure on mixed data : Training step Input: Training set x train and associated labels z train Output: Learned parameters ˜ train Initialise Ÿ 0 , fi 0 , " 0 , ÷ 0 , µ 0 , 0 , p 0 , r 0 , s 0 and q 0 for iter = 1 to itermax do Update -jkc , -jkc ,

rx miss cj jkc , μx miss q jkc , ˜ x miss qj k Update Ÿk , ÷kc , "k , pkc , rkc , skc , qkc , fkc , μkc , ˜ k Calculate the lower bound L if L iter ≠ L iter≠1 AE tol ◊ L iter≠1 then return ˜ train = 1 Ÿk , ÷kc , "k , pkc , rkc , skc , qkc , fkc , μkc , ˜ k 2 (k,c)oeK◊Cq

end if end for

For the classification experiment, results are shown in Figure 3.8 where classification performance are exhibited for the 3 datasets. Without missing data, both algorithms cannot perfectly classify the 55 radar emitters for the 3 datasets. Indeed, both algorithms reach accuracies of 90% for the continuous dataset, 75% for the categorical dataset and 98% for the mixed dataset. These performance can be explained by the non total separability of continuous and categorical datasets since the 55 emitters share 42 combinations of categorical features (Table 3.2) and (RF,PRI,PW) intervals as shown in Figure 3.7. Nonetheless when mixed data are taken into consideration, the dataset becomes more separable leading to higher performance of both algorithms. When the proportion of missing values increases, the proposed model outperforms comparisons algorithms for each dataset. It achieves accuracies of 80%, 55% and 95% for 90% of deleted continuous, categorical and mixed values whereas accuracies of comparison algorithms are lower than 65%, 50% and 75% with missing data imputation from standard methods. As in the previous chapter, these higher performance of the proposed model reveal that the proposed method embeds a more e cient inference method than other imputation methods. That result is confirmed on Figure 3.8 when comparison algorithms are applied on data reconstructed by the proposed model. Indeed when the proposed inference is chosen, comparison algorithms share the same performance than the proposed model and manage to handle missing even for 90% of deleted values. Finally, this e ciency is shown on Figure 3.9 where data reconstructed by the proposed model exhibit lower mean-squared errors and Jaccard distances for missing data imputation than the standard imputation methods. Indeed, the lowest mean-squared errors and Jaccard distances are obtained by the proposed model reconstruction on the mixed dataset, which demonstrate that missing data imputation is even more e cient when both continuous and categorical are jointly modeled. Furthermore, a correlation between higher performance of the proposed model and the quality of its reconstructions can be noticed for any percentage of missing values. Then, e↵ectiveness of the proposed model can be explained by the fact that missing data imputation methods can create outliers that deteriorate performance of classification algorithms whereas the inference on missing data and labels prediction are jointly estimated in the proposed model. Indeed, embedding the inference procedure into the model framework allows properties of the model, such as outliers handling, to counterbalance drawbacks of imputation methods such as outlier creation.

Clustering experiment

The clustering experiment is composed of two experiments that aim to exhibit the clustering ability of each algorithm according to an a priori number of clusters K oe {K min , . . . , K max }. As developed in previous chapters, the clustering algorithm is decomposed into two parts. First, a semi-supervised classification is performed for each K ranges from K min to K max to estimate variational parameters of q( , H) and labels in a mixture of K components. Then, the value of K that maximizes the lower bound is retained as the posterior number of clusters as well as its associated parameters.

Procedure 3.3 Semi-supervised classification procedure on mixed data Input: Unlabelled dataset x and number of classes K Output: Labels z and parameters ˜ Initialise Ÿ 0 , fi 0 , " 0 , ÷ 0 , µ 0 , 0 , p 0 , r 0 , s 0 and q 0 for iter = 1 to itermax do Update -jkc , -jkc , Ÿ 0 fi 0 ÷ 0 " 0 p 0 r 0 q 0 s 0 0.5 0.5 100 4 0.9 1 1 1 order to evaluate the impact of the a priori number of clusters on data clustering. Parameters of DBSCAN are set to Minpts = 5 and eps = 0.01 by using an heuristic proposed in the original paper [EKS + 96]. A supervised initialisation is retained for the proposed model due to its sensitivity to initialisation. It consists in initialising prior component means µ 0c from results of a k-means algorithm and prior component covariance matrices 0 from diagonal matrices whose diagonal elements are variances of observed features. Other hyper-parameters are initialised as in Table 3.4. Since comparison algorithms do not handle observations with missing values and do not provide a clustering result for them, missing data are reconstructed through the mean, the KNN and the proposed model imputation methods before running these algorithms.

rx miss cj jkc , μx miss q jkc , ˜ x miss qj k , rjk Update Ÿk , ÷kc , "k , pkc , rkc , skc , qkc , fkc , μkc , ˜ k Calculate the lower bound L if L iter ≠ L iter≠1 AE tol ◊ L iter≠1 then return ˜ = 1 Ÿk , ÷kc , "k , pkc , rkc , skc , qkc , fkc , μkc , ˜ k 2 (k,c
The first clustering experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K oe {K min , . . . , K max }. Performance are evaluated through the Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF] that compares estimated partitions of data with the ground-truth. Results of the first experiment on the 3 datasets are shown in Figures 3.10 and 3.11. Without the presence of missing values, performance of DBSCAN, k-means and the proposed model are similar with ARIs of 62%, 72% and 88% for the continuous, categorical and mixed datasets. As in the classification experiment, these performance are explained by the non total separability of continuous and categorical datasets since the 55 emitters share 42 combinations of categorical features (Table 3.2) and (RF,PRI,PW) intervals as shown in Figure 3.7. Once again when mixed data are taken into consideration, the dataset becomes more separable leading to higher performance of both algorithms (ARI = 88%). When the proportion of missing values increases, the proposed model outperforms both DBSCAN and k-means and achieves ARIs of 35% ,40% and 58% on continuous, categorical and mixed datasets for 90% of deleted values whereas the ARIs of comparison algorithms with standard missing data imputation are lower than 35% on each dataset. As in the classification experiment, these higher performance reveal that the proposed method embeds a more e cient inference method than other imputation methods. That result is confirmed on both Figure 3.10 and Figure 3.11 where DBSCAN and k-means are applied on data reconstructed by the proposed model. Indeed, performance of both algorithms increase up to performance of the proposed model for any percentages of deleted values when the proposed inference is chosen.

The second experiment tests the ability of each algorithm to find the true number of clusters K among {K min , . . . , K max }. The lower bound (3.22) and the average Silhouette score [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF] are criteria used to select the optimal number of clusters for the proposed model and the k-means algorithm. Indeed, the ARI cannot be used since it requires the ground-truth and DBSCAN automatically selects a number of clusters for a given dataset. Results of the second experiment on the 3 datasets are visible on Figures Figure 3.12 and 3.13. Figure 3.12 presents numbers of clusters selected by the lower bound and average Silhouette scores for the proposed model and k-means algorithm according to di↵erent proportions of missing values and imputation methods. Without missing data, the correct number of clusters (K=55) is selected by the two criteria for the k-means algorithm and the proposed model when continuous and mixed data are clustered. As for categorical data, both criteria select 45 as the optimal number of clusters since the 42 combinations of categorical features (Table 3.2) shared by the 55 emitters constitute 42 distinct clusters. In presence of missing values, the average Silhouette score mainly selects K = 65 when the k-means algorithm is run on the 3 datasets completed by standard imputation methods. When, the k-means algorithm performs clustering on the posterior reconstructions, the average Silhouette score correctly selects K = 55 until 60% of missing values for continuous data and 40% of missing values for mixed data. Eventually when the proposed model does clustering, the two criteria select the correct number of clusters K = 55 until 70% of missing values for continuous and mixed data. These results show two main advantages of the proposed model. As previously, the proposed model provides a more robust inference on missing data since the average Silhouette score chooses more representative number of clusters when the k-means algorithm is run on the posterior reconstructions than on data completed by standard imputation methods. Furthermore, since the lower bound criterion also selects the correct number of clusters as the average Silhouette score, it can be used as a valid criterion for selecting the optimal number of clusters and does not require extra computational costs as the Silhouette score since it is computed during the model parameter estimation. Finally, the proposed approach provides a more robust inference on missing data and a criterion for selecting the optimal number of clusters without extra computations. As for the Figure 3.13, it shows the evolution of the number of clusters estimated by DBSCAN according to di↵erent proportions of missing values and imputation methods. Since DBSCAN automatically estimates the number of clusters and manages outliers by creating new clusters, results on Figure 3.13 can be used to evaluate performance of imputations methods. For mean, mode and k-NN imputation methods, DBSCAN estimates a number of clusters greater than the number estimated for the proposed model according to any proportion of missing values. These performance indicate that the proposed approach creates less outliers than other imputation methods by providing a more robust inference on missing data since DBSCAN localizes less outliers in the posterior reconstructions (3.23-3.24) than in standard imputation methods.

Conclusion

In this chapter, various modulations of radar emitter patterns have been presented. These modulations can be used as categorical data in the classification and clustering tasks. Hence, a mixture model handling both continuous data and categorical data has been developed. An approach based on the Location Mixture Model has been investigated by establishing conditional relations between continuous and categorical data. Benefiting from a dependence structure designed for mixed data, the proposed model shows its e ciency for inferring on missing data, performing classification and clustering tasks and selecting the correct number of clusters. Since the posterior distribution is intractable, model learning is processed through a variational Bayes inference where variational posterior distributions are proposed for continuous and categorical missing values. Experiments show that the proposed approach handles mixed data even in presence of missing values and can outperform standard algorithms in clustering tasks. Indeed the main advantage of our approach is that it enables the counterbalance of imputation methods drawbacks by embedding the inference procedure into the model framework.

Chapter 4

Temporal Evolution Data

Continuous data describing radar emitters waveforms such as the Carrier Frequency, the Pulse Width and the Pulse Repetition Interval have been previously taken into account in order to cluster radar emitters. Nonetheless, the Pulse Description World gathers other features such the Amplitude whose relation with the Time of Arrival reflects the scanning behaviour of a radar emitter. Therefore, this temporal relation can be exploited to cluster radar emitters. Depending on the scanning type, this relation can be represented by either a parabola or a piecewise parabola. These two relations have to be included into the mixture distribution to take advantage of the temporal behaviour of each radar emitter. This chapter contains three sections which focus on three di↵erent cases. The first section deals with data where only emitters having a parabolic scanning behaviour are observed. The second section introduces the case where temporal evolution data are only distributed according to piecewise parabolic relations. As for the last section, it is about the case where any type of scanning behaviours can be observed in data. Each section presents the model integrating radar temporal evolution data and its inference procedure before proposing a more complete model taking into consideration temporal evolution data and mixed data. Eventually, experiments are carried out to exhibit performance of the proposed approach. In this chapter, radar temporal evolution data consist of J pulses gathering J amplitudes x t = (x tj ) joeJ and J times of arrival t = (t j ) joeJ from K distinct emitters.

Parabolic data

In this section, K emitters presenting parabolic scanning behaviours are considered. Therefore, the main objective is to develop a mixture model which can build K distinct clusters formed by K parabolas and to define an inference procedure for parameter estimation. Then, the proposed model is enhanced with the mixture model designed for mixed data in Chapter 3 in order to improve clustering performance. Finally, experiments on synthetic and real data are carried out to exhibit performance of the proposed approach.

Model

Before introducing a mixture model that handles parabolic data, the parabolic relation between amplitudes and times of arrival is defined. Then, the mixture model is developed into a Bayesian framework.

Parabola Equation

The parabolic relation between the amplitude x tj of the j th pulse and its time of arrival t j , visible on Figure 4.1, can be described by the following parabolic equation

x tj = at 2 j + bt j + c + ' (4.1)
where

' ≥ N (0, ‡ 2 )
x tj = (t j ) T Ê + ' (4.2)
with (t j ) = (t 2 j , t j , 1) T the vector containing polynomial transformations of t j and Ê = (a, b, c) T the vector of regression parameters. Since the measurement error is assumed to be Gaussian, the amplitude x tj is distributed according to a normal distribution centered in (t j ) T Ê with variance ‡ 2

x tj ≥ N

1 x tj | (t j ) T Ê, ‡ 2 2 . (4.3)

Mixture model

Since each radar emitter has its own scanning behaviour, K unique parabolas exist in data and they are configured with K regression parameters Ê = (Ê k ) koeK . Then, each amplitude x tj belongs to one of these parabolas which is related to a specific emitter. In other words, conditionally to its label z j and its time of arrival t j , the amplitude x tj is distributed according to (4.3) such that the component distribution is defined by

x tj |t j , z j = k ≥ N 1 x tj | (t j ) T Ê k , ‡ 2 2 (4.4)
Recalling that p(z j = k) = a k where a = (a k ) koeK are the weights related to component distributions, the mixture model is obtained from (4.4) such that

'j oe J , p(x tj |t j , ) = ÿ koeK a k N 1 x tj | (t j ) T Ê k , ‡ 2 2 (4.5)
where = (a, Ê, ‡ 2 ) is the set of parameters.

Bayesian framework

As in chapters 2 and 3, a Bayesian framework is used to estimate parameters . Assuming datasets (x t , t) of i.i.d observations (x tj , t j ) joeJ and independent labels z = (z j ) joeJ , the complete likelihood associated to (4.5) is defined by

p(x t , z|t, , K) = Ÿ joeJ Ÿ koeK 1 a k N 1 x tj | (t j ) T Ê k , ‡ 2 22 " k z j
Eventually, the prior distribution required for is chosen as

p( |K) = p(a|K)p(Ê| ‡ 2 , K)p( ‡ 2 )
where a follows a Dirichlet distribution, each Ê k follows a Normal distribution and ‡ 2 follows an Inverse Gamma distribution such that

Y _ _ _ _ ] _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) , p(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k |Ê 0 , ‡ 2 0 2 , p( ‡ 2 ) = IG( ‡ 2 |› 0 1 , › 0 2 ) .
The resulting mixture model is shown on Figure 4.2.

Inference

The Variational Bayes (VB) procedure is derived to estimate parameters of the mixture model defined in (4.5). Variational posterior distributions are obtained from the VB Expectation (VBE) and VB Maximization (VBM) steps and a lower bound on the log evidence is defined to master the convergence of the VB procedure.

K J Ê 0 0 Ê k x tj (t j ) z j a Ÿ 0 ‡ 2 › 0 1 › 0 2 Figure 4.2:
Graphical representation of the proposed mixture model handling parabolic data. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x tj , t j ) decomposed into the amplitude x tj and the polynomial transformation (t j ) and the indicator variables z j . Known quantities, respectively unknown quantities, are in blue, respectively in red.

Variational posterior distributions

As previously, a factorized posterior distribution q(z, |K) = q(z|K)q( |K) is chosen as an approximation of the intractable posterior joint distribution p(z, |x t , t, K) such that latent variables z and parameters are a posteriori independent and q( |K) = q(a|K)q(Ê| ‡ 2 K)q( ‡ 2 ).

According to VB assumptions, the following conjugate variational posterior distributions are obtained from the VB procedure

Y _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ [ q(z|K) = Ÿ joeJ Cat(z j |r j ) , q(a|K) = D(a|Ÿ) , q(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k | Ễk , ‡ 2 ˜ k 2 , q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 ) . (4.6)
Their respective parameters are estimated during the VBE and VBM steps.

VBE-step

The VBE-step consists in deriving the following expectation

E [log p(x t , z|t, , K)] = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 D 44 = ÿ joeJ ÿ koeK " k z j log fl jk (4.7)
where

log fl jk = E [log a k ] ≠ 1 2 A log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 DB . (4.8)
Hence, a categorical distribution for labels z is deduced from (4.7) such that

q(z|K) = Ÿ joeJ Cat(z j |r j )
and their parameters (r j ) joeJ are obtained from (4.8) as follows 'j oe J , 'k oe K, rjk = fl jk q koeK fl jk .

VBM-step

The VBM-step consists in deriving the following expectation

E z [log p(x t , z, |t, K)] =E h [log p(x t , z|t, , K)] + log p( |K) = ÿ joeJ ÿ koeK E h Ë " k z j È 3 log a k ≠ 1 2 3 log 2fi + log ‡ 2 + (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 44 + ÿ koeK (Ÿ 0 k ≠ 1) log a k ≠ 1 2 A 3(log 2fi + log ‡ 2 ) + log | 0 | + (Ê k ≠ Ê 0 ) T ≠1 0 ‡ 2 (Ê k ≠ Ê 0 ) B ≠ (› 0 1 + 1) log ‡ 2 ≠ › 0 2 ‡ 2 + log c D (Ÿ 0 ) + log c IG (› 0 1 , › 0 2 ) . (4.9)
By factorizing terms related to a in (4.9), the following Dirichlet distribution is obtained

q(a|K) = D(a|Ÿ)
where

'k oe K, Ÿk = Ÿ 0 k + ÿ joeJ E z Ë " k z j È .
By aggregating terms related to each Ê k in (4.9), a Normal distribution is obtained for each

Ê k such that q(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k | Ễk , ‡ 2 ˜ k 2
where 'k oe K,

˜ k = Q a ÿ joeJ E z Ë " k z j È (t j ) (t j ) T + ≠1 0 R b ≠1 , Ễk = ˜ k Q a ÿ joeJ E z Ë " k z j È x tj (t j ) + ≠1 0 Ê 0 R b .
Eventually, an Inverse Gamma distribution is deduced from (4.9) such that

q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 )
where

›1 = › 0 1 + J 2 , ›2 = › 0 2 + 1 2 ÿ koeK Q a ÿ joeJ E z Ë " k z j È x 2 tj + Ê T 0 ≠1 0 Ê 0 ≠ ỄT k ˜ ≠1 k Ễk R b .

Lower bound

Recalling that the lower bound on the log evidence is given by

L(q|K) = E z, [log p(x t , z, |t, K)] ≠ E z, [log q(z, |K)]
where E z, [log p(x t , z, |t, K)] is the free energy and E z, [log q(z, |K)] is the entropy of the approximate posterior q(z, |K). The free energy can be developed as

E z, [log p(x t , z, |t, K)] = E z, [log p(x t , z|t, , K)] + E [log p( |K)]
where

E z, [log p(x t , z|t, , K)] = ÿ joeJ ÿ koeK E z Ë " k z j È log fl jk and E [log p( |K)] = ÿ koeK (Ÿ 0 k ≠ 1)E [log a k ] ≠ 1 2 3 3(log 2fi + E Ë log ‡ 2 È ) + log | 0 | + E 5 1 ‡ 2 6 (E [Ê k ] ≠ Ê 0 ) T ≠1 0 (E [Ê k ] ≠ Ê 0 ) + Trace 1 ˜ k ≠1 0 2 4 ≠ (› 0 1 + 1)E Ë log ‡ 2 È ≠ › 0 2 E 5 1 ‡ 2 6 + log c D (Ÿ 0 ) + log c IG (› 0 1 , › 0 2 ) .
As for the entropy term, the following decompositon is obtained

E z, [log q(z, |K)] = E z [log q(z|K)] + E [log q( |K)]
where

E z [log q(z|K)] = ÿ joeJ ÿ koeK E z Ë " k z j È log rjk
and

E [log q( |K)] = ÿ koeK (Ÿ k ≠ 1)E [log a k ] ≠ 1 2 1 3(1 + log 2fi + E Ë log ‡ 2 È ) + log | ˜ k | 2 ≠ ( ›1 + 1)E Ë log ‡ 2 È ≠ ›2 E 5 1 ‡ 2 6 + log c D (Ÿ) + log c IG ( ›1 , ›2 ) .

Expectations

Expectations developed in variational calculations are derived from properties of variational posterior distributions and are obtained as follows. Categorical distribution properties lead to 'j oe J , 'k oe K :

E z Ë " k z j È = rjk .

Dirichlet distribution properties lead to

'k oe K :

E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b ,
where Â(•) is the digamma function. Normal distribution properties lead to

'k oe K : E [Ê k ] = Ễk , E Ë Ê k Ê T k È = V [Ê k ] + E [Ê k ] E [Ê k ] T = ‡ 2 ˜ k + Ễk ỄT k ,
Inverse Gamma distribution properties lead to

E 5 1 ‡ 2 6 = ›1 ›2 , E Ë log ‡ 2 È = log ›2 ≠ Â( ›1 ) .
Using all these properties, the following expectation can be calculated as 'j oe J , 'k oe K :

E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 D = ›1 1 x tj ≠ (t j ) T Ễk 2 2 ›2 + Trace 1 (t j ) T ˜ k (t j ) 2

Complete model

A model integrating parabolic data and mixed data is now presented. By taking into consideration any types of available data, the resulting model can fit data better and can estimate more accurate clusters. First, data formalism and assumptions are detailed. Then, the resulting mixture model and its inference procedure are developed.

Data and assumptions

In this part, data consist of J pulses gathering J amplitudes x t = (x tj ) joeJ associated to J times of arrival t = (t j ) joeJ , J continuous features x q = (x qj ) joeJ and J categorical features x c = (x cj ) joeJ from K distinct emitters. Let x j = (x qj , x cj , x tj ) the j th observation vector of mixed variables where

• x qj oe R d is a vector of d continuous radar features such as the Radio Frequency, the Pulse Width, the Azimuth or the Pulse Repetition Interval,

• x cj = ! x cj 0 , . . . , x cj q≠1 "
oe C q is a vector of q categorical radar modulations such as intrapulse modulations or pulse-to-pulse modulations,

• x tj oe R is a continuous variable modeling the Amplitude.

For each pulse j, the temporal evolution variable x tj and mixed variables (x qj , x cj ) are assumed to be independent conditionally to each cluster k oe K

'j oe J , (x q , x c )|z j = k |= x t |z j = k . ( 4
x cj = A x miss cj x obs cj B
with (x miss cj , x obs cj ) oe C q miss j ◊ C q obs j and q miss j + q obs j = q .

Mixture model

According to the independence assumption (4.10), the distribution of mixed data (Chapter 3) and the parabolic relation between temporal evolution data, the component distribution results in

'j oe J , p(x qj , x cj , x tj |z j = k) = p(x qj , x cj |z j = k)p(x tj |z j = k)
where

'j oe J , p(x qj , x cj |z j = k) = Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 22 " c x cj , p(x tj |t j , z j = k) = N 1 x tj | (t j ) T Ê k , ‡ 2 2 (4.11)
with

• u = (u j ) joeJ the scale latent variables handling outliers for quantitative data x q and distributed according to a Gamma distribution with shape and rate parameters (-, -) = (kc ,kc ) (k,c)oeK◊Cq conditionally to categorical data x c and labels z = (z j ) joeJ ,

• (µ, ) = ((µ kc ) coeCq , k ) koeK the mean and the variance parameters of quantitative data

x q for each cluster, • fi = (fi k ) koeK the weights of the multivariate Categorical distribution of categorical data x c for each cluster, • Ê = (Ê k ) koeK the regression parameters for temporal evolution data x t for each cluster,

• ‡ 2 the variance of the measurement noise related to temporal evolution data x t .

Recalling that p(z j = k) = a k where a = (a k ) koeK are the weights related to component distributions, the mixture model is obtained from (4.11) such that 'j oe J , fi,µ, ) is the set of parameters.

p(x j , u j |t j , ) = ÿ koeK a k N 1 x tj | (t j ) T Ê k , ‡ 2 2 Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c x cj (4.12) where = (a, Ê, ‡ 2 ,

Bayesian Framework

As in chapters 2 and 3, a Bayesian framework is used to estimate parameters . Assuming datasets (x = (x q , x c , x t ), t) of i.i.d observations (x j = (x qj , x cj , x tj ), t j ) joeJ , independent labels z = (z j ) joeJ and scale latent variables u = (u j ) joeJ , the complete likelihood associated to (4.12) is defined by The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x qj , x cj , x tj , t j ) decomposed into temporal evolution data (x tj , t j ) and mixed-type data (x qj , x cj ), the scale variables u j and the indicator variables z j . Known quantities, respectively unknown quantities, are in blue, respectively in red.

p(x, z, u|t, , K) = Ÿ joeJ Ÿ koeK Q a a k N 1 x tj | (t j ) T Ê k , ‡ 2 2 Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c x cj R b " k z j K J " 0 0 µ 0 ÷ 0 k µ k Ê 0 0 Ê k › 0 1 › 0 2 ‡ 2 fi 0 fi k -k -k u j x miss qj x obs qj x miss cj x obs cj z j a Ÿ 0 p 0 r 0 q 0 s 0 x tj (t j )
Eventually, the prior distribution required for is chosen as

p( |K) = p(a|K)p(Ê| ‡ 2 , K)p( ‡ 2 )

p(fi|K)p(-, -|K)p(µ, |K)

where

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) , p(fi|K) = Ÿ koeK D(fi k |fi 0 ) , p(µ, |K) = Ÿ koeK Ÿ coeCq N 1 µ kc |µ 0 , ÷ ≠1 0 k 2 IW( k |" 0 , 0 ) , p(-, -|K) = Ÿ koeK Ÿ coeCq p(-kc , -kc |p 0 , q 0 , s 0 , r 0 ) , p(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k |Ê 0 , ‡ 2 0 2 , p( ‡ 2 ) = IG( ‡ 2 |› 0 1 , › 0 2
) . Graphical representation of the proposed model is shown in Figure 4.3.

Inference

As previously, a factorized posterior distribution q(x miss q , u, x miss c , z, |K) = q(x miss q , u, x miss c , z|K)q( |K) is chosen as an approximation of the intractable posterior joint distribution p(x miss q , u, x miss c , z, |x, t, K) such that latent variables h = (x miss q , u, x miss c , z) and parameters are a posteriori independent and q(h|K) = q(x miss q |u, x miss c , z, K)q(u|x miss c , z, K)q(x miss c |z, K)q(z|K) , q( |K) = q(a|K)q(Ê| ‡ 2 , K)q( ‡ 2 )q(fi|K)q(-, -|K)q(µ, |K) .

According to VB assumptions, the following conjugate variational posterior distributions are obtained from the VB procedure

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ q(x miss q |u, x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC N 3 x miss qj | μx miss q jkc , u ≠1 j ˜ x miss q k 4 " c x cj " k z j , q(u|x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC G 1 u j |-jkc , -jkc 2 " c x cj " k z j , q(x miss c |z, K) = Ÿ joeJ Ÿ koeK MC(x miss cj |r x miss c jk ) " k z j q(z|K) = Ÿ joeJ Cat(z j |r j ) , q(a|K) = D(a|Ÿ) , q(fi|K) = Ÿ koeK D(fi| fk ) , q(µ, |K) = Ÿ koeK Ÿ coeC N 1 µ kc | μkc , ÷≠1 kc k 2 IW( k |" k , ˜ k ) , q(-, -|K) = Ÿ koeK p(-kc , -kc |p k , qk , sk , rk ) , q(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k | Ễk , ‡ 2 ˜ k 2 , q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 ) .
Their

E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b , log fl t jk = ≠ 1 2 A log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 DB , log fl qc jk = log ÿ c miss oeC q miss j fl x miss cj kc miss .
The red term E [log a k ] is deduced from properties of the Dirichlet distribution, the blue term log fl t jk is deduced from (4.7) and (4.8) in subsection 4.1.2 and the green term log fl qc jk has been detailed in Chapter 3. Hence, z is distributed a posteriori according to a product of Categorical distributions parametrized by r = (r jk ) (j,k)oeJ ◊K given by

'j oe J , 'k oe K, rjk = fl jk q koeK fl jk (4.16)
The lower bound on the log evidence is still required to master the VB inference and can be also decomposed into terms related to temporal evolution data (blue terms), mixed data (green terms) and labels z (red terms). This decomposition is obtained as follows

L(q|K) = E h, [log p(x, u, z, |t, K)] ≠ E h, Ë log q(x miss q , x miss c , u, z, |K) È
where the free energy can be developed as a,fi,µ,|K)] and the entropy as fi,µ

E h, [log p(x, u, z, |t, K)] = E h, [log p(x, u, z|t, , K)] + E Ë log p(Ê, ‡ 2 |K) È + E [log p(
E h, Ë log q(x miss q , x miss c , u, z, |K) È = E h, Ë log q(Ê, ‡ 2 |K) È + E h, Ë log q(x miss q , x miss c , u|z, K) È + E h, [log q(z|K)] + E h, [log q(a,
, |K)] .
Blue terms, respectively green terms, have been previously detailed in subsection 4.1.2 , respectively in Chapter 3. As for red terms, they are detailed below :

E h, [log p(x, u, z|t, , K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log fl jk , E h, [log q(z|K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log rjk .

Experiments

Two experiments are carried out to evaluate clustering performance with respect to a set of synthetic data and a set of real data. In the first experiment, only temporal evolution data are taken into consideration in the clustering procedure. Then, both temporal evolution data and quantitative data are considered in the second one. For comparison, the spectral clustering [START_REF] Von | A tutorial on spectral clustering[END_REF] and the k-means algorithm from [START_REF] John | Algorithm AS 136: A k-means clustering algorithm[END_REF] are also evaluated. First, characteristics of data, comparison algorithms and evaluation metrics are detailed. Then, both experiments are described and performance are shown to exhibit the e↵ectiveness of the proposed model.

Data, algorithms and metrics

Both synthetic and real data are composed of temporal evolution data related to amplitudes which are distributed according to a parabolic relation and quantitative data related to continuous radar features which are jointly distributed according to a multivariate normal distribution.

In synthetic data, temporal evolution data are generated by sampling a set of data from four parabolas directed by include unlabeled pulses that are mainly described by their Pulse Description World (PDW) composed of

Ê = Q c a ≠1 ≠2 ≠3 ≠4 1 2 3 4 1 2 3 4
• Time Of Arrival (TOA) ,

• Amplitude (A) ,

• Radio Frequency (RF) ,

• Pulse Width (PW) ,

• Azimuth (Az) .

Therefore, temporal evolution data are pairs (TOA,A) and quantitative data are triplets (RF,PW,Az). Unfortunately, these real data are classified and values of PDW cannot be released. Hence, axes of figures related to real data are not displayed. Three cases are obtained from real recordings and they are visible on Figure 4.6. Di↵erent numbers of parabolas and (RF,PW,Az) clusters can be observed according to a chosen real case. These di↵erences result in di↵erent numbers of emitters for the 3 cases such that 5, 4 and 2 emitters are identified in the cases 1, 2 and 3. Eventually, synthetic and real data are linearly transformed by a min-max normalization to meet algorithms requirements.

Except for the k-means algorithm, an initialisation is required for clustering algorithms that are involved in these experiments. The similarity graph required for the spectral clustering is obtained from a k-nearest neighbor graph as suggested in [START_REF] Von | A tutorial on spectral clustering[END_REF] where the number of neighbors k is chosen as the product of the log number of observations and the number of clusters. As for the proposed model, a supervised initialisation is retained due to its sensitivity to initialisation. First, prior hyperparameters › 0 1 and › 0 2 are initialised such that the prior mean

E[ 1 ‡ 2 ] = › 0 1 › 0 2
of the variance parameter ‡ 2 is equal to the inverse of the determinant of the covariance matrix of temporal evolution data points. This choice is motivated by the fact that the determinant of the covariance matrix can be interpreted as the generalized variance that reflects the overall spread of the data. Setting › 0 2 = 1, › 0 1 is initialised as the inverse of the generalized variance of the sample of temporal evolution data. Then, prior component means µ 0 , respectively covariance matrices 0 , are initialised from results of a k-means algorithm on quantitative data, respectively from diagonal matrices whose diagonal elements are variances of quantitative data. Other hyper-parameters are initialised as in Table 4.1. ÷ 0 " 0 p 0 r 0 q 0 s 0 (0, 0, 0) T I 3 0.5 100 1 1 1 1 1 the k-means algorithm (ARI = 0.33) can be explained by the fact that the k-means algorithm creates convex and isotropic clusters that cannot handle the parabolic structure of the generated data. This limitation is emphasized by higher Silhouette Coe cients of the k-means algorithm since the Silhouette Coe cient is generally higher for convex clusters. Moreover, the lower Silhouette Coe cients of the ground-truth for ‡ 2 oe {0.0001, 0.01} confirm the non-convexity of the data. Even if all algorithms poorly perform when data are embedded in noise ( ‡ 2 = 0.25), the proposed algorithm estimates clusters with a more parabolic shape than other algorithms which build more isotropic clusters (Subfigures (f), (i) and (l) in 4.3. Interpretation of algorithm performance through Silhouette coe cients and visual representations is complex since the Silhouette coe cient enhances algorithms that create convex clusters whereas the visual representations of estimated clusters tend to choose algorithms that create clusters with a parabolic shape. As in the example of real case 1, the proposed model succeeds in finding the radar emitter whose scanning behaviour is described by the red parabola (Subfigure (d) in Figure 4.8) whereas spectral clustering and k-means find that this parabola belongs to many emitters (Subfigures (g) and (j) in Figure 4.8). Nonetheless, spectral clustering and k-means provide higher Silhouette Coe cient (S = 0.46 and S = 0.57) than the proposed model (S = 0.05).

The second experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K when all types of data are taken into consideration. The number of clusters K is still set to 4 for synthetic data and to 5, 4 and 2 for the three real cases. Results of the second experiment on synthetic data are shown in Figure 4.9 and in Table 4.4. All algorithms succeed in clustering synthetic data for ‡ 2 oe {0.0001, 0.01, 0.25} since the ground-truth partition is recovered in Figure 4.9 with an ARI equals to 1 visible on Table 4.4. Adding quantitative information enables algorithms to recover the ground-truth for any value of ‡ 2 . Results of the second experiment on real data are shown in Figure 4.10 and in Table 4.5. The proposed model perfectly estimates clusters for the three cases whereas spectral clustering and k-means cannot manage to recover the correct clusters in real case 1 (Subfigures (d), (g) and (j)). Indeed, the proposed model finds the five di↵erent emitters by exploiting quantitative features while spectral clustering and k-means cannot identify the emitter whose temporal evolution features are distributed according to the parabola which is slowly increasing (red parabola in Subfigure (d)). Nonetheless, k-means and spectral clustering have higher Silhouette coe cients than the proposed model since they provide more convex clusters than the proposed method.

Piecewise parabolic data

In this section, K emitters presenting piecewise parabolic scanning behaviours are considered. Therefore, the main objective is to develop a mixture model which can build K distinct clusters formed by K piecewise parabolas with P piecewises and to define an inference procedure for parameter estimation. Then, the proposed model is enhanced with the mixture model designed for mixed data in Chapter 3 in order to improve clustering performance. Finally, experiments on synthetic data are carried out to exhibit performance of the proposed approach.

Model

Before introducing a mixture model that handles piecewise parabolic data, the piecewise parabolic relation between amplitudes and times of arrival is defined. Then, the mixture model is developed into a Bayesian framework.

Piecewise parabola equation

The piecewise parabolic relation between amplitudes (x tj ) joeJ and times of arrival (t j ) joeJ gives form to a set of P piecewises of constant amplitudes (µ t p ) poeP that are linked by a parabolic relation (Figure 4.11). Each piecewise µ t p gathers pulses (x tj , t j ) joeJp whose amplitude x tj is equal to µ t p and where J p is the set of indexes of pulses that belong to the p th piecewise. These P sets J p of pulses are disjoint and constitute a partition of J such t p J p = J . Finally, the piecewises (µ t p ) poeP belong to a parabola parametrized by Ê and ! min joeJp t j " poeP which are the times of the first pulses belonging to the piecewises. That definition can be translated into the following system

'j oe J p , Y _ ] _ [ x tj = µ t p + ' µ t p = 3 min joeJp t j 4 T Ê (4.17)
where ' ≥ N (0, ‡ 2 ) is a measurement noise introduced to model defects of materials, J p = {j oe J , |x tj ≠ µ t p | AE ‡} and (•), respectively Ê, are the polynomial transformation, respectively the regression parameter, defined in (4.2). Since the measurement error is assumed to be Gaussian, amplitude (x tj ) joeJp are distributed according to a normal distribution centered in µ t p with variance ‡ 2 'j oe

J p , x tj ≥ N 1 x tj |µ t p , ‡ 2 2 . ( 4.18)
If (J p ) poeP and (µ t p ) poeP are known, the regression parameter Ê is the solution of a linear problem given by µ t = (X ) T Ê where (X ) is a 3 ◊ P matrix whose columns are the P polynomial transformations ! ! min joeJp t j "" poeP and µ t = (µ t 1 , . . . , µ t P ) T . Then the regression parameter Ê is obtained such that

Ê = 1 (X ) (X ) T 2 ≠1 (X )µ t .

Mixture model

Since each radar emitter has its own scanning behaviour, K unique piecewise parabolas exist in data and they are configured with K regression parameters Ê = (Ê k ) koeK and K sets of piecewises

µ t = (µ t kp ) (p,k)oeP◊K .
Then, each amplitude x tj belongs to one of these sets of piecewises which is related to a specific emitter. In other words, conditionally to its label z j and its a liations to one of the piecewises, the amplitude x tj is distributed according to (4.18) such that the component distribution is defined by

'j oe J p , x tj |z j = k ≥ N 1 x tj |µ t kp , ‡ 2 2 (4.19)
In order to model its a liations to one of the piecewises, a latent discrete variable y j is introduced such that (4.19) becomes

'j oe J , x tj |y j = p, z j = k ≥ N 1 x tj |µ t kp , ‡ 2 2 (4.20)
where y j belongs to P and follows, conditionally to z j = k, a categorical distribution with weights b k = (b k1 , . . . , b kP ). Therefore the initial component distribution (4.19) can be reformulated as a mixture model such that

p(x tj |z j = k, , K) = ÿ poeP b kp N 1 x tj |µ t kp , ‡ 2 2
Recalling that p(z j = k) = a k where a = (a k ) koeK are the weights related to component distributions, the proposed mixture model is a mixture of mixture models given by

'j oe J , p(x tj | ) = ÿ koeK a k ÿ poeP b kp N 1 x tj |µ t kp , ‡ 2 2 (4.21) where = (a, b, µ t , ‡ 2
) is the set of parameters.

Bayesian framework

As in previous chapters, a Bayesian framework is used to estimate parameters . Assuming datasets (x t , t) of i.i.d observations (x tj , t j ) joeJ and independent labels z = (z j ) joeJ and y = (y j ) joeJ for clusters and piecewises, the complete likelihood associated to (4.21) is defined by

p(x t , h| , K) = Ÿ joeJ Ÿ koeK Q a a k Ÿ poeP 1 b kp N 1 x tj |µ t kp , ‡ 2 22 " p y j R b " k z j
where h = (y, z) is the set of latent variables. Eventually, the prior distribution required for is chosen as

p( |K) = p(a|K)p(b|K)p(µ t | ‡ 2 , K)p( ‡ 2 )
where a and b follow a Dirichlet distribution, each µ t kp follows a Normal distribution and ‡ 2 follows an Inverse Gamma distribution such that

Y _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) , p(b|K) = Ÿ koeK D(b k |o 0 ) , p(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |µ t 0 , • ≠1 0 ‡ 2 2 , p( ‡ 2 ) = IG( ‡ 2 |› 0 1 , › 0 2 ) .
The resulting mixture model is shown on Figure 4.12. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed amplitudes x j and the indicator variables (z j , y j ). Known quantities, respectively unknown quantities, are in blue, respectively in red.

K J › 0 2 › 0 1 ‡ 2 µ t 0 • 0 µ t k o 0 b k y j x tj z j a Ÿ 0

Inference

The Variational Bayes (VB) procedure is derived to estimate parameters of the mixture model defined in (4.21). Variational posterior distributions are obtained from the VB Expectation (VBE) and VB Maximization (VBM) steps and a lower bound on the log evidence is defined to master the convergence of the VB procedure.

Variational posterior distributions

As previously, a factorized posterior distribution q(h, |K) = q(h|K)q( |K) is chosen as an approximation of the intractable posterior joint distribution p(h, |x t , K) such that latent variables h and parameters are a posteriori independent and q(h|K) = q(y|z, K)q(z|K) , q( |K) = q(a|K)q(b|K)q(µ t | ‡ 2 K)q( ‡ 2 ) .

According to VB assumptions, the following conjugate variational posterior distributions are obtained from the VB procedure

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ q(y|z, K) = Ÿ joeJ Ÿ koeK Cat(y j |r y jk ) " k z j , q(z|K) = Ÿ joeJ Cat(z j |r j ) , q(a|K) = D(a|Ÿ) , q(b|K) = Ÿ koeK D(b k |õ k ) , q(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |μ t kp , • ≠1 kp ‡ 2 2 , q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 ) .
Their respective parameters are estimated during the VBE and VBM steps.

VBE-step

The VBE-step consists in deriving the following expectation

E [log p(x t , h|t, , K)] = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + ÿ poeP " p y j 3 E [log b kp ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ µ t kp ) 2 ‡ 2 D 444 = ÿ joeJ ÿ koeK " k z j Q a E [log a k ] + ÿ poeP " p y j log fl y jkp R b (4.22)
where

log fl y jkp = E [log b kp ] ≠ 1 2 A log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ µ t kp ) 2 ‡ 2 DB . ( 4.23)
Hence, a categorical distribution for piecewise labels y is deduced from (4.22) conditionally to cluster labels z such that q(y|z, K) =

VBM-step

The VBM-step consists in deriving the following expectation

E h [log p(x t , h, |t, K)] = E h [log p(x t , h|t, , K)] + log p( |K) = ÿ joeJ ÿ koeK E h Ë " k z j È 3 log a k + ÿ poeP E h Ë " p y j È 3 E [log b kp ] ≠ 1 2 3 log 2fi + log ‡ 2 + (x tj ≠ µ t kp ) 2 ‡ 2 444 + ÿ koeK (Ÿ 0 k ≠ 1) log a k + log c D (Ÿ 0 ) ≠ (› 0 1 + 1) log ‡ 2 ≠ › 0 2 ‡ 2 + log c IG (› 0 1 , › 0 2 ) + ÿ koeK ÿ poeP (o 0 kp ≠ 1) log b kp ≠ 1 2 3 log 2fi + log ‡ 2 + • 0 ‡ 2 (µ t kp ≠ µ t 0 ) 2 4 + ÿ koeK log c D (o 0 k ) . (4.24)
By factorizing terms related to a in (4.24), the following Dirichlet distribution is obtained

q(a|K) = D(a|Ÿ)
where

'k oe K, Ÿk = Ÿ 0 k + ÿ joeJ E h Ë " k z j È .
Following the same reasoning, b is distributed according to a product of Dirichlet distributions given by q(b|K)

= Ÿ koeK D(b k |õ k )
where 'k oe K, 'p oe P, õkp = o 0 kp +

ÿ joeJ E h Ë " k z j È E h Ë " p y j È .
By aggregating terms related to each µ t kp in (4.24), a Normal distribution is obtained for each

µ t kp such that q(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |μ t kp , • ≠1 kp ‡ 2 2
where 'k oe K and 'p oe P

•kp = ÿ joeJ E h Ë " k z j È E h Ë " p y j È + • 0 , μt kp = q joeJ E h Ë " k z j È E h Ë " p y j È x tj + • 0 µ t 0 •kp .
Eventually, an Inverse Gamma distribution is deduced from (4.24) such that

q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 )
where

›1 = › 0 1 + J 2 , ›2 = › 0 2 + 1 2 ÿ koeK ÿ poeP Q a ÿ joeJ E h Ë " k z j È E h Ë " p y j È x 2 tj + • 0 µ t 0 2 ≠ •kp (μ t kp ) 2 R b .

Lower bound

Recalling that the lower bound on the log evidence is given by

L(q|K) = E h, [log p(x t , h, |t, K)] ≠ E h, [log q(h, |K)]
where E h, [log p(x t , h, |t, K)] is the free energy and E h, [log q(h, |K)] is the entropy of the approximate posterior q(h, |K). The free energy can be developed as

E h, [log p(x t , h, |t, K)] = E h, [log p(x t , h|t, , K)] + E [log p( |K)]
where

E h, [log p(x t , h|t, , K)] = ÿ joeJ ÿ koeK E h Ë " k z j È Q a E [log a k ] + ÿ poeP E h Ë " p y j È log fl y jkp R b
and

E [log p( |K)] = ÿ koeK (Ÿ 0 k ≠ 1)E [log a k ] ≠ (› 0 1 + 1)E Ë log ‡ 2 È ≠ › 0 2 E 5 1 ‡ 2 6 + log c D (Ÿ 0 ) + log c IG (› 0 1 , › 0 2 ) + ÿ koeK ÿ poeP (o 0 kp ≠ 1)E [log b kp ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E 5 1 ‡ 2 6 • 0 (E Ë µ t kp È ≠ µ t 0 ) 2 + • 0 • ≠1 kp 4 + ÿ koeK log c D (o 0 k ) .
As for the entropy term, the following decompositon is obtained

E h, [log q(h, |K)] = E h [log q(y, z|K)] + E [log q( |K)] = E h [log q(y|z, K)] + E h [log q(z|K)] + E [log q( |K)]
where

E h [log q(y|z, K)] = ÿ joeJ ÿ koeK E h Ë " k z j È ÿ poeP E h Ë " p y j È log ry jkp , E h [log q(z|K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log rjk .
and

E [log q( |K)] = ÿ koeK (Ÿ k ≠ 1)E [log a k ] ≠ ( ›1 + 1)E Ë log ‡ 2 È ≠ ›2 E 5 1 ‡ 2 6 + log c D (Ÿ) + log c IG ( ›1 , ›2 ) + ÿ koeK ÿ poeP (õ kp ≠ 1)E [log b kp ] ≠ 1 2 1 log 2fi + E Ë log ‡ 2 È + 1 2 + ÿ koeK log c D (õ k ) .

Expectations

Expectations developed in variational calculations are derived from properties of variational posterior distributions and are obtained as follows. Categorical distribution properties lead to 'j oe J , 'k oe K, 'p oe P :

E h Ë " p y j È = ry jkp , E h Ë " k z j È = rjk .

Dirichlet distribution properties lead to

'k oe K, 'p oe P :

E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b , E [log b kp ] = Â(õ kp ) ≠ Â Q a ÿ poeP õkp R b .
where Â(•) is the digamma function. Normal distribution properties lead to 'k oe K, 'p oe P :

E Ë µ t kp È = μt kp , E Ë (µ t kp ) 2 È = V Ë µ t kp È + E Ë µ t kp È 2 = ‡ 2 • ≠1 kp + (μ t kp ) 2 ,
Inverse Gamma distribution properties lead to

E 5 1 ‡ 2 6 = ›1 ›2 , E Ë log ‡ 2 È = log ›2 ≠ Â( ›1 ) .
Using all these properties, the following expectation can be calculated as 'j oe J , 'k oe K, 'p oe P :

E C (x tj ≠ µ t kp ) 2 ‡ 2 D = ›1 1 x tj ≠ μt kp 2 2 ›2 + • ≠1 kp .

Complete model

A model integrating piecewise parabolic data and mixed data is now presented. By taking into consideration any types of available data, the resulting model can fit data better and can estimate more accurate clusters. First, data formalism and assumptions are detailed. Then, the resulting mixture model and its inference procedure are developed.

Data and assumptions

In this part, data consist of J pulses gathering J amplitudes x t = (x tj ) joeJ associated to J times of arrival t = (t j ) joeJ , J continuous features x q = (x qj ) joeJ and J categorical features x c = (x cj ) joeJ from K distinct emitters. Let x j = (x qj , x cj , x tj ) the j th observation vector of mixed variables where

• x qj oe R d is a vector of d continuous radar features such as the Radio Frequency, the Pulse Width, the Azimuth or the Pulse Repetition Interval,

• x cj = ! x cj 0 , . . . , x cj q≠1 "
oe C q is a vector of q categorical radar modulations such as intrapulse modulations or pulse-to-pulse modulations,

• x tj oe R is a continuous variable modeling the Amplitude.

For each pulse j, the temporal evolution variable x tj and mixed variables (x qj , x cj ) are assumed to be independent conditionally to each cluster k oe K 'j oe J , (x q , x c )|z j = k |= x t |z j = k . (4.25) with z j the latent variable modeling the label of the j th observation vector x j = (x qj , x cj , x tj ).

Moreover, the temporal evolution data (x tj , t j ) joeJ are distributed according to a piecewise parabolic relation and the quantitative data with (x miss cj , x obs cj ) oe C q miss j ◊ C q obs j and q miss j + q obs j = q .

(

Mixture model

According to the independence assumption (4.25), the distribution of mixed data (Chapter 3) and the piecewise parabolic relation between temporal evolution data, the component distribution results in

'j oe J , p(x qj , x cj , x tj |u j , z j = k) = p(x qj , x cj |u j , z j = k)p(x tj |z j = k)
where

'j oe J , p(x qj , x cj |u j , z j = k) = Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 22 " c x cj , p(x tj |z j = k) = Ÿ poeP 1 b kp N 1 x tj |µ t kp , ‡ 2
22 " p y j (4.26) with

• u = (u j ) joeJ the scale latent variables handling outliers for quantitative data x q and distributed according to a Gamma distribution with shape and rate parameters (-, -) = (kc ,kc ) (k,c)oeK◊Cq conditionally to categorical data x c and labels z = (z j ) joeJ ,

• (µ, ) = ((µ kc ) coeCq , k ) koeK the mean and the variance parameters of quantitative data

x q for each cluster,

• fi = (fi k ) koeK the weights of the multivariate Categorical distribution of categorical data

x c for each cluster,

• y = (y j ) joeJ the latent variables indicating the p th piecewise temporal evolution data x t belong to,

• b = ((b kp ) poeP ) koeK the weights of the Categorical distribution of latent variables y,

• µ t = ((µ t kp ) poeP ) koeK the set of piecewises for temporal evolution data x t for each cluster,

• ‡ 2 the variance of the measurement noise related to temporal evolution data x t . The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x qj , x cj , x tj , t j ) decomposed into temporal evolution data (x tj , t j ) and mixed-type data (x qj , x cj ), the scale variables u j and the indicator variables (y j , z j ). Known quantities, respectively unknown quantities, are in blue, respectively in red.

K J " 0 0 µ 0 ÷ 0 k µ k µ t 0 • 0 µ t k › 0 1 › 0 2 ‡ 2 fi 0 fi k -k -k u j x miss qj x obs qj x miss cj x obs cj z j y j b o 0 a Ÿ 0 p 0 r 0 q 0 s 0 x tj (t j )
Recalling that p(z j = k) = a k where a = (a k ) koeK are the weights related to component distributions, the mixture model is obtained from (4.26) such that 'j oe J , p(x j , u j , y j | ) =

Bayesian framework

As in previous chapters 2 and 3, a Bayesian framework is used to estimate parameters . Assuming datasets (x = (x q , x c , x t ), t) of i.i.d observations (x j = (x qj , x cj , x tj ), t j ) joeJ , independent labels z = (z j ) joeJ , piecewise indicators y = (y j ) joeJ and scale latent variables u = (u j ) joeJ , the complete likelihood associated to (4.27) is defined by

p(x, y, u, z| , K) = Ÿ joeJ Ÿ koeK 3 a k Ÿ poeP 1 b kp N 1 x tj |µ t kp , ‡ 2 22 " p y j ◊ Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c x cj 4 " k z j
Eventually, the prior distribution required for is chosen as

p( |K) = p(a|K)p(b|K)p(µ t | ‡ 2 , K)p( ‡ 2 )p(fi|K)p(-, -|K)p(µ, |K)
where

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) , p(b|K) = Ÿ koeK D(b k |o 0 ) , p(fi|K) = Ÿ koeK D(fi k |fi 0 ) , p(µ, |K) = Ÿ koeK Ÿ coeCq N 1 µ kc |µ 0 , ÷ ≠1 0 k 2 IW( k |" 0 , 0 ) , p(-, -|K) = Ÿ koeK Ÿ coeCq p(-kc , -kc |p 0 , q 0 , s 0 , r 0 ) , p(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |µ t 0 , • ≠1 0 ‡ 2 2 , p( ‡ 2 ) = IG( ‡ 2 |› 0 1 , › 0 2 ) .
Graphical representation of the proposed model is shown in Figure 4.13.

As previously, a factorized posterior distribution q(x miss q , u, x miss c , y, z, |K) = q(x miss q , u, x miss c , y, z|K)q( |K) is chosen as an approximation of the intractable posterior joint distribution p(x miss q , u, x miss c , y, z, |x obs , K) such that latent variables h = (x miss q , u, x miss c , y, z) and parameters are a posteriori independent and q(h|K) = q(x miss q |u, x miss c , z, K)q(u|x miss c , z, K)q(x miss c |z, K)q(y|z, K)q(z|K) , q( |K) = q(a|K)q(b|K)q(µ t | ‡ 2 , K)q( ‡ 2 )q(fi|K)q(-, -|K)q(µ, |K) .

According to VB assumptions, the following conjugate variational posterior distributions are obtained from the VB procedure The lower bound on the log evidence is still required to master the VB inference and can be also decomposed into terms related to temporal evolution data (blue terms), mixed data (green terms) and labels z (red terms). This decomposition is obtained as follows

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ q(x miss q |u, x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC N 3 x miss qj | μx miss q jkc , u ≠1 j ˜ x miss q k 4 " c x cj " k z j , q(u|x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC G 1 u j |-jkc , -jkc 2 " c x cj " k z j , q(x miss c |z, K) = Ÿ joeJ Ÿ koeK MC(x miss cj |r x miss c jk ) " k z j q(y|z, K) = Ÿ joeJ Ÿ koeK Cat(y j |r y jk ) " k z j , q(z|K) = Ÿ joeJ Cat(z j |r j ) , q(a|K) = D(a|Ÿ) , q(fi|K) = Ÿ koeK D(fi| fk ) , q(µ, |K) = Ÿ koeK Ÿ coeC N 1 µ kc | μkc , ÷≠1 kc k 2 IW( k |" k , ˜ k ) , q(-, -|K) = Ÿ koeK p(-kc , -kc |p k , qk , sk , rk ) , q(b|K) = Ÿ koeK D(b k |õ k ) , q(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |μ t kp , • ≠1 kp ‡ 2 2 , q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 ) . ( 4
E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b , log fl t jk = log
L(q|K) = E h, [log p(x, u, y, z, |K)] ≠ E h, Ë log q(x miss q , x miss c , u, y, z, |K) È
where the free energy can be developed as fi,µ,|K)] and the entropy as

E h, [log p(x, u, y, z, |K)] = E h, [log p(x, u, y, z| , K)] + E Ë log p(b, µ t , ‡ 2 |K) È + E [log p(a,
E h, Ë log q(x miss q , x miss c , u, y, z, |K) È = E h, Ë log q(b, µ t , ‡ 2 |K) È + E h, [log q(y|z, K)] + E h, Ë log q(x miss q , x miss c , u|z, K) È + E h, [log q(z|K)] + E h, [log q(a, fi, -, -, µ, |K)] .
Blue terms, respectively green terms, have been previously detailed in subsection 4.2.2 , respectively in Chapter 3. As for red terms, they are detailed below :

E h, [log p(x, u, y, z|t, , K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log fl jk , E h, [log q(z|K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log rjk .

Experiments

Two experiments are carried out to evaluate clustering performance with respect to a set of synthetic data. In the first experiment, only temporal evolution data are taken into consideration in the clustering procedure. Then, both temporal evolution data and quantitative data are considered in the second one. For comparison, the spectral clustering [START_REF] Von | A tutorial on spectral clustering[END_REF] and the k-means algorithm from [START_REF] John | Algorithm AS 136: A k-means clustering algorithm[END_REF] are also evaluated. First, characteristics of data, comparison algorithms and evaluation metrics are detailed. Then, both experiments are described and performance are shown to exhibit the e↵ectiveness of the proposed model.

Data, algorithms and metrics

Synthetic data are composed of temporal evolution data related to amplitudes which are distributed according to a piecewise parabolic relation and quantitative data related to continuous radar features which are jointly distributed according to a multivariate normal distribution. Temporal evolution data are generated by sampling a set of data from four piecewise parabolas directed by

Ê = Q c a ≠1 ≠2 ≠3 ≠4 1 2 3 4 1 2 3 4 R d b .
For each piecewise parabola, p = 4 piecewises are obtained by dividing the time interval in p equal subintervals and assigning to each piecewise the value of the parabola at the minimal time of its related time subinterval. Quantitative data are generated by sampling a set of data from four well-separated bivariate clusters with centers [0, 0] T , [1, 0] T ,[0, 1] T and [1, 1] T and identity covariance matrices. Three synthetic datasets are generated with respect to a range of values of ‡ 2 and are linearly transformed by a min-max normalization to meet algorithms requirements. These datasets are shown in Figures 4.14 Except for the k-means algorithm, an initialisation is required for clustering algorithms that are involved in these experiments. The similarity graph required for the spectral clustering is obtained from a k-nearest neighbor graph as suggested in [START_REF] Von | A tutorial on spectral clustering[END_REF] where the number of neighbors k is chosen as the product of the log number of observations and the number of clusters. As for the proposed model, a supervised initialisation is retained due to its sensitivity to initialisation. First, prior hyperparameters › 0 1 and › 0 2 are initialised such that the prior mean

E[ 1 ‡ 2 ] = › 0 1 › 0 2
of the variance parameter ‡ 2 is equal to the inverse of the determinant of the covariance matrix of temporal evolution data points. This choice is motivated by the fact that the determinant of the covariance matrix can be interpreted as the generalized variance that reflects the overall spread of the data. Setting › 0 2 = 1, › 0 1 is initialised as the inverse of the generalized variance of the sample of temporal evolution data. In addition, prior piecewise means µ t 0 are initialised from results of a k-means algorithm on temporal evolution data. Then, prior component means • 0 Ÿ 0 ÷ 0 " 0 p 0 r 0 q 0 s 0 1 0.5 100 1 1 1 1 1 Performance on synthetic data are evaluated through the Adjusted Rand Index (ARI) [HA85] that compares estimated partitions of data with the ground-truth and the Silhouette Coe cient [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF] which does not require the ground-truth and provides a higher score when clusters are dense and well separated.

Experiments and results

The first experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K when only temporal evolution data are taken into consideration. According to datasets visualised in Figure 4.14, K is set to 4 for synthetic data. Results of the first experiment on synthetic data are shown in Figure 4.16 and in Table 4.7. When ‡ 2 oe {0.0001, 0.01}, the proposed model and the spectral clustering have similar performance in clustering synthetic data with an ARI equals to 0.73 while the k-means algorithm has the lowest performance (ARI = 0.35) in creating convex and isotropic clusters that cannot handle the piecewise parabolic structure of the generated data. This limitation is emphasized by higher Silhouette Coe cients of the k-means algorithm whereas the non-convexity of the data is confirmed by the lower Silhouette Coe cients of the ground-truth. Even if all algorithms poorly perform when data are embedded in noise ( ‡ 2 = 0.25), the proposed algorithm estimates clusters with a more parabolic shape than other algorithms which build more isotropic clusters (Subfigures (f), (i) and (l) in The second experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K when all types of data are taken into consideration. The number of clusters K is still set to 4 for synthetic data.. Results of the second experiment on synthetic data are shown in Figure 4.17 and in Table 4.8. All algorithms succeed in clustering synthetic data for ‡ 2 oe {0.0001, 0.01, 0.25} since the ground-truth partition is recovered in Figure 4.17 with an ARI equals to 1 visible on Table 4.8. Adding quantitative information enables algorithms to recover the ground-truth for any value of ‡ 2 .

Parabolic and piecewise parabolic data

In this section, we consider that both scanning behaviours are observed among the K emitters. Hence, parabolic and piecewise parabolic relations are observed in data and have to be taken into account in the clustering procedure by developing a mixture model that can build K distinct clusters formed by either parabolas or piecewise parabolas. Then, the proposed model is enhanced with the mixture model designed for mixed data in Chapter 3 in order to improve clustering performance. Finally, experiments on synthetic data are carried out to exhibit performance of the proposed approach.

Model

Definitions of parabolic relation and piecewise parabolic have been previously introduced in Sections 4.1 and 4.2. Now, data can be modeled either by a parabolic relation (4.1) or a piecewise parabolic relation (4.17). Since the measurement noise ' is still Gaussian, amplitudes (x tj ) joeJ are normally distributed according to (4.3) and (4.18) such that

'j oe J , x tj |t j ≥ Y _ ] _ [ N 1 x tj | (t j ) T Ê, ‡ 2 2 if (x

Mixture model

Since each radar emitter has its own scanning behaviour, radar emitters can be distinguished into two groups where K 0 unique parabolas exist in the first group and K 1 piecewise parabolas exist in the second one such that K = K 0 + K 1 . Therefore K 0 from K regression parameters Ê = (Ê k ) koeK and K 1 from K sets of piecewises µ t = (µ t kp ) (p,k)oeP◊K have to be estimated. Then, each amplitude x tj belongs to one of these sets which is related to a specific emitter. In other words, conditionally to its label z j , the amplitude x tj is distributed according to (4.33) such that the component distribution is defined by

'j oe J , x tj |t j , z j = k ≥ Y _ ] _ [ N 1 x tj | (t j ) T Ê k , ‡ 2 2 if k oe K 0 . ÷p oe P, N 1 x tj |µ t kp , ‡ 2 2 if k oe K 1 . (4.34)
In order to model its a liations to one of the groups, a latent discrete variable w j is introduced such that (4.34) becomes

'j oe J , x tj |t j , w j = 0, z j = k ≥ N 1 x tj | (t j ) T Ê k , ‡ 2 2
x tj |y j = p, w j = 1,

z j = k ≥ N 1 x tj |µ t kp , ‡ 2 2 (4.35)
where w j oe {0, 1} follows, conditionally to z j = k, a categorical distribution with weights c k = (c k0 , c k1 ) and y j is the latent variable defined in (4.20) that follows, conditionally to z j = k and w j = 1, a categorical distribution with weights b k = (b k1 , . . . , b kP ). Therefore the initial component distribution (4.34) can be reformulated as

p(x tj |t j , z j = k, , K) = c k0 N 1 x tj | (t j ) T Ê k , ‡ 2 2 + c k1 ÿ poeP b kp N 1 x tj |µ t kp , ‡ 2 2
Recalling that p(z j = k) = a k where a = (a k ) koeK are the weights related to component distributions, the proposed mixture model is a mixture of mixture models given by 'j oe J , p(x tj |t j , ) =

ÿ koeK a k Q a c k0 N 1 x tj | (t j ) T Ê k , ‡ 2 2 + c k1 ÿ poeP b kp N 1 x tj |µ t kp , ‡ 2 2 R b (4.36)
where = (a, b, c, Ê, µ t , ‡ 2 ) is the set of parameters.

Bayesian framework

As in previous chapters, a Bayesian framework is used to estimate parameters . Assuming datasets (x t , t) of i.i.d observations (x tj , t j ) joeJ and independent labels z = (z j ) joeJ , w = (w j ) joeJ and y = (y j ) joeJ for clusters, scanning types and piecewises, the complete likelihood associated to (4.36) is defined by

p(x t , h| , K) = Ÿ joeJ Ÿ koeK Q c aa k 1 c k0 N 1 x tj | (t j ) T Ê k , ‡ 2 22 " 0 w j Q a c k1 Ÿ poeP 1 b kp N 1 x tj |µ t kp , ‡ 2 22 " p y j R b " 1 w j R d b " k z j
.

where h = (y, w, z) is the set of latent variables. Eventually, the prior distribution required for is chosen as The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x tj , t j ) decomposed into the amplitude x j and the polynomial transformation (t j ) and the indicator variables (y j , w j , z j ). Known quantities, respectively unknown quantities, are in blue, respectively in red.

p( |K) = p(a|K)p(b|K)p(c|K)p(Ê| ‡ 2 , K)p(µ t | ‡ 2 , K)p( ‡ 2 ) K J › 0 2 › 0 1 ‡ 2 µ t 0 • 0 µ t k ' 0 c k o 0 b k Ê 0 0 Ê k y j x j (t j ) w j z j a Ÿ 0
where a, b k and c k follow a Dirichlet distribution, each Ê k and µ t kp follow a Normal distribution and ‡ 2 follows an Inverse Gamma distribution such that

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) , p(b|K) = Ÿ koeK D(b k |o 0 ) , p(c|K) = Ÿ koeK D(c k |' 0 ) , p(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k |Ê 0 , ‡ 2 0 2 , p(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |µ t 0 , • ≠1 0 ‡ 2 2 , p( ‡ 2 ) = IG( ‡ 2 |› 0 1 , › 0 2 ) .
The resulting mixture model is shown on Figure 4.19.

Inference

The Variational Bayes (VB) procedure is derived to estimate parameters of the mixture model defined in (4.36). Variational posterior distributions are obtained from the VB Expectation (VBE) and VB Maximization (VBM) steps and a lower bound on the log evidence is defined to master the convergence of the VB procedure.

Variational posterior distributions

As previously, a factorized posterior distribution q(h, |K) = q(h|K)q( |K) is chosen as an approximation of the intractable posterior joint distribution p(h, |x t , K) such that latent variables h and parameters are a posteriori independent and q(h|K) = q(y|w, z, K)q(w|z, K)q(z|K) , q( |K) = q(a|K)q(b|K)q(c|K)q(Ê| ‡ 2 , K)q(µ t | ‡ 2 , K)q( ‡ 2 ) .

According to VB assumptions, the following conjugate variational posterior distributions are obtained from the VB procedure

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ q(y|w, z, K) = Ÿ joeJ Ÿ koeK Cat(y j |r y jk ) " 1 w j " k z j , q(w|z, K) = Ÿ joeJ Ÿ koeK Cat(y j |r w jk ) " k z j , q(z|K) = Ÿ joeJ Cat(z j |r j ) , q(a|K) = D(a|Ÿ) , q(b|K) = Ÿ koeK D(b k |õ k ) , q(c|K) = Ÿ koeK D(c k | 'k ) , q(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k | Ễk , ‡ 2 ˜ k 2 , q(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |μ t kp , • ≠1 kp ‡ 2 2 , q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 ) .
Their respective parameters are estimated during the VBE and VBM steps.

VBE-step

The VBE-step consists in deriving the following expectation where where where

E [log p(x t , h|t, , K)] = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + " 0 w j 3 E [log c k0 ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 D 44 + " 1 w j 3 E [log c k1 ] + ÿ poeP " p y j 3 E [log b kp ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ µ t kp ) 2 ‡ 2 D 4444 = ÿ joeJ ÿ koeK " k z j 3 E [log a k ] + " 0 w j 3 E [log c k0 ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 D 44 + " 1 w j 3 E [log c k1 ] +
log fl y jkp = E [log b kp ] ≠ 1 2 A log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ µ t kp ) 2 ‡ 2 DB . ( 4
log fl w jk0 = E [log c k0 ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2
fl jk = E [log a k ] + log ÿ ioe{0,1}
fl w jki .

VBM-step

The VBM-step consists in deriving the following expectation

E h [log p(x t , h, |t, K)] = E h [log p(x t , h|t, , K)] + log p( |K) = ÿ joeJ ÿ koeK E h Ë " k z j È 3 log a k + E h Ë " 0 w j È 3 log c k0 ≠ 1 2 3 log 2fi + log ‡ 2 + (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 44 + E h Ë " 1 w j È 3 log c k1 + ÿ poeP E h Ë " p y j È 3 E [log b kp ] ≠ 1 2 3 log 2fi + log ‡ 2 + (x tj ≠ µ t kp ) 2 ‡ 2 4444 + ÿ koeK (Ÿ 0 k ≠ 1) log a k + log c D (Ÿ 0 ) ≠ 1 2 A 3(log 2fi + log ‡ 2 ) + log | 0 | + (Ê k ≠ Ê 0 ) T ≠1 0 ‡ 2 (Ê k ≠ Ê 0 ) B ≠ (› 0 1 + 1) log ‡ 2 ≠ › 0 2 ‡ 2 + log c IG (› 0 1 , › 0 2 ) + ÿ koeK ÿ poeP (o 0 kp ≠ 1) log b kp ≠ 1 2 3 log 2fi + log ‡ 2 + • 0 ‡ 2 (µ t kp ≠ µ t 0 ) 2 4 + ÿ koeK ÿ ioe{0,1} (' 0 ki ≠ 1) log c ki + ÿ koeK log c D (o 0 k ) + log c D (' 0 k ) . ( 4.41)
By factorizing terms related to a in (4.41), the following Dirichlet distribution is obtained

q(a|K) = D(a|Ÿ) where 'k oe K, Ÿk = Ÿ 0 k + ÿ joeJ E h Ë " k z j È .
Following the same reasoning, b and c are distributed according to a product of Dirichlet distributions given by

q(b|K) = Ÿ koeK D(b k |õ k ) , q(c|K) = Ÿ koeK D(c k | 'k ) , where 'k oe K, 'p oe P, õkp = o 0 kp + ÿ joeJ E h Ë " k z j È E h Ë " 1 w j È E h Ë " p y j È and 'k oe K, 'i oe {0, 1}, 'ki = ' 0 ki + ÿ joeJ E h Ë " k z j È E h Ë " i w j È .
By aggregating terms related to each µ t kp and Ê k in (4.41), a Normal distribution is obtained for each µ t kp and Ê k such that

q(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |μ t kp , • ≠1 kp ‡ 2 2 , q(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k | Ễk , ‡ 2 ˜ k 2 ,
where 'k oe K and 'p oe P

•kp = ÿ joeJ E h Ë " k z j È E h Ë " 1 w j È E h Ë " p y j È + • 0 , μt kp = q joeJ E h Ë " k z j È E h Ë " 1 w j È E h Ë " p y j È x tj + • 0 µ t 0 •kp , ˜ k = Q a ÿ joeJ E h Ë " k z j È E h Ë " 0 w j È (t j ) (t j ) T + ≠1 0 R b ≠1 , Ễk = ˜ k Q a ÿ joeJ E h Ë " k z j È E h Ë " 0 w j È x tj (t j ) + ≠1 0 Ê 0 R b .
Eventually, an Inverse Gamma distribution is deduced from (4.41) such that

q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 )
where

›1 = › 0 1 + J 2 , ›2 = › 0 2 + 1 2 ÿ koeK ÿ joeJ 3 ÿ poeP E h Ë " k z j È E h Ë " 1 w j È E h Ë " p y j È x 2 tj + • 0 (µ t 0 ) 2 ≠ •kp (μ t kp ) 2 + E h Ë " k z j È E h Ë " 0 w j È x 2 tj + Ê T 0 ≠1 0 Ê 0 ≠ ỄT k ˜ ≠1 k Ễk 4 .

Lower bound

Recalling that the lower bound on the log evidence is given by

L(q|K) = E h, [log p(x t , h, |t, K)] ≠ E h, [log q(h, |K)]
where E h, [log p(x t , h, |t, K)] is the free energy and E h, [log q(h, |K)] is the entropy of the approximate posterior q(h, |K). The free energy can be developed as

E h, [log p(x t , h, |t, K)] = E h, [log p(x t , h|t, , K)] + E [log p( |K)]
where

E h, [log p(x t , h|t, , K)] = ÿ joeJ ÿ koeK E h Ë " k z j È 3 E [log a k ] + E h Ë " 0 w j È log fl w jk0 + E h Ë " 1 w j È 3 E [log c k1 ] + ÿ poeP E h Ë " p y j È log fl y jkp 44 and E [log p( |K)] = ÿ koeK (Ÿ 0 k ≠ 1)E [log a k ] ≠ (› 0 1 + 1)E Ë log ‡ 2 È ≠ › 0 2 E 5 1 ‡ 2 6 + log c D (Ÿ 0 ) + log c IG (› 0 1 , › 0 2 ) ≠ 1 2 3 3(log 2fi + E Ë log ‡ 2 È ) + log | 0 | + Trace 1 ˜ k ≠1 0 2 + E 5 1 ‡ 2 6 (E [Ê k ] ≠ Ê 0 ) T ≠1 0 (E [Ê k ] ≠ Ê 0 ) 4 + ÿ koeK ÿ poeP (o 0 kp ≠ 1)E [log b kp ] ≠ 1 2 3 log 2fi + E Ë log ‡ 2 È + E 5 1 ‡ 2 6 • 0 (E Ë µ t kp È ≠ µ t 0 ) 2 + • 0 • ≠1 kp 4 + ÿ koeK ÿ ioe{0,1} (' 0 ki ≠ 1)E [log c ki ] + ÿ koeK log c D (o 0 k ) + log c D (' 0 k ) .
As for the entropy term, the following decompositon is obtained

E h, [log q(h, |K)] = E h [log q(y, w, z|K)] + E [log q( |K)] = E h [log q(y|w, z, K)] + E h [log q(w|z, K)] + E h [log q(z|K)] + E [log q( |K)]
where

E h [log q(y|w, z, K)] = ÿ joeJ ÿ koeK E h Ë " k z j È E h Ë " 1 w j È ÿ poeP E h Ë " p y j È log ry jkp , E h [log q(w|z, K)] = ÿ joeJ ÿ koeK E h Ë " k z j È ÿ ioe{0,1} E h Ë " i w j È log rw jki , E h [log q(z|K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log rjk and E [log q( |K)] = ÿ koeK (Ÿ k ≠ 1)E [log a k ] ≠ ( ›1 + 1)E Ë log ‡ 2 È ≠ ›2 E 5 1 ‡ 2 6 + log c D (Ÿ) + log c IG ( ›1 , ›2 ) ≠ 1 2 1 3(1 + log 2fi + E Ë log ‡ 2 È ) + log | ˜ k | 2 + ÿ koeK ÿ poeP (õ kp ≠ 1)E [log b kp ] ≠ 1 2 1 log 2fi + E Ë log ‡ 2 È + 1 2 + ÿ koeK log c D (õ k ) + c D ( 'k ) .

Expectations

Expectations developed in variational calculations are derived from properties of variational posterior distributions and are obtained as follows. Categorical distribution properties lead to 'j oe J , 'k oe K, 'p oe P, 'i oe {0, 1} :

E h Ë " p y j È = ry jkp , E h Ë " i w j È = rw jki , E h Ë " k z j È = rjk .
Dirichlet distribution properties lead to 'k oe K, 'p oe P, 'i oe {0, 1} :

E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b , E [log b kp ] = Â(õ kp ) ≠ Â Q a ÿ poeP õkp R b , E [log c ki ] = Â( 'ki ) ≠ Â Q a ÿ ioe{0,1} 'ki R b .
where Â(•) is the digamma function. Normal distribution properties lead to 'k oe K, 'p oe P :

E Ë µ t kp È = μt kp , E Ë (µ t kp ) 2 È = V Ë µ t kp È + E Ë µ t kp È 2 = ‡ 2 • ≠1 kp + (μ t kp ) 2 , E [Ê k ] = Ễk , E Ë Ê k Ê T k È = V [Ê k ] + E [Ê k ] E [Ê k ] T = ‡ 2 ˜ k + Ễk ỄT k ,
and Inverse Gamma distribution properties lead to

E 5 1 ‡ 2 6 = ›1 ›2 , E Ë log ‡ 2 È = log ›2 ≠ Â( ›1 ) .
Using all these properties, the following expectations can be calculated as 'j oe J , 'k oe K, 'p oe P :

E C (x tj ≠ µ t kp ) 2 ‡ 2 D = ›1 1 x tj ≠ μt kp 2 2 ›2 + • ≠1 kp , E C (x tj ≠ (t j ) T Ê k ) 2 ‡ 2 D = ›1 1 x tj ≠ (t j ) T Ễk 2 2 ›2 + Trace 1 (t j ) T ˜ k (t j ) 2 .

Complete model

A model integrating parabolic data, piecewise parabolic data and mixed data is now presented. By taking into consideration any types of available data, the resulting model can fit data better and can estimate more accurate clusters. First, data formalism and assumptions are detailed.

Then, the resulting mixture model and its inference procedure are developed.

Data and assumptions

In this part, data consist of J pulses gathering J amplitudes x t = (x tj ) joeJ associated to J times of arrival t = (t j ) joeJ , J continuous features x q = (x qj ) joeJ and J categorical features x c = (x cj ) joeJ from K distinct emitters. Let x j = (x qj , x cj , x tj ) the j th observation vector of mixed variables where

• x qj oe R d is a vector of d continuous radar features such as the Radio Frequency, the Pulse Width, the Azimuth or the Pulse Repetition Interval,

• x cj = ! x cj 0 , . . . , x cj q≠1 "
oe C q is a vector of q categorical radar modulations such as intrapulse modulations or pulse-to-pulse modulations,

• x tj oe R is a continuous variable modeling the Amplitude.

For each pulse j, the temporal evolution variable x tj and mixed variables (x qj , x cj ) are assumed to be independent conditionally to each cluster k oe K with (x miss cj , x obs cj ) oe C q miss j ◊ C q obs j and q miss j + q obs j = q .

'j oe J , (x q , x c )|z j = k |= x t |z j = k . ( 4

Mixture model

According to the independence assumption (4.42), the distribution of mixed data (Chapter 3) and the parabolic and piecewise parabolic relations between temporal evolution data, the component distribution results in

'j oe J , p(x qj , x cj , x tj |z j = k) = p(x qj , x cj |z j = k)p(x tj |z j = k)
where

'j oe J , p(x qj , x cj |z j = k) = Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 22 " c x cj , p(x tj |t j , z j = k) = 1 c k0 N 1 (t j ) T Ê k , ‡ 2 22 " 0 w j Q a c k1 Ÿ poeP 1 b kp N 1 x tj |µ t kp , ‡ 2 22 " p y j R b " 1 w j (4.43)
with

• u = (u j ) joeJ the scale latent variables handling outliers for quantitative data x q and distributed according to a Gamma distribution with shape and rate parameters (-, -) = (kc ,kc ) (k,c)oeK◊Cq conditionally to categorical data x c and labels z = (z j ) joeJ ,

• (µ, ) = ((µ kc ) coeCq , k ) koeK the mean and the variance parameters of quantitative data

x q for each cluster,

• fi = (fi k ) koeK the weights of the multivariate Categorical distribution of categorical data

x c for each cluster,

• y = (y j ) joeJ the latent variables indicating the p th piecewise temporal evolution data x t belong to,

• w = (w j ) joeJ the latent variables indicating if temporal evolution data x t are distributing according to a parabolic relation or a piecewise parabolic relation,

• b = ((b kp ) poeP ) koeK the weights of the Categorical distribution of latent variables y,

• c = (c k0 , c k1 ) koeK the weights of the Categorical distribution of latent variables w,

• µ t = ((µ t kp ) poeP ) koeK the set of piecewises for temporal evolution data x t for each cluster,

• Ê = (Ê k ) koeK the regression parameters for temporal evolution data x t for each cluster,

• ‡ 2 the variance of the measurement noise related to temporal evolution data x t .

Recalling that p(z j = k) = a k where a = (a k ) koeK are the weights related to component distributions, the mixture model is obtained from (4.43) such that 'j oe J , p(x j , u j , y j |t j , ) =

ÿ koeK a k 1 c k0 N 1 x tj | (t j ) T Ê k , ‡ 2 22 " 0 w j Q a c k1 Ÿ poeP 1 b kp N 1 x tj |µ t kp , ‡ 2 22 " p y j R b " 1 w j ◊ Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c x cj (4.44)
where = (a, b, c, µ t , Ê, ‡ 2 , fi, -, -, µ, ) is the set of parameters.

Bayesian framework

As in previous chapters, a Bayesian framework is used to estimate parameters . Assuming datasets (x = (x q , x c , x t ), t) of i.i.d observations (x j = (x qj , x cj , x tj ), t j ) joeJ , independent labels z = (z j ) joeJ , indicators (y, w) = (y j , w j ) joeJ and scale latent variables u = (u j ) joeJ , the complete likelihood associated to (4.44) is defined by

p(x, y, u, z|t, , K) = Ÿ joeJ Ÿ koeK 3 a k 1 c k0 N 1 x tj | (t j ) T Ê k , ‡ 2 22 " 0 w j Q a c k1 Ÿ poeP 1 b kp N 1 x tj |µ t kp , ‡ 2 22 " p y j R b " 1 w j ◊ Ÿ coeCq 1 fi kc N 1 x qj |µ kc , u ≠1 j k 2 G(u j |-kc , -kc ) 2 " c x cj 4 " k z j
.

Eventually, the prior distribution required for is chosen as

p( |K) = p(a|K)p(b|K)p(c|K)p(µ t | ‡ 2 , K)p(Ê| ‡ 2 , K)p( ‡ 2 )p(fi|K)p(-, -|K)p(µ, |K) K J " 0 0 µ 0 ÷ 0 k µ k Ê 0 0 Ê k µ t 0 • 0 µ t k › 0 1 › 0 2 ‡ 2 fi 0 fi k -k -k u j x miss qj x obs qj x miss cj x obs cj z j w j y j c ' 0 b o 0 a Ÿ 0 p 0 r 0 q 0 s 0 x tj (t j )
Figure 4.20: Graphical representation of the proposed model integrating temporal evolution data and mixed-type data. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x qj , x cj , x tj , t j ) decomposed into temporal evolution data (x tj , t j ) and mixed-type data (x qj , x cj ), the scale variables u j and the indicator variables (y j , w j , z j ). Known quantities, respectively unknown quantities, are in blue, respectively in red.

where

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ p(a|K) = D(a|Ÿ 0 ) , p(b|K) = Ÿ koeK D(b k |o 0 ) , p(c|K) = Ÿ koeK D(c k |' 0 ) , p(fi|K) = Ÿ koeK D(fi k |fi 0 ) , p(µ, |K) = Ÿ koeK Ÿ coeCq N 1 µ kc |µ 0 , ÷ ≠1 0 k 2 IW( k |" 0 , 0 ) , p(-, -|K) = Ÿ koeK Ÿ coeCq p(-kc , -kc |p 0 , q 0 , s 0 , r 0 ) , p(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |µ t 0 , • ≠1 0 ‡ 2 2 , p(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k |Ê 0 , ‡ 2 0 2 , p( ‡ 2 ) = IG( ‡ 2 |› 0 1 , › 0 2 ) .
Graphical representation of the proposed model is shown in Figure 4.20.

As previously, a factorized posterior distribution q(x miss q , u, x miss c , y, w, z, |K) = q(x miss q , u, x miss c , y, w, z|K)q( |K) is chosen as an approximation of the intractable posterior joint distribution p(x miss q , u, x miss c , y, w, z, |x obs , K) such that latent variables h = (x miss q , u, x miss c , y, w, z) and parameters are a posteriori independent and q(h|K) = q(x miss q |u, x miss c , z, K)q(u|x miss c , z, K)q(x miss c |z, K)q(y|w, z, K)q(w|z, K)q(z|K) , q( |K) = q(a|K)q(b|K)q(c|K)q(µ t | ‡ 2 , K)q(Ê| ‡ 2 , K)q( ‡ 2 )q(fi|K)q(-, -|K)q(µ, |K) .

According to VB assumptions, the following conjugate variational posterior distributions are obtained from the VB procedure The lower bound on the log evidence is still required to master the VB inference and can be also decomposed into terms related to temporal evolution data (blue terms), mixed data (green terms) and labels z (red terms). This decomposition is obtained as follows

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ q(x miss q |u, x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC N 3 x miss qj | μx miss q jkc , u ≠1 j ˜ x miss q k 4 " c x cj " k z j , q(u|x miss c , z, K) = Ÿ joeJ Ÿ koeK Ÿ coeC G 1 u j |-jkc , -jkc 2 " c x cj " k z j , q(x miss c |z, K) = Ÿ joeJ Ÿ koeK MC(x miss cj |r x miss c jk ) " k z j q(y|z, K) = Ÿ joeJ Ÿ koeK Cat(y j |r y jk ) " 1 w j " k z j , q(w|z, K) = Ÿ joeJ Ÿ koeK Cat(w j |r w jk ) " k z j , q(z|K) = Ÿ joeJ Cat(z j |r j ) , q(a|K) = D(a|Ÿ) , q(fi|K) = Ÿ koeK D(fi| fk ) , q(µ, |K) = Ÿ koeK Ÿ coeC N 1 µ kc | μkc , ÷≠1 kc k 2 IW( k |" k , ˜ k ) , q(-, -|K) = Ÿ koeK p(-kc , -kc |p k , qk , sk , rk ) , q(b|K) = Ÿ koeK D(b k |õ k ) , q(c|K) = Ÿ koeK D(c k | 'k ) , q(µ t | ‡ 2 , K) = Ÿ koeK Ÿ poeP N 1 µ t kp |μ t kp , • ≠1 kp ‡ 2 2 , q(Ê| ‡ 2 , K) = Ÿ koeK N 1 Ê k | Ễk , ‡ 2 ˜ k 2 , q( ‡ 2 ) = IG( ‡ 2 | ›1 , ›2 ) . ( 4
'j oe J , 'k oe K, E [log a k ] = Â(Ÿ k ) ≠ Â Q a ÿ koeK Ÿk R b , log fl t jk = log
L(q|K) = E h, [log p(x, u, y, w, z, |t, K)] ≠ E h, Ë log q(x miss q , x miss c , u, y, w, z, |K) È
where the free energy can be developed as fi,µ,|K)] and the entropy as fi,µ

E h, [log p(x, u, y, w, z, |t, K)] = E h, [log p(x, u, y, w, z|t, , K)] + E Ë log p(b, c, Ê, µ t , ‡ 2 |K) È + E [log p(a,
E h, Ë log q(x miss q , x miss c , u, y, w, z, |K) È = E h, Ë log q(b, c, Ê, µ t , ‡ 2 |K) È + E h, [log q(y, w|z, K)] + E h, Ë log q(x miss q , x miss c , u|z, K) È + E h, [log q(z|K)] + E h, [log q(a,
, |K)] .
Blue terms, respectively green terms, have been previously detailed in subsection 4.3.2 , respectively in Chapter 3. As for red terms, they are detailed below :

E h, [log p(x, u, y, w, z|t, , K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log fl jk , E h, [log q(z|K)] = ÿ joeJ ÿ koeK E h Ë " k z j È log rjk .

Experiments

Two experiments are carried out to evaluate clustering performance with respect to a set of synthetic data. In the first experiment, only temporal evolution data are taken into consideration in the clustering procedure. Then, both temporal evolution data and quantitative data are considered in the second one. For comparison, the spectral clustering [START_REF] Von | A tutorial on spectral clustering[END_REF] and the k-means algorithm from [START_REF] John | Algorithm AS 136: A k-means clustering algorithm[END_REF] are also evaluated. First, characteristics of data, comparison algorithms and evaluation metrics are detailed. Then, both experiments are described and performance are shown to exhibit the e↵ectiveness of the proposed model.

Data, algorithms and metrics

Synthetic data are composed of temporal evolution data related to amplitudes which are distributed according to either a parabolic relation or a piecewise parabolic relation and quantitative data related to continuous radar features which are jointly distributed according to a multivariate normal distribution. Temporal evolution data are generated by sampling a set of data from 2 parabolas and 2 piecewise parabolas directed by

Ê = Q c a ≠1 ≠2 ≠3 ≠4 1 2 3 4 1 2 3 4 R d b .
For each piecewise parabola, p = 4 piecewises are obtained by dividing the time interval in p equal subintervals and assigning to each piecewise the value of the parabola at the minimal time of its related time subinterval. Quantitative data are generated by sampling a set of data from Except for the k-means algorithm, an initialisation is required for clustering algorithms that are involved in these experiments. The similarity graph required for the spectral clustering is obtained from a k-nearest neighbor graph as suggested in [START_REF] Von | A tutorial on spectral clustering[END_REF] where the number of neighbors k is chosen as the product of the log number of observations and the number of clusters. As for the proposed model, a supervised initialisation is retained due to its sensitivity to initialisation. First, prior hyperparameters › 0 1 and › 0 2 are initialised such that the prior mean Quantitative feature 2 • 0 Ÿ 0 ÷ 0 " 0 p 0 r 0 q 0 s 0 (0, 0, 0) T I 3 1 0.5 100 1 1 1 1 1 the variance parameter ‡ 2 is equal to the inverse of the determinant of the covariance matrix of temporal evolution data points. This choice is motivated by the fact that the determinant of the covariance matrix can be interpreted as the generalized variance that reflects the overall spread of the data. Setting › 0 2 = 1, › 0 1 is initialised as the inverse of the generalized variance of the sample of temporal evolution data. In addition, prior piecewise means µ t 0 are initialised from results of a k-means algorithm on temporal evolution data. Then, prior component means µ 0 , respectively covariance matrices 0 , are initialised from results of a k-means algorithm on quantitative data, respectively from diagonal matrices whose diagonal elements are variances of quantitative data. Other hyper-parameters are initialised as in Table 4.9.

E[ 1 ‡ 2 ] = › 0 1 › 0 2 of -0.
Performance on synthetic data are evaluated through the Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF] that compares estimated partitions of data with the ground-truth and the Silhouette Coe cient [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF] which does not require the ground-truth and provides a higher score when clusters are dense and well separated.

Experiments and results

The first experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K when only temporal evolution data are taken into consideration. According to datasets visualised in Figure 4.21, K is set to 4 for synthetic data. Results of the first experiment on synthetic data are shown in Figure 4.23 and in Table 4.10. When ‡ 2 oe {0.0001, 0.01}, the proposed model and the spectral clustering have similar performance in clustering synthetic data with an ARI equals to 0.73 while the k-means algorithm has the lowest performance (ARI = 0.35) in creating convex and isotropic clusters that cannot handle the parabolic and piecewise parabolic structures of the generated data. This limitation is emphasized by higher Silhouette Coe cients of the k-means algorithm whereas the non-convexity of the data is confirmed by the lower Silhouette Coe cients of the ground-truth. Even if all algorithms poorly perform when data are embedded in noise ( ‡ 2 = 0.25), the proposed algorithm estimates clusters with a more parabolic shape than other algorithms which build more isotropic clusters (Subfigures (f), (i) and (l) in The second experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K when all types of data are taken into consideration. The number of clusters K is still set to 4 for synthetic data.. Results of the second experiment on synthetic data are shown in Figure 4.24 and in Table 4.11. All algorithms succeed in clustering synthetic data for ‡ 2 oe {0.0001, 0.01, 0.25} since the ground-truth partition is recovered in Figure 4.24 with an ARI equals to 1 visible on Table 4.11. Adding quantitative information enables algorithms to recover the ground-truth for any value of ‡ 2 .

Conclusion

In this chapter, two types of scanning behaviours have been presented. They can be observed in data when amplitudes and times of arrival shared either a parabolic relation or a piecewise parabolic relation. Since scanning behaviours fully characterise radar emitters, they can be taken into consideration to cluster radar emitters. Hence, mixture models handling parabolic data and piecewise parabolic data have been developed. Three approaches have been investigated : the first one assumes that temporal evolution data are only distributed according to parabolic relations, the second one mainly focuses on piecewise parabolic data and the last one handles both types of relations. Moroever in each approach, the proposed mixture model is enhanced with the mixture model designed for mixed data in Chapter 3 in order to improve clustering performance. Then, parameter estimation has been derived from the Variational Bayesian inference and experiments on real and synthetic data have exhibited the e↵ectiveness of these three approaches.

Chapter 5

Conclusion and perspectives

Radar emitter identification is a crucial function of ESM systems since it prevents enemy forces from surprise attacks by detecting enemy radar signals and it improves military databases by analyzing unknown signals. Depending on radar emitter function and geopolitical context, radar emitters can emit complex signals based on pulse-to-pulse modulation patterns. Radar signal patterns can be decomposed into continuous features given by the continuous parameters of radar pattern pulses and categorical features that represent modulations of pulse sequences. Furthermore, radar signals are often partially observed in the electromagnetic environment due to failures of deinterleaving techniques or sensors deficiencies. Therefore, a framework handling any types of data has been developed in this work to perform classification and clustering of radar emitters even in presence of outliers and missing data.

State-of-the-art algorithms have been reviewed in Chapter 1. They perform either classification or clustering by learning boundaries that separate data into heterogeneous clusters or by learning underlying structure of data to constitute clusters. However, neither of these algorithms provides an internal framework that infers on missing data and copes with any degrees of supervision or any types of data. Therefore, an approach based on mixture models has been proposed. Theoretical aspects of mixture models have been introduced by detailing ways of modeling and estimating them for a generic type of data. Indeed, mixture models benefit from a flexible and probabilistic framework to handle outliers and missing data by introducing a latent space where each latent variable focuses on a specific constraint. However, the resulting model is not tractable and model learning is processed through Variational Bayes Approximation. Eventually whatever degree of supervision is required, the number of classes K and parameters can be estimated to perform classification and clustering tasks. Nonetheless, the Variational Bayes Appromixation tends to under-estimate uncertainties related to variational estimation. To this end, the Expectation propagation algorithm [START_REF] Thomas | Expectation propagation for approximate Bayesian inference[END_REF], which minimizes the reversed Kullback-Leibler divergence, can be implemented in future works to improve estimation accuracy.

Chapter 2 has focused on classification and clustering of continuous data with a scale mixture of Normal distributions accounting for missing data and outliers. Benefiting from Gaussian properties and the introduction of latent variables, the proposed model has shown its e ciency for inferring on missing data, performing classification and clustering tasks and selecting the correct number of clusters in a dataset obtained from an experimental protocol generating realistic data. A major contribution in this chapter is the incorporation of latent variables handling missing data and provided with a variational posterior distribution that leads to a more e↵ective inference on missing data. As pointed out in both experiments, the e↵ectiveness of the proposed model results from the fact that standard missing data imputation methods can create outliers that deteriorate performance of classification and clustering algorithms whereas in the proposed model inference on missing data and labels prediction are jointly estimated. Indeed, embedding the inference procedure into the model framework allows properties of the model, such as outliers handling, to counterbalance drawbacks of imputation methods such as outlier creation. As for outliers handling, a full Bayesian approach has been adopted to avoid using the deterministic variable ‹ k parametrizing the distribution of the latent variable u modeling outliers. This full Bayesian treatment enables a less restrictive modeling of data since parameters of the latent variable u are estimated into the Variational Bayes Approximation framework instead of being updated via an optimization procedure. Despite these advantages, the proposed model has a higher computational cost than comparison algorithms especially during the model learning step. Hence, a parallelization of the proposed model would be useful in order to reduce its computational burden.

As for Chapter 3, it has presented the general case where both continuous and categorical data are used for classification and clustering tasks. This chapter has precisely focused on modeling dependencies between continuous and categorical data in order to infer on missing data while performing classification and clustering. To this end, an approach based on the Location Mixture Model has been investigated on by establishing conditional relations between continuous and categorical data to tackle issues related to outliers and missing data. The developed approach has exhibited its e↵ectiveness for inferring on missing data, performing classification and clustering tasks and selecting the correct number of clusters even for high proportions of missing values. As pointed out in Chapter 2, the proposed approach enables joint estimation of missing components and labels by embedding the inference procedure into the model framework. Indeed, estimating the missing components conditionally to their labels proves to be more e↵ective than standard imputation methods which do not take into consideration information related to the data partition. Moreover, experiments have pointed out that using continuous and categorical data can really improve classification and clustering performance than considering either only continuous data or only categorical data. Indeed, higher performance on mixed data lie in a more relevant separation of clusters obtained by taking advantage of the more complex structure of mixed-type data. In this work, the Location Mixture Model assumption has been naturally considered since radar pattern designers use to define pulse modulation sequences (categorical radar data) before assigning pulse values (continuous radar data). However in the Electronic Warfare context, continuous radar data are first measured by sensors before being processed to deduce categorical data related to modulation patterns. Therefore, it would be interesting to investigate the Underlying Response Variable approach [START_REF] Everitt | A finite mixture model for the clustering of mixed-mode data[END_REF][START_REF] Brian | The clustering of mixed-mode data: a comparison of possible approaches[END_REF] to assess assumption on conditioning categorical radar data by continuous radar data.

In Chapter 4, two types of scanning behaviours have been presented. They can be observed in radar data when amplitudes and times of arrival of radar emitters share either a parabolic or a piecewise parabolic relation. Since scanning behaviours fully characterise radar emitters, they have been taken into consideration to cluster radar emitters. To this end, both types of relation have been integrated into the mixture model framework by modeling the parabolic relation as a Bayesian regression and the piecewise parabolic relation as a mixture of normal distributions. As in Chapters 2 and 3, estimation of parameters has been processed through the Variational Bayesian inference and experiments on real and synthetic data have exhibited the e↵ectiveness of the proposed model. Indeed, the resulting model enables creation of non isotropic clusters which better fit parabolic data than standard algorithms as the k-means algorithm. Nonetheless, when radar emitters share similar scanning behaviours, namely when amplitude parabolas of radar emitters intersect during their scanning period, the proposed model cannot perfectly identify radar emitters. Therefore, the proposed mixture model has been enhanced with the mixture model designed for mixed data in Chapter 3 in order to improve clustering performance. Then, the complete model has managed to separate radar emitters in real operational cases by taking advantage of the whole available information related to radar emitters. In this work, only radar emitters presenting parabolic scanning behaviours have been clustered in real data cases. Hence, piecewise parabolic scanning behaviours remain to be evaluated on real data cases in future practical studies. Moreover, a radar signal pattern could be modeled as a Markovian process through its temporal evolution, continuous and categorical features since it is mainly defined as a pattern of pulses whose features share sequential and conditional relations. To this end, a clustering method, that integrates a Markovian process while handling missing data and outliers, could be developed in future works.

To conclude, the di↵erent mixture models developed in this work have focused on performing classification and clustering on various types of real and simulated data while handling outliers and missing values. These models have managed to reach high performance in classification and clustering tasks even in presence of large proportions of missing data and have proposed an e↵ective inference procedure to reconstruct missing features. These models have been assessed on datasets gathering around 50 radar emitters according to the di↵erent experiments. Nonetheless, real databases may contain thousands of radar emitters. Therefore, the proposed models have to be tested on larger databases before being integrated into industrial and operational products.

Dans le contexte de la Guerre Electronique, l'identification des signaux radar est un atout majeur de la prise de décisions tactiques liées au théâtre d'opérations militaires. En fournissant des informations sur la présence de menaces, la classification et le partitionnement des signaux radar ont alors un rôle crucial assurant un choix adapté des contre-mesures dédiées à ces menaces et permettant la détection de signaux radar inconnus pour la mise à jour des bases de données. Les systèmes de Mesures de Soutien Electronique enregistrent la plupart du temps des mélanges de signaux radar provenant de di↵érents émetteurs présents dans l'environnement électromagnétique. Le signal radar, décrit par un motif de modulations impulsionnelles, est alors souvent partiellement observé du fait de mesures manquantes et aberrantes. Le processus d'identification se fonde sur l'analyse statistique des paramètres mesurables du signal radar qui le caractérisent tant quantitativement que qualitativement. De nombreuses approches mêlant des techniques de fusion de données et d'apprentissage statistique ont été développées. Cependant, ces algorithmes ne peuvent pas à la fois e↵ectuer la classification ainsi que le partionnement des émetteurs radar et gérer les données manquantes. A cette fin, des méthodes de substitution de données sont requises en amont de la classification et du partitionnement mais leur utilisation entraine l'apparition de nouvelles valeurs aberrantes. L'objectif principal de cette thèse est alors de définir un modèle de classification et partitionnement intégrant la gestion des valeurs aberrantes et manquantes présentes dans tout type de données. Une approche fondée sur les modèles de mélange de lois de probabilité est proposée dans cette thèse. Les modèles de mélange fournissent un formalisme mathématique flexible favorisant l'introduction de variables latentes permettant la gestion des données aberrantes et la modélisation des données manquantes dans les problèmes de classification et de partitionnement. L'apprentissage du modèle ainsi que la classification et le partitionnement sont réalisés dans un cadre d'inférence bayésienne où une méthode d'approximation variationnelle est introduite afin d'estimer la loi jointe a posteriori des variables latentes et des paramètres. Des expériences sur diverses données montrent que la méthode proposée fournit de meilleurs résultats que les algorithmes standards.

Le premier chapitre de cette thèse présente les di↵érents algorithmes de l'état de l'art en matière de classification et de partitionnement. Tant par l'apprentissage de frontières au sein des données que par l'apprentissage d'une structure sous-jacente des données, ces algorithmes répondent aux problématiques de classification et partitionnement en créant des groupes hétérogènes d'observations. Cependant, aucun de ces algorithmes ne peut à la fois intégrer des données continues et catégorielles, gérer les contraintes liées aux données manquantes et s'adapter à di↵érents degrés de supervision. L'approche fondée sur les modèles de mélange de lois de probabilité est alors introduite afin de palier ces divers problèmes. En e↵et, les modèles de mélange fournissent un cadre probabiliste favorisant l'introduction de variables latentes permettant l'intégration de données de tout type, la gestion de données aberrantes ainsi que la modélisation des données manquantes dans les problèmes de classification et de partitionnement. Néanmoins, cette approche requiert l'utilisation d'une méthode d'approximation bayésienne appelée Variational Bayes et dont les aspects théoriques sont détaillés dans ce chapitre.

Le deuxième chapitre introduit un modèle de mélange de lois gaussiennes afin de prendre en compte les données continues représentant les paramètres physiques des impulsions. Ce modèle de mélange gaussien est ensuite mis à jour via l'introduction de variables latentes modélisant les valeurs aberrantes et manquantes afin d'obtenir un modèle robuste à ces deux contraintes. L'inférence est menée au travers de la méthode Variational Bayes permettant d'obtenir une approximation de la distribution jointe a posteriori des paramètres et des variables latentes du modèle. Enfin, ce modèle est testé sur des données acquises à l'aide d'un protocole expérimental fournissant des données réalistes intégrant les contraintes des systèmes d'acquisition opérationnels.

Le troisième chapitre intègre les données categorielles au modèle précédent en conditionnant les variables continues d'une observation par ses variables catégorielles. Le modèle obtenu est alors un mélange de lois gaussiennes conditionnelles intégrant également des variables latentes modélisant les valeurs manquantes propres aux observations catégorielles. L'inference est à nouveau faite par le biais de l'approximation variationelle bayésienne afin d'obtenir la distribution jointe a posterior des paramètres et variables latentes du modèle. Les performances du modèle proposé sont ensuite testées sur des données générées à partir d'une base de données réelles comportant 55 émetteurs radar avec des motifs impulsionnels variés.

Enfin, le dernier chapitre se focalise sur le caractère temporel des données impulsionnelles. En e↵et, l'évolution temporelle des amplitudes liées aux impulsions d'un émetteur radar présente une forme parabolique qui peut être exploitée afin d'améliorer la classification et le partionnement des émetteurs radar. Dans ce but, cette relation parabolique est modélisée par le biais d'une régression parabolique bayésienne intégrée au modèle de mélange. Les paramètres du modèle sont alors estimés par le biais de la précédente méthode d'approximation variationnelle et le modèle résultant est testé sur des données synthétiques et réelles provenant de di↵érents cas opérationnels.

En conclusion, les di↵érents modèles développés dans cette thèse ont permis la classification et le partionnement d'émetteurs radar caractérisés par des motifs impulsionnels présentant des valeurs manquantes et aberrantes tant continues que catégorielles. Ces modèles ont démontré leur e cacité en réalisant de bonnes performances sur des bases de données synthétiques et réelles même en présence d'une grande proportion de valeurs manquantes. Néanmoins, les bases de données réelles peuvent contenir des milliers d'émetteurs radar et les modèles proposés doivent alors être mis à l'épreuve sur de plus grandes bases de données avant d'être intégrés dans de futurs produits industriels. Here, an approach based on mixture models is developed since mixture models provide a mathematically based, flexible and meaningful framework for the wide variety of classification and clustering requirements. The proposed approach focuses on the introduction of latent variables that give us the possibility to handle sensitivity of the model to outliers and to allow a less restrictive modelling of missing data. A Bayesian treatment is adopted for model learning, supervised classification and clustering and inference is processed through a variational Bayesian approximation since the joint posterior distribution of latent variables and parameters is untractable. Some numerical experiments on synthetic and real data show that the proposed method provides more accurate results than standard algorithms.

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

Introduction 1 1

 1 State of the art and selected approach 5 1.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.1 Supervised learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.2 Unsupervised learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Selected approach : mixture models . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

•

  DBSCAN : Density-Based Spatial Clustering of Applications with Noise • GMM : Gaussian Mixture Model • KM : k-means algorithm • KNN or k-nn : k-nearest neighbors algorithm • LDA : Linear Discriminant Analysis • NN : Neural Networks • RdF : Random Forests • SC : Spectral Clustering • SVM : Support Vector Machines

•

  A : Amplitude • PRI : Pulse Repetition Interval vii • PW : Pulse Width The abbreviations related to Bayesian theory are : • MAP : Maximum A Posteriori • VB : Variational Bayes • VBA : Variational Bayes Approximation • VBE : Variational Bayes Expectation • VBM : Variational Bayes Maximization viii List of Figures 1.1 Linear Discriminant Analysis and Quadratic Discriminant Analysis performed on Iris dataset [Lic13] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Impact of the choice of the number of neighbors on the classification rule . . . . 1.3 A classification rule di↵erently learned by four standard algorithms . . . . . . . . 1.4 Clusters centers are moving as and when iterations of the k-means algorithm are progressing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Clustering performance of state-of-the-art algorithms . . . . . . . . . . . . . . . .2.1 Four pulses from a radar emitter . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Diagram of the acquisition system . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Acquired pulses from a radar emitter where the three features (PRI,PW,RF) are shown on the figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Outliers formation during primary parameters measurement on real data . . . . . 2.5 Presence of outliers in observations of a radar emitter. . . . . . . . . . . . . . . . 2.6 Graphical representation of the standard Gaussian mixture model . . . . . . . . 2.7 Graphical representation of the Gaussian mixture model handling missing data . 2.8 Graphical representation of the Gaussian mixture model handling outliers . . . . 2.9 Graphical representation of the Gaussian mixture model handling outliers with hyper-parameters (-, -) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Graphical representation of the proposed model for continuous data . . . . . . . 2.11 Dataset gathering 6300 observations from 42 radar emitters . . . . . . . . . . . . 2.12 Classification performance on continuous data . . . . . . . . . . . . . . . . . . . . 2.13 Mean-squared errors of missing data imputation methods on continuous data . . 2.14 Evolution of computing times taken by model learning . . . . . . . . . . . . . . . 2.15 Performance of the proposed model compared with DBSCAN on continuous data 2.16 Performance of the proposed model compared with k-means algorithm on continuous data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Estimation of the number of clusters for continuous data . . . . . . . . . . . . . . 2.18 Estimation of the number of clusters by DBSCAN on continuous data . . . . . . 2.19 Evolution of computing times for clustering algorithms . . . . . . . . . . . . . . . 3.1 Acquired pulses from a radar emitter where the three features (PRI,PW,RF) are shown on the figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Examples of di↵erent modulations of parameter values . . . . . . . . . . . . . . . 3.3 Linear Frequency Modulation on a pulse . . . . . . . . . . . . . . . . . . . . . . . 3.4 Phase coding generated from a Zado↵ code. . . . . . . . . . . . . . . . . . . . . . 3.5 Di↵erent scanning types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Graphical representation of the proposed model integrating mixed data . . . . . 3.7 Dataset gathering 5500 continuous observations from 55 radar emitters . . . . . . ix 3.8 Classification performance on mixed data . . . . . . . . . . . . . . . . . . . . . . 3.9 Evaluation of imputation methods and posterior reconstructions on mixed data . 3.10 Performance of the proposed model compared with DBSCAN on mixed data . . 3.11 Performance of the proposed model compared with k-means algorithm on mixed data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Estimation of the number of clusters on mixed data . . . . . . . . . . . . . . . . 3.13 Estimation of the number of clusters by DBSCAN on mixed data . . . . . . . . . 4.1 Simulated data where amplitudes x t and times of arrival t are distributed according to a parabolic relation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Graphical representation of the proposed mixture model handling parabolic data 4.3 Graphical representation of the proposed model integrating temporal evolution data and mixed-type data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Synthetic parabolic data generated from di↵erent values of the variance parameter ‡ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Synthetic quantitative data generated from 4 multivariate normal distributions . 4.6 Real data obtained from 3 operational cases . . . . . . . . . . . . . . . . . . . . . 4.7 Results on synthetic data during the first experiment when only temporal evolution data are considered . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Results on real parabolic data during the first experiment when only temporal evolution data are considered . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Results on synthetic parabolic data when all types of data are considered . . .. 4.10 Results on real parabolic data when any types of data are considered . . . . . . . 4.11 Simulated data where amplitudes x t and times of arrival t are distributed according to a piecewise parabolic relation defined with P = 4 piecewises. . . . . . . . . 4.12 Graphical representation of the proposed mixture model handling piecewise parabolic data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Graphical representation of the proposed model integrating temporal evolution data and mixed-type data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Synthetic piecewise parabolic data generated from di↵erent values of the variance parameter ‡ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Synthetic quantitative data generated from 4 multivariate normal distributions . 4.16 Results on synthetic parabolic piecewise data when only temporal evolution data are considered . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Results on synthetic data when all types of data are considered . . . . . . . . . . 4.18 Data generated from two distinct emitters presenting a parabolic scanning behaviour and a piecewise parabolic scanning behaviour . . . . . . . . . . . . . . . 4.19 Graphical representation of the proposed mixture model handling parabolic and piecewise parabolic data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Graphical representation of the proposed model integrating temporal evolution data and mixed-type data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Synthetic parabolic and piecewise parabolic data generated from di↵erent values of the variance parameter ‡ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Synthetic quantitative data generated from 4 multivariate normal distributions . 4.23 Results on synthetic data when only temporal evolution data are considered . . . 4.24 Results on synthetic data when any types of data are considered . . . . . . . . . x List of Tables 2.1 Initialisation of hyper-parameters values for classification on continuous data . . 2.2 Initialisation of hyper-parameters values for clustering on continuous data . . . . 3.1 All known binary Barker codes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Categorical features observed in categorical data . . . . . . . . . . . . . . . . . . 3.3 Initialisation of hyper-parameters values for classification on mixed data . . . . . 3.4 Initialisation of hyper-parameters values for clustering on mixed data . . . . . . . 4.1 Initialisation of hyper-parameters values for clustering on parabolic data . . . . . 4.2 Clustering performance on synthetic parabolic data . . . . . . . . . . . . . . . . . 4.3 Clustering performance on real parabolic data . . . . . . . . . . . . . . . . . . . . 4.4 Clustering performance on synthetic parabolic and mixed data . . . . . . . . . . 4.5 Clustering performance on real parabolic and mixed data . . . . . . . . . . . . . 4.6 Initialisation of hyper-parameters values for clustering on piecewise parabolic data 4.7 Clustering performance on synthetic piecewise parabolic data . . . . . . . . . . . 4.8 Clustering performance on synthetic piecewise parabolic and mixed data . . . . . 4.9 Initialisation of hyper-parameters values for clustering on parabolic and piecewise parabolic data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Clustering performance on synthetic parabolic and piecewise parabolic data . . . 4.11 Clustering performance on synthetic parabolic, piecewise parabolic and mixed data xi

Contents 1 . 1

 11 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.1 Supervised learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.2 Unsupervised learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 Selected approach : mixture models . . . . . . . . . . . . . . . . . . . 15 1.2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.2.2 Latent variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.2.3 Bayesian framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.2.4 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.2.5 Classification and clustering . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Figure 1 . 1 :

 11 Figure 1.1: Linear Discriminant Analysis and Quadratic Discriminant Analysis performed on Iris dataset [Lic13]

(1. 4 )

 4 k-nearest neighbors k-nearest neighbors (k-NN) algorithm belongs to the family of instance-based learning algorithms [AKA91] which are non-parametric general algorithms that classify a new unlabelled observation

Figure 1 . 2 :

 12 Figure 1.2: Impact of the choice of the number of neighbors on the classification rule

Figure 1

 1 Figure 1.3: A classification rule di↵erently learned by four standard algorithms

  Procedure 1.2 K-means algorithm Input: Unlabeled dataset x oe R dúN and number of clusters K Output: Partition of z and cluster centroids µ 1 , µ 2 , . . . , µ K oe R d Initialise cluster centroids µ 1 , µ 2 , . . . , µ K randomly repeat Assign each data point to its closest cluster centroid : 'n oe {1, . . . , N}, z n = arg min k ||x n ≠ µ k ||

Figure 1 . 4 :

 14 Figure 1.4: Clusters centers are moving as and when iterations of the k-means algorithm are progressing.

  Figure 1.5: Clustering performance of k-means algorithm in Figure (a), DBSCAN in Figure (b) and spectral clustering in Figure (c) for various datasets by using Scikit-learn library [PVG + 11].

Figure 2 . 1 :Figure 2 . 2 :

 2122 Figure 2.1: Four pulses from a radar emitter whose related amplitudes and times are shown on the vertical axis and the horizontal axis. Then, the continuous features PRI and RF are obtained by delimiting each pulse according to their TOA and PW. Radar Waveform Emitter Signal Receiver Laptop

Figure 2 . 4 :

 24 Figure 2.4: Outliers formation during primary parameters measurement on real data. (a), (c) and (e) are respectively exact TOA, PW and PRI. (b), (d) and (f) are outliers for TOA, PW and PRI.

Figure 2

 2 Figure 2.5: Presence of outliers in observations of a radar emitter.

  modeled by a latent variable x miss j oe R d miss j such that 1 AE d obs j AE d and d miss j = d ≠ d obs .

Figure 2 . 11 :

 211 Figure 2.11: Dataset gathering 6300 observations from 42 radar emitters. Some clusters are completely separable whereas some others share features and can not be linearly separated.

Procedure 2 . 1

 21 Classification procedure on continuous data : Training step Input: Training set x train and associated labels z train Output: Learned parameters ˜ train

Figure 2 .Figure 2 .Figure 2 . 13 :

 22213 Figure2.12: Classification performance are presented for the proposed model (PM) in blue, the NN in red, the RnF in green and the KNN in cyan. The solid lines represent the average accuracies with discarded observations for the NN, the RnF and the KNN, the dashed lines stands for the average accuracies with mean imputation for the NN, the RnF and the KNN whereas the doted lines shows average accuracies with KNN imputation. Shaded error regions represent standard deviations of accuracies.

Figure 2 .

 2 Figure 2.14: Evolution of computing times taken by model learning for Random Forests (RF), K nearest neighbors algorithm (KNN), Neural Network (NN) and the proposed model (PM).

Figure 2 .Figure 2 .Figure 2 .

 222 Figure 2.15: Performance of the proposed model compared with DBSCAN for K = 42 according to di↵erent proportions of missing values and imputation methods on realistic data.

Figure 2 .

 2 Figure 2.18: Estimation of the number of clusters by DBSCAN according to mean imputation, k-NN imputation and posterior reconstruction of the proposed model.

Figure 2 .

 2 Figure 2.19: Evolution of computing times for DBSCAN, k-means algorithm (KM) and the proposed model (PM).

Contents 3. 1

 1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.1.1 Assumptions on continuous data . . . . . . . . . . . . . . . . . . . . . . 52 3.1.2 Assumptions on categorical data . . . . . . . . . . . . . . . . . . . . . . 53 3.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.2 Assumptions on mixed data . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.2.3 Proposed model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.3.1 Variational posterior distributions . . . . . . . . . . . . . . . . . . . . . 65 3.3.2 VBE-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.3.3 VBM-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3.4 Lower Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.3.5 Expectations from variational distributions . . . . . . . . . . . . . . . . 72 3.4 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.4.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.4.2 Classification experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.4.3 Clustering experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

Figure 3

 3 Figure 3.1: Acquired pulses from a radar emitter where the three features (PRI,PW,RF) are shown on the figure.

j=µ obs 2 ,

 2 and d miss j + d obs j = d , where R d miss j and R d obs j , are disjoint subsets of R d embedding missing features x miss qj and observed features x obs qj . Then, properties of the multivariate normal distribution leads to obtain two normal distributions for observed and missing features such that 'j oe J , µ miss + cov obs ≠1 1 x obs qj ≠

Figure 3 . 2 :

 32 Figure 3.2: Examples of di↵erent modulations of parameter values. Figure (a) shows no modulation of parameter values. Figure (b) introduces sliding values on sequences of 5 pulses.Figure (c) exhibits a sinusoidal wobulated emission. Figure (d) presents a dwelled emission shaped with 5 piecewises. Figure (e) shows a staggered emission composed of 5 moments distributed according to a sequence mapped on 10 pulses. At last, Figure (f) presents a jittered emission where parameter values are normally distributed.

  Time domain representation of a Linear Frequency Modulation signal. (dB/Hz) (b) Frequency domain representation of a Linear Frequency Modulation signal.

Figure 3

 3 Figure 3.3: Linear Frequency Modulation on a pulse. The Figure (a), respectively Figure (b), shows a time domain representation, respectively a frequency domain representation, of a linear frequency modulation signal.

  Figure 3.4 exhibits an example of pulse compression with a Zado↵ code.

Figure 3 . 4 :

 34 Figure 3.4: Phase coding generated from a Zado↵ code.

Figure 3

 3 Figure 3.5: Di↵erent scanning types from [AWT18]. Sketches on the left illustrate scanning behaviours over the scanning period. Graphics on the right represent the evolution of pulse amplitudes over the scanning period.

  miss ,c obs x miss cj ,x obs cj c miss ,c obs where " c miss ,c obs x miss cj ,x obs cj = I 1 if x miss cj = c miss and x obs cj = c obs 0 otherwise J = " c miss x miss cj ◊ " c obs x obs cj . (3.2) By using the previous equality (3.2), the multivariate categorical distribution in (3.1) becomes p(x miss cj , x obs cj ) = Ÿ c miss oeC q miss j marginal distribution for observed features x obs cj and a conditional distribution for missing features x miss cj are obtained from (3.3) such that

  , a multivariate categorical distribution for missing features x miss cj conditionally to observed features x obs cj is deduced from (3.4) where '(c miss , c obs ) oe C q miss j ◊ C q obs j : p(x miss cj = c miss |x obs cj = c obs ) = fi c miss ,c obs q c miss oeC q miss j fi c miss ,c obs

  miss qj , x obs qj ) oe R d miss j ◊ R d obs j

≠1 .

 ≠1 Regarding missing categorical data, they are handled as in subsection 3.1.2 such that missing features x miss cj follow a multivariate categorical distribution conditionally to observed features x obs cj where '(c miss , c obs ) oe C q miss j ◊ C q obs j : p(x miss cj = c miss |x obs cj = c obs ) = fi c miss ,c obs q c miss oeC q miss j fi c miss ,c obs with fi c miss ,c obs the joint probability fi c defined in (3.5) and (3.6) for c = (c miss , c obs ).

Figure 3 . 6 :

 36 Figure3.6: Graphical representation of the proposed model integrating mixed data. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x qj , x cj ) joeJ decomposed into quantitative data (x qj ) joeJ and categorical data (x cj ) joeJ , the scale variables u j and the indicator variables z j . Known quantities, respectively unknown quantities, are in blue, respectively in red

s

  p -≠1 e ≠q--s- (-) r I ->0 I ->0 ˆ-ˆ-. Graphical representation of the proposed model is shown in Figure 3.6.

  j, k, c miss ) oe J ◊ K ◊ C q miss j : log fl x miss cj kc miss = ÿ c obs oeC q obs j " c obs x obs cj log fl x cj kc miss c obs (3.17) Hence, a multivariate categorical distribution is deduced for x miss c from (3.16) conditionally to labels z such that qj,k)oeJ ◊K are obtained from (3.17) where '(j, k, c) oe J ◊ K ◊ C q miss j posterior categorical distributions are obtained for labels z by marginalising over x miss c in (3.16) such that ⁄ E [log p(x, u, z| , K)] ˆxmiss q oe J , k oe K, log fl jk = E [log a k ] + log ÿ c miss oeC q miss j

Figure 3 . 7 :

 37 Figure3.7: Dataset gathering 5500 continuous observations from 55 radar emitters. Some clusters are completely separable whereas some others share features and cannot be linearly separated.

  miss qj , u j , x cj z j )ˆu j ˆxcj ˆzj data, reconstruction performance are evaluated through the comparison of Jaccard distances between original categorical data and imputation methods against Jaccard distances [Jac01] between original categorical data and the variational posterior marginal mode of missing categorical data given by 'j oe J , xmiss cj = arg max c miss oeC q miss j ⁄ q(x miss cj , z j )dz j

Procedure 3. 2

 2 Classification procedure on mixed data: Prediction step Input: Unlabelled dataset x pred and learned parameters ˜ train Output: Predicted labels zpred Update -jkc , -jkc , rx miss cj jkc , μx miss q jkc , ˜ x miss qj k , rjk return zpred such that each zpred j = arg max koeK rjk

  Classification on mixed data.

Figure 3 . 8 :

 38 Figure 3.8: Classification performance are presented for the proposed model (PM) in blue, the NN in red, the RnF in green and the KNN in cyan.Figure (a) exhibits classification performance when only continuous data are taken into consideration. Figure (b) exhibits classification performance when only categorical data are taken into consideration. Figure (c) exhibits classification performance when both continuous and categorical data are taken into consideration. For each figure, the solid lines represent accuracies with a posteriori reconstructed missing data for the NN, the RnF and the KNN, the doted dashed lines stands for accuracies with mean/mode imputation for the NN, the RnF and the KNN whereas the dashed lines shows accuracies with KNN imputation.

  Clustering performance for categorical data.

  Clustering performance for mixed data.

Figure 3 .

 3 Figure 3.10: Performance of the proposed model compared with DBSCAN according to di↵erent proportions of missing values and imputation methods. The number of clusters K is fixed at 45 for categorical data and 55 for continuous and mixed data.

  Clustering performance for mixed data.

Figure 3. 11 :

 11 Figure 3.11: Performance of the proposed model compared with k-means algorithm according to di↵erent proportions of missing values and imputation methods. The number of clusters K is fixed at 45 for categorical data and 55 for continuous and mixed data.

Figure 3 .

 3 Figure 3.12: Estimation of the number of clusters using the lower bound (LB) and the silhouette score (S) for the proposed model and only the silhouette score (S) for the k-means algorithm.

Figure 3 . 13 :

 313 Figure 3.13: Estimation of the number of clusters by DBSCAN according to mean imputation, k-NN imputation and posterior reconstruction of the proposed model.

Figure 4 . 3 :

 43 Figure 4.3: Graphical representation of the proposed model integrating temporal evolution data and mixed-type data. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x qj , x cj , x tj , t j ) decomposed into temporal evolution data (x tj , t j ) and mixed-type data (x qj , x cj ), the scale variables u j and the indicator variables z j . Known quantities, respectively unknown

Figure 4 . 4 :

 44 Figure 4.4: Synthetic parabolic data generated from di↵erent values of the variance parameter ‡ 2 . Figures (a), (b) and (c) present unlabeled data where 4 parabolas are generated. Ground-truth are visible on Figures (d), (e) and (f).

Figure 4 . 5 :

 45 Figure 4.5: Synthetic quantitative data generated from 4 multivariate normal distributions. Figure (a) shows unlabeled data and Figure (b) exhibits the ground-truth.

Figure 4 . 2 )

 42 . Indeed the Silhouette Coe cient of the proposed model (S = 0.14) is closer to the Silhouette Coe cient of the ground-truth (S = 0.10) than Silhouette Coe cients of spectral clustering (S = 0.53) and k-means (S = 0.56). Results of the first experiment on real data are shown in Figure 4.8 and in Table

Figure 4 . 7 :

 47 Figure 4.7: Results on synthetic data during the first experiment when only temporal evolution data are considered. Figures (a), (b) and (c) show synthetic data generated with di↵erent values of the variance ‡ 2 . Figures (d), (e) and (f) show clustering results of the proposed model (PM). Figures (g), (h) and (i) show clustering results of the spectral clustering (SC). Figures (j), (k) and (l) show clustering results of the k-means algorithm (KM).

Figure 4 . 8 :

 48 Figure 4.8: Results on real parabolic data during the first experiment when only temporal evolution data are considered. Figures (a), (b) and (c) show real data in di↵erent cases. Figures (d), (e) and (f) show clustering results of the proposed model (PM). Figures (g), (h) and (i) show clustering results of the spectral clustering (SC). Figures (j), (k) and (l) show clustering results of the k-means algorithm (KM).

Figure 4 . 9 :

 49 Figure 4.9: Results on synthetic parabolic data when all types of data are considered. Figures (a), (b) and (c) show synthetic data generated with di↵erent values of the variance ‡ 2 . Figures (d), (e) and (f) show clustering results of the proposed model (PM). Figures (g), (h) and (i) show clustering results of the spectral clustering (SC). Figures (j), (k) and (l) show clustering results of the k-means algorithm (KM).

Figure 4 . 11 :

 411 Figure 4.11: Simulated data where amplitudes x t and times of arrival t are distributed according to a piecewise parabolic relation defined with P = 4 piecewises.

Figure 4 .

 4 Figure 4.12: Graphical representation of the proposed mixture model handling piecewise parabolic data.The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed amplitudes x j and the indicator variables (z j , y j ). Known quantities, respectively unknown quantities, are in blue, respectively in red.

.

  (r y jk ) joeJ are obtained from (4.23) as follows 'j oe J , 'k oe K, 'p oe P, ry jkp = Marginalising over y in (4.22), a categorical is obtained for cluster labels z such that q(z|K) = Ÿ joeJ Cat(z j |r j ) with 'j oe J , 'k oe K, rjk = fl jk q koeK fl jk . where fl jk = E [log a k ] + log ÿ poeP fl y jkp .

Figure 4 . 13 :

 413 Figure 4.13: Graphical representation of the proposed model integrating temporal evolution data and mixed-type data.The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x qj , x cj , x tj , t j ) decomposed into temporal evolution data (x tj , t j ) and mixed-type data (x qj , x cj ), the scale variables u j and the indicator variables (y j , z j ). Known quantities, respectively

x qj |µ kc , u ≠1 j k 2 G

 2 (u j |kc ,kc ) , b, µ t , ‡ 2 ,fi, -, -, µ,) is the set of parameters.

  The red term E [log a k ] is deduced from properties of the Dirichlet distribution, the blue term log fl t jk is deduced from (4.22) and (4.23) in subsection 4.2.2 and the green term log fl qc jk has been detailed in Chapter 3. Hence, z is distributed a posteriori according to a product of Categorical distributions parametrized by r = (r jk ) (j,k)oeJ ◊K given by 'j oe J , 'k oe K, rjk = fl jk q koeK fl jk(4.32)

and 4 .

 4 15 where each radar emitter is represented by a piecewise parabola (Figure4.14) and a Gaussian cluster (Figure4.15).

Figure 4 .

 4 Figure 4.14: Synthetic piecewise parabolic data generated from di↵erent values of the variance parameter ‡ 2 . Figures (a), (b) and (c) present unlabeled data where 4 piecewise parabolas are generated. Groundtruth are visible on Figures (d), (e) and (f).

Figure 4 .

 4 Figure 4.15: Synthetic quantitative data generated from 4 multivariate normal distributions. Figure (a) shows unlabeled data and Figure (b) exhibits the ground-truth.

Figure 4 .

 4 7). Indeed the Silhouette Coe cient of the proposed model (S = 0.14) is closer to the Silhouette Coe cient of the ground-truth (S = 0.10) than Silhouette Coe cients of spectral clustering (S = 0.53) and k-means (S = 0.56).

Figure 4 .

 4 Figure 4.19: Graphical representation of the proposed mixture model handling parabolic and piecewise parabolic data. The arrows represent conditional dependencies between the random variables. The K-plate represents the K mixture components and the J-plate the independent identically distributed observations (x tj , t j ) decomposed into the amplitude x j and the polynomial transformation (t j ) and the indicator

D 4 ,.

 4 log fl w jk1 = E [log c k1 ] + log ÿ poeP fl y jkp . (4.40)Then, a categorical distribution for scanning type labels w is deduced from (4.39) conditionally to cluster labels z such that q(w|z, K) = their parameters (r y jk ) joeJ are obtained from (4.40) as follows 'j oe J , 'k oe K, 'i oe {0, 1}, rw jki =Eventually by marginalising over w in (4.39) a categorical distribution is obtained for cluster labels z such thatq(z|K) = Ÿ joeJ Cat(z j |r j )with 'j oe J , 'k oe K, rjk = fl jk q koeK fl jk .

  The red term E [log a k ] is deduced from properties of the Dirichlet distribution, the blue term log fl t jk is deduced from (4.37) and (4.38) in subsection 4.3.2 and the green term log fl qc jk has been detailed in Chapter 3. Hence, z is distributed a posteriori according to a product of Categorical distributions parametrized by r = (r jk ) (j,k)oeJ ◊K given by 'j oe J , 'k oe K, rjk = fl jk q koeK fl jk(4.49)

Figure 4 .

 4 Figure 4.21: Synthetic parabolic and piecewise parabolic data generated from di↵erent values of the variance parameter ‡ 2 . Figures (a), (b) and (c) present unlabeled data where 4 parabolas are generated. Ground-truth are visible on Figures (d), (e) and (f).

Figure 4 .

 4 Figure 4.22: Synthetic quantitative data generated from 4 multivariate normal distributions. Figure (a) shows unlabeled data and Figure (b) exhibits the ground-truth.

Figure 4 .

 4 10). Indeed the Silhouette Coe cient of the proposed model (S = 0.14) is closer to the Silhouette Coe cient of the ground-truth (S = 0.10) than Silhouette Coe cients of spectral clustering (S = 0.53) and k-means (S = 0.56).

Titre:

  Gestion des incertitudes en identification des modes radar Mots cl és : émetteurs radar ,classification,partitionnement,valeurs aberrantes,donn ées manquantes,mod èles de m élange R ésum é : En Guerre Electronique, l'identification des signaux radar est un atout majeur de la prise de d écisions tactiques li ées au th é âtre d'op érations militaires. En fournissant des informations sur la pr ésence de menaces, la classification et le partitionnement des signaux radar ont alors un r ôle crucial assurant un choix adapt é des contre-mesures d édi ées à ces menaces et permettant la d étection de signaux radar inconnus pour la mise à jour des bases de donn ées. Les syst èmes de Mesures de Soutien Electronique enregistrent la plupart du temps des m élanges de signaux radar provenant de diff érents émetteurs pr ésents dans l'environnement électromagn étique. Le signal radar, d écrit par un motif de modulations impulsionnelles, est alors souvent partiellement observ é du fait de mesures manquantes et aberrantes. Le processus d'identification se fonde sur l'analyse statistique des param ètres mesurables du signal radar qui le caract érisent tant quantitativement que qualitativement. De nombreuses approches m êlant des techniques de fusion de donn ées et d'apprentissage statistique ont ét é développ ées. Cependant, ces algorithmes ne peuvent pas g érer les donn ées manquantes et des m éthodes de substitution de donn ées sont requises afin d'utiliser ces derniers. L'objectif principal de cette th èse est alors de d éfinir un mod èle de classification et partitionnement int égrant la gestion des valeurs aberrantes et manquantes pr ésentes dans tout type de donn ées. Une approche fond ée sur les mod èles de m élange de lois de probabilit é est propos ée dans cette th èse. Les mod èles de m élange fournissent un formalisme math ématique flexible favorisant l'introduction de variables latentes permettant la gestion des donn ées aberrantes et la mod élisation des donn ées manquantes dans les probl èmes de classification et de partitionnement. L'apprentissage du mod èle ainsi que la classification et le partitionnement sont r éalis és dans un cadre d'inf érence bay ésienne o ù une m éthode d'approximation variationnelle est introduite afin d'estimer la loi jointe a posteriori des variables latentes et des param ètres. Des exp ériences sur diverses donn ées montrent que la m éthode propos ée fournit de meilleurs r ésultats que les algorithmes standards. Title : Uncertainty in radar emitter classification and clustering Keywords : radar emitter,classification,clustering,outliers,missing data,mixture models Abstract : In Electronic Warfare, radar signals identification is a supreme asset for decision making in military tactical situations. By providing information about the presence of threats, classification and clustering of radar signals have a significant role ensuring that countermeasures against enemies are wellchosen and enabling detection of unknown radar signals to update databases. Most of the time, Electronic Support Measures systems receive mixtures of signals from different radar emitters in the electromagnetic environment. Hence a radar signal, described by a pulse-to-pulse modulation pattern, is often partially observed due to missing measurements and measurement errors. The identification process relies on statistical analysis of basic measurable parameters of a radar signal which constitute both quantitative and qualitative data. Many general and practical approaches based on data fusion and machine learning have been developed and traditionally proceed to feature extraction, dimensionality reduction and classification or clustering. However, these algo-rithms can not handle missing data and imputation methods are required to generate data to use them. Hence, the main objective of this work is to define a classification/clustering framework that handles both outliers and missing values for any types of data.

  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Continuous radar features . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Realistic data acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Standard Gaussian mixture models . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Gaussian mixture models with missing data . . . . . . . . . . . . . . . . . 2.2.4 Gaussian mixture models with outliers . . . . . . . . . . . . . . . . . . . .

	2 Continuous data	21
	2.1

. . 1.2.2 Latent variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Bayesian framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.4 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.5 Classification and clustering . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 2.2.5 Proposed mixture model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Variational posterior distributions . . . . . . . . . . . . . . . . . . . . . . 2.3.2 VBE-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 VBM-step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii

Table 2 .

 2 1: Initialisation of hyper-parameters values for classification on continuous data

  Classification procedure on continuous data : Prediction step Input: Unlabelled dataset x pred and learned parameters ˜ train Output: Predicted labels zpred Update -jk , -jk , μx miss jk , ˜ x miss

	Procedure 2.2
	koeK
	end if
	end for

, rr , sk , qk , μk , ˜ k 2 k , rjk return zpred such that each zpred j = arg max koeK rjk

  Pulse-to-pulse modulations consist in modulating parameters of a group of pulses to minimize ambiguities related to range and velocity. They are mostly applied on RF and PRI parameters through deterministic and random patterns. These various patterns are defined below for parameter values (p i ) 1AEiAEn of a group of n pulses and are also visible in Figure3.2.

Constant modulation When there is no modulation, all values (p i ) 1AEiAEn are identical such that 'i oe {1, . . . , n}, p i = v where v is a constant. Slide modulation When parameter values (p i ) 1AEiAEn are sliding, they are linearly modulated around a nominal v value such that 'i oe {1, . . . , n}, p i

Table 3

 3

	.1: All known binary Barker codes
	Code length	Code
	2	11 or 10
	3	110
	4	1110 or 1101
	5	11101
	7	1110010
	11	11100010010
	13	1111100110101

rate of frequency change. That non linear variation results in spending more time at frequencies that need to be enhanced and in avoiding high autocorrelation sidelobes.

Table 3 .

 3 2: Categorical features observed in categorical data where 42 di↵erent combinations of these features are shared by the 55 emitters.

	Intrapulse	Pulse-to-Pulse PRI Pulse-to-pulse RF	Scan
	Barker 13	Dwell	Agility	Circular
	Barker 7	Triangular Dwell	Agility Burst	Sector
	Chirp	Complex	Diversity	None
	Diversity	HFR	FMICW
	Double Chirp Increasing Wobble	None
	FMCW	Jitter
	Trapeze Chirp	Sinus Jitter
	Phase Code	Stagger
	S Law	Wobble
	Trapeze	None
	None

Table 3 .

 3 3: Initialisation of hyper-parameters values for classification on mixed data

  According to the dataset visualised in Figure3.7, K min and K max are set to 35 and 85 in Procedure 3.4 Clustering procedure on mixed data Input: Unlabelled dataset x and a priori range of numbers of clusters K oe {K min , . . . , K max } Output: Labels z, parameters ˜ and optimal number of clusters Kfor K = K min to K max do Perform semi-supervised classification with K classes Stock labels zK , parameters ˜ K and L K

	end for return ˜ K and	z	K such that K = arg max K	L K
	Table 3.4: Initialisation of hyper-parameters values for clustering on mixed data
	zj = arg max koeK rjk			)oeK◊Cq	and z such that each
	end if
	end for

  jk is deduced from red, blue and green expectations in (4.13) as follows 'j oe J , 'k oe K, log fl jk = E [log a k ] + log fl t jk + log fl qc jk .

	(4.15)
	with
	'j oe J , 'k oe K,
	(4.14)
	the VBE (4.13) and VBM (4.14) steps can be independently derived for latent variables and parameters related to temporal evolution data x t and mixed data (x q , x c ). Therefore, variational posterior distributions of latent variables (x miss q , u, x miss c ) and parameters (fi, -, -, µ, ) related to mixed data (x q , x c ) are obtained as in Chapter 3 by deriving green expectations in (4.13) and (4.14). As for (Ê, ‡ 2 ), their variational posterior distribution are obtained as in subsection 4.1.2
	by developing blue expectations in (4.13) and (4.14). As in subsection 4.1.2 or in Chapter 3, the

respective parameters are estimated during the VBE and VBM steps by developing expectations E [log p(x, u, z|t, , K)] and E h [log p(x, u, z, |t, K)]. Noting that

E [log p(x, u, z|t, , K)] =E Ë log p(x t |t, z, Ê, ‡ 2 , K) È + E [log p(x q , u, x c |z,

fi, -, -, µ, , K)

]

+ E [log p(z|a, K)] , (4.13) and E h [log p(x, u, z, |t, K)] =E h Ë log p(x t , Ê, ‡ 2 |t, z, K) È + E h [log p(x q , u, x c ,

fi, -, -, µ

, |z, K)] + E h [log p(z, a|K)] , Dirichlet posterior distribution of a is deduced from the red expectation in (4.14). Eventually, the variational distribution of labels z is obtained by marginalising over latent variables in the green expectation and developing both blue and red expectations in (4.13) such that ⁄ E [log p(x, u, z|t, , K)] ˆxmiss q ˆuˆx miss c = ÿ joeJ ÿ koeK " k z j log fl jk where log fl

Table 4 .

 4 1: Initialisation of hyper-parameters values for clustering on parabolic data

	Ê 0	0	Ÿ 0

Table 4 .

 4 2: Adjusted Rand Index (ARI) and Silhouette coe cient (S) values for the proposed model (PM), the spectral clustering (SC) and the k-means algorithm (KM) during the first experiment on synthetic data when only temporal evolution data are considered.

		ARI	S
	‡ 2 = 0.0001 ‡ 2 = 0.01 ‡ 2 = 0.25	PM SC KM Data PM SC KM 1 1 0.34 0.32 0.32 0.32 0.63 1 1 0.33 0.27 0.27 0.27 0.63 0.42 0.23 0.27 0.10 0.14 0.53 0.56

Table 4 .

 4 3: Silhouette coe cients of the proposed model (PM), the spectral clustering (SC) and the kmeans algorithm (KM) during the first experiment on real data when only temporal evolution data are considered.

		Silhouette Coe cient
		PM SC	KM
		Case 1 0.05 0.46	0.57
		Case 2 0.12 0.28	0.57
		Case 3 0.26 0.61	0.61
	Table 4.4: Adjusted Rand Index (ARI) and Silhouette coe cient (S) values for the proposed model (PM),
	the spectral clustering (SC) and the k-means algorithm (KM) during the experiment on synthetic data
	when all types of data are considered.
		ARI	S
	‡ 2 = 0.0001 ‡ 2 = 0.01 ‡ 2 = 0.25	PM SC KM Data PM SC KM 1 1 1 0.76 0.76 0.76 0.76 1 1 1 0.75 0.75 0.75 0.75 1 1 1 0.73 0.73 0.73 0.73

Table 4 .

 4 5: Silhouette coe cients of the proposed model (PM), the spectral clustering (SC) and the k-means algorithm (KM) during the second experiment on real data when all types of data are considered.

	Silhouette Coe cient
	PM SC	KM
	Case 1 0.77 0.75	0.81
	Case 2 0.57 0.32	0.61
	Case 3 0.71 0.73	0.74

  x qj ) joeJ are normally distributed conditionally to categorical data (x cj ) joeJ . Both quantitative and categorical data (x qj , x cj ) joeJ can be partially observed. Hence (x qj , x cj ) joeJ are decomposed into observed features (x obs

	missing features (x miss qj , x miss cj ) joeJ such that		qj , x obs cj ) joeJ and
	'j oe J ,	x qj = x cj =	A x miss qj x obs qj A cj x obs cj x miss	B B	with (x miss qj , x obs qj ) oe R d miss j	j ◊ R d obs	and d miss j	+ d obs j	= d ,

  VBE (4.29) and VBM (4.30) steps can be independently derived for latent variables and parameters related to temporal evolution data x t and mixed data (x q , x c ). Therefore, variational posterior distributions of latent variables (x miss , x c ) are obtained as in Chapter 3 by deriving green expectations in (4.29) and (4.30). As for (y, b, µ t , ‡ 2 ), their variational posterior distribution are obtained as in subsection 4.2.2 by developing blue expectations in (4.29) and (4.30). As in subsection 4.2.2 or in Chapter 3, the Dirichlet posterior distribution of a is deduced from the red expectation in (4.30). Eventually, the variational distribution of labels z is obtained by marginalising over latent variables in the green expectation and developing both blue and red expectations in (4.29) such that

	.28) Their respective parameters are estimated during the VBE and VBM steps by developing expec-tations E [log p(x, u, y, z| , K)] and E h [log p(x, u, y, z, |K)]. Noting that E [log p(x, u, y, z| , K)] =E Ë log p(x t , y|z, b, µ t , ‡ 2 , K) È + E [log p(x q , u, x c |z, fi, -, -, µ, , K)] + E [log p(z|a, K)] , (4.29) and , u, x miss c ) and parameters (fi, -, -, µ, ) related to mixed data (x q ⁄ E [log p(x, u, y, z| , K)] ˆxmiss q ˆyˆuˆx miss c = ÿ joeJ ÿ koeK " k z j log fl jk where log fl jk is deduced from red, blue and green expectations in (4.29) as follows 'j oe J , 'k oe K, log fl jk = E [log a k ] + log fl t jk + log fl qc jk . (4.31) with E q 'j oe J , 'k oe K,

h [log p(x, u, y, z, |K)] =E h Ë log p(x t , y, b, µ t , ‡ 2 |z, K) È + E h [log p(x q , u, x c ,

fi, -, -, µ, |z, K)

]

+ E h [log p(z, a|K)] ,

(4.30) the

Table 4 .

 4 6: Initialisation of hyper-parameters values for clustering on piecewise parabolic data

Table 4

 4 , respectively covariance matrices 0 , are initialised from results of a k-means algorithm on quantitative data, respectively from diagonal matrices whose diagonal elements are variances of quantitative data. Other hyper-parameters are initialised as in Table4.6.

	.7: Adjusted Rand Index (ARI) and Silhouette coe cient (S) values for the proposed model (PM),
	the spectral clustering (SC) and the k-means algorithm (KM) during the first experiment on synthetic
	data when only temporal evolution data are considered.
	ARI	S
	PM SC KM Data PM SC KM ‡ 2 = 0.0001 0.73 0.73 0.35 0.30 0.35 0.35 0.64 ‡ 2 = 0.01 0.76 0.84 0.34 0.26 0.31 0.33 0.64 ‡ 2 = 0.25 0.54 0.26 0.27 0.10 0.14 0.53 0.56
	µ 0

  tj , t j ) have a parabolic relation . tj , t j ) have a piecewise parabolic relation . (4.33) Figure4.18 presents data where amplitudes x t and times of arrival t are distributed according to a parabolic relation and a piecewise parabolic relation from two distinct emitters. Figure 4.18: Data generated from two distinct emitters presenting a parabolic scanning behaviour and a piecewise parabolic scanning behaviour. Hence, amplitudes x t and times of arrival t are distributed according to a parabolic relation and a piecewise parabolic relation defined with P = 4 piecewises.

	÷p oe P, N 0.5 0.55 0.6 0.65 0.85 0.9 0.95 1 if (x 0 1 x tj |µ t p , ‡ 2 2 0.2 0.4 0.7 0.75 0.8 Amplitudes : x	0.6	0.8	1
	Times of Arrival : t

t

  .42) with z j the latent variable modeling the label of the j th observation vector x j = (x qj , x cj , x tj ). Moreover, the temporal evolution data (x tj , t j ) joeJ are distributed according to either a parabolic relation or a piecewise parabolic relation and the quantitative data (x qj ) joeJ are normally distributed conditionally to categorical data (x cj )

	x qj =	A x miss qj x obs qj	B	with (x miss qj , x obs qj ) oe R d miss j	◊ R d obs j	and d miss j	+ d obs j	= d ,
	x cj =

joeJ . Both quantitative and categorical data (x qj , x cj ) joeJ can be partially observed. Hence (x qj , x cj ) joeJ are decomposed into observed features (x obs qj , x obs cj ) joeJ and missing features (x miss qj , x miss cj ) joeJ such that 'j oe J , A x miss cj x obs cj B

  .45) Their respective parameters are estimated during the VBE and VBM steps by developing expec-‡ 2 ), their variational posterior distribution are obtained as in subsection 4.3.2 by developing blue expectations in (4.46) and (4.47). As in subsection 4.3.2 or in Chapter 3, the Dirichlet posterior distribution of a is deduced from the red expectation in (4.47). Eventually, the variational distribution of labels z is obtained by marginalising over latent variables in the green expectation and developing both blue and red expectations in (4.46) such that

						(4.46)
	+ E [log p(z|a, K)] ,
	and
	E h [log p(x, u, y, w, z, |t, K)] =E h					(4.47)
	+ E h [log p(z, a|K)] ,
	the VBE (4.46) and VBM (4.47) steps can be independently derived for latent variables and
	parameters related to temporal evolution data x q	ˆyˆwˆuˆx miss c	=	ÿ	ÿ	" k z
				joeJ	koeK

tations E [log p(x, u, y, w, z|t, , K)] and E h [log p(x, u, y, w, z, |t, K)]. Noting that

E [log p(x, u, y, w, z|t, , K)] =E Ë log p(x t , y, w|z, t, b, c, Ê, µ t , ‡ 2 , K) È + E [log p(x q , u, x c |z,

fi, -, -, µ, , K

)] Ë log p(x t , y, w, b, c, Ê, µ t , ‡ 2 |z, t, K) È + E h [log p(x q , u, x c ,

fi, -, -, µ, |z, K

)] t and mixed data (x q , x c ). Therefore, variational posterior distributions of latent variables (x miss q , u, x miss c ) and parameters (fi, -, -, µ, ) related to mixed data (x q , x c ) are obtained as in Chapter 3 by deriving green expectations in (4.46) and (4.47). As

for (y, w, b, c, Ê, µ

t , ⁄ E [log p(x, u, y, w, z|t, , K)] ˆxmiss j log fl jk where log fl jk is deduced from red, blue and green expectations in (4.46) as follows 'j oe J , 'k oe K, log fl jk = E [log a k ] + log fl t jk + log fl qc jk . (4.48) with

Table 4 .

 4 9: Initialisation of hyper-parameters values for clustering on parabolic and piecewise parabolic data

	Ê 0	0

Table 4 .

 4 10: Adjusted Rand Index (ARI) and Silhouette coe cient (S) values for the proposed model (PM), the spectral clustering (SC) and the k-means algorithm (KM) during the first experiment on synthetic data when only temporal evolution data are considered.

	ARI	S
	PM SC KM Data PM SC KM ‡ 2 = 0.0001 0.73 0.73 0.35 0.30 0.35 0.35 0.64 ‡ 2 = 0.01 0.76 0.84 0.34 0.26 0.31 0.33 0.64 ‡ 2 = 0.25 0.54 0.26 0.27 0.10 0.14 0.53 0.56

Table 4 .

 4 11: Adjusted Rand Index (ARI) and Silhouette coe cient (S) values for the proposed model (PM), the spectral clustering (SC) and the k-means algorithm (KM) during the first experiment on synthetic data when only temporal evolution data are considered. ARI S PM SC KM Data PM SC KM ‡ 2 = 0.0001 0.73 0.73 0.35 0.30 0.35 0.35 0.64 ‡ 2 = 0.01 0.76 0.84 0.34 0.26 0.31 0.33 0.64 ‡ 2 = 0.25 0.54 0.26 0.27 0.10 0.14 0.53 0.56

Acknowledgements

I would like to thank Technical Department of Thales Defense Mission Systems (DMS) and Inverse Problems Group (GPI) from Laboratory of Signals and Systems (L2S) for giving me the opportunity to realise a PhD.

Into the Technical Department, I mostly thank my advisors Cyrille Enderli and Jean-François Grandin for introducing me to Electronic Warfare domain and guiding me during these 3 years. I have really appreciated their motivation and their scientific open-mindedness.

Into the Inverse Problems Group, I would also like to truly thank my academic supervisors Ali Mohammad-Djafari and Charles Soussen for their concrete and wise advice, their full educational support and their time throughout the past 3 years.

As for corrections related to this manuscript, I am extremely grateful to my referees Nathalie Peyrard and Charles Bouveyron for their careful proofreading and their relevant remarks.

At last, I would like to lovely thank my beloved family and friends for their emotional support that has helped me to achieve this work. i

Performance on synthetic data are evaluated through the Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF] that compares estimated partitions of data with the ground-truth and the Silhouette Coe cient [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF] which does not require the ground-truth and provides a higher score when clusters are dense and well separated. Performance on real data are only evaluated through the Silhouette score since the ground-truth is not available for each case.

Experiments and results

The first experiment aims to determine the ability of each algorithm to restore the true clusters according to an a priori number of clusters K when only temporal evolution data are taken into consideration. According to datasets visualised in Figure 4.4 and Figure 4.6, K is set to 4 for synthetic data and to 5, 4 and 2 for the three real cases. Results of the first experiment on synthetic data are shown in Figure 4.7 and in Table 4.2. The proposed model and the spectral clustering succeed in clustering synthetic data for ‡ 2 oe {0.0001, 0.01} since the ground-truth partition is recovered in Figure 4.7 with an ARI equals to 1 visible on Table 4.2. The lower performance of

Publications

Journal Paper