Etienne Krönert
email: etienne.kronert@worldline.com

Alain Célisse
email: alain.celisse@univ-paris1.fr

Dalila Hattab
email: dalila.hattab@worldline.com

FDR control for Online Anomaly Detection

Keywords: Anomaly detection, Time series, Multiple testing

The goal of anomaly detection is to identify observations generated by a process that is different from a reference one. An accurate anomaly detector must ensure low false positive and false negative rates. However in the online context such a constraint remains highly challenging due to the usual lack of control of the False Discovery Rate (FDR). In particular the online framework makes it impossible to use classical multiple testing approaches such as the Benjamini-Hochberg (BH) procedure. Our strategy overcomes this difficulty by exploiting a local control of the "modified FDR" (mFDR). An important ingredient in this control is the cardinality of the calibration set used for computing empirical p-values, which turns out to be an influential parameter. It results a new strategy for tuning this parameter, which yields the desired FDR control over the whole time series. The statistical performance of this strategy is analyzed by theoretical guarantees and its practical behavior is assessed by simulation experiments which support our conclusions.

Introduction

Context

By observing indicators along the time to check the system health, anomaly detection aims at raising an alarm if abnormal patterns are detected [START_REF] Ahmed | A survey of network anomaly detection techniques[END_REF][START_REF] Li | An Intelligent Framework for Timely, Accurate, and Comprehensive Cloud Incident Detection[END_REF]. A motivation for automatic anomaly detection is to reduce the workload of operations teams by allowing them to prioritize their efforts where necessary. This is usually made possible by using statistical and machine learning models [START_REF] Chandola | Anomaly detection: A survey[END_REF][START_REF] Braei | Anomaly detection in univariate timeseries: A survey on the state-of-the-art[END_REF][START_REF] Blázquez-García | A review on outlier/anomaly detection in time series data[END_REF]. However when badly calibrated an anomaly detector leads to alarm fatigue. An overwhelming number of alarms desensitizes the people tasked responding to them, leading to missed or ignored alarms or delayed responses [START_REF] Cvach | Monitor alarm fatigue: an integrative review[END_REF][START_REF] Blum | Alarms in the intensive care unit: too much of a good thing is dangerous: is it time to add some intelligence to alarms?[END_REF]. One of the reasons for alarm fatigue is the high number of false positives which take time to be managed [START_REF] Solet | Managing alarm fatigue in cardiac care[END_REF][START_REF] Lewandowska | Determining Factors of Alarm Fatigue among Nurses in Intensive Care Units-A Polish Pilot Study[END_REF]. The main goal of the present work is to design a new (theoretically grounded) strategy allowing to control the number of false positives when performing automatic anomaly detection in sequential context.

Related works

Anomaly Detection in time series According to [START_REF] Hawkins | Identification of outliers 11[END_REF] an anomaly is "An observation which deviates so much from other observations as to arouse suspicions that it was generated by a different mechanism." Actually the features of what is usually called an "anomaly" is more diverse depending on the context. For example an anomaly can refer to a value that appears larger than the mean of other values, or to a sequence of values having a smaller variance than other neighboring sequences of observations. Anomaly detection in time series is applied in various domains such as fault [START_REF] Notaro | A survey of aiops methods for failure management[END_REF] or attack [START_REF] Pang | ANDEA: anomaly and novelty detection, explanation, and accommodation[END_REF] detection in information systems, fault detection in industrial equipment [START_REF] Li | MAD-GAN: Multivariate anomaly detection for time series data with generative adversarial networks[END_REF][START_REF] Miele | Deep anomaly detection in horizontal axis wind turbines using graph convolutional autoencoders for multivariate time series[END_REF], fault detection in vehicles [START_REF] Hundman | Detecting spacecraft anomalies using lstms and nonparametric dynamic thresholding[END_REF][START_REF] Han | ADS-lead: Lifelong anomaly detection in autonomous driving systems[END_REF], medical diagnostics [START_REF] Malhotra | Long Short Term Memory Networks for Anomaly Detection in Time Series[END_REF][START_REF] Li | A survey of heart anomaly detection using ambulatory electrocardiogram (ECG)[END_REF][START_REF] Wulsin | Semi-supervised anomaly detection for EEG waveforms using deep belief nets[END_REF] and astronomy [START_REF] Hönes | Automatically detecting anomalous exoplanet transits[END_REF]. The question of designing anomaly detectors has been considered in different contexts as explained by [START_REF] Chandola | Anomaly detection: A survey[END_REF]: supervised, unsupervised and semi-supervised. The present work focuses on the unsupervised framework which is the most widespread one since labeling is often difficult and costly. Sometimes, specific pattern of anomalies are known, and specific detectors can be designed to detect them [START_REF] Bos | Towards software-based signature detection for intrusion prevention on the network card[END_REF][START_REF] Hindy | A taxonomy and survey of intrusion detection system design techniques, network threats and datasets[END_REF][START_REF] Khraisat | Survey of intrusion detection systems: techniques, datasets and challenges[END_REF]. But most of the time the patterns have to be discovered as well. Review papers [START_REF] Chandola | Anomaly detection: A survey[END_REF][START_REF] Blázquez-García | A review on outlier/anomaly detection in time series data[END_REF][START_REF] Ruff | A unifying review of deep and shallow anomaly detection[END_REF] reveal the diversity of existing approaches. Some of the most important categories of anomaly detectors are: distance-based anomaly detectors identifying anomalies far from other past observations [START_REF] Li | Network anomaly detection based on TCM-KNN algorithm[END_REF], or probability-based anomaly detectors identifying anomalies within a low probability area of a statistical model fitted on past observations [START_REF] Ahmed | Online anomaly detection using KDE[END_REF][START_REF] Laxhammar | Anomaly detection in sea traffic-a comparison of the gaussian mixture model and the kernel density estimator[END_REF]. Prediction-based anomaly detector aims at building a forecast model on training data. Anomalies are defined with high prediction error on streaming data [START_REF] Roberts | Control chart tests based on geometric moving averages[END_REF], [START_REF] Pincombe | Anomaly detection in time series of graphs using arma processes[END_REF], [START_REF] Buda | Deepad: A generic framework based on deep learning for time series anomaly detection[END_REF], [START_REF] Munir | Deep-AnT: A deep learning approach for unsupervised anomaly detection in time series[END_REF]. And reconstructionbased anomaly detection uses dimension reduction and anomalies correspond to high reconstruction errors on streaming data [START_REF] Kundu | A3D: Attention-based auto-encoder anomaly detector for false data injection attacks[END_REF], [START_REF] Miele | Deep anomaly detection in horizontal axis wind turbines using graph convolutional autoencoders for multivariate time series[END_REF], [START_REF] Vahdat | NVAE: A deep hierarchical variational autoencoder[END_REF].

The present work describes a versatile anomaly detector that is sensitive to all categories of anomalies depending on the underlying score function (see Section 2.1).

A high diversity of abnormality score functions can be used to detect different abnormality patterns [START_REF] Foorthuis | On the nature and types of anomalies: a review of deviations in data[END_REF]. Abnormality scores are often not easily interpretable if the score distribution is unknown. Therefore, it is impossible to make a judicious choice of the detection threshold. The Conformal Anomaly Detection was introduced to alleviate this issue.

Conformal Anomaly Detection Conformal Anomaly Detection (CAD) introduced in [START_REF] Laxhammar | Conformal anomaly detection: Detecting ab-normal trajectories in surveillance applications[END_REF] is a method derived from Conformal Prediction [START_REF] Angelopoulos | A gentle introduction to conformal prediction and distribution-free uncertainty quantification[END_REF]. The goal of CAD is to give a probabilistic interpretation of the score using empirical p-values. Inductive Conformal Anomaly Detection (ICAD) introduced in [START_REF] Laxhammar | Conformal anomaly detection: Detecting ab-normal trajectories in surveillance applications[END_REF] improves the CAD linear complexity in time and adapts it for Online Anomaly Detection by introducing the concept of calibration set. CAD can be used with a wide variety of anomaly score functions. For instance [START_REF] Safin | Conformal kernel expected similarity for anomaly detection in time-series data[END_REF] presents an anomaly detector based kernels combined with CAD. [START_REF] Burnaev | Conformalized density-and distance-based anomaly detection in time-series data[END_REF] combined distance and density based scoring function with CAD. CAD gives the opportunity to control the expected number of false positives within a time period. But its main limitation is that it yields no control over the false alarm rate on the whole time series that is, the proportion of false positive among all detections. By contrast the present work aims at having a control over it (to avoid alarm fatigue [START_REF] Blum | Alarms in the intensive care unit: too much of a good thing is dangerous: is it time to add some intelligence to alarms?[END_REF][START_REF] Cvach | Monitor alarm fatigue: an integrative review[END_REF]), more precisely on the False Discovery Rate (FDR).

FDR Control Benjamini-Hochberg (BH) procedure [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF][START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] is a multiple testing procedure that controls the proportions of false positives among rejections that is the False Discovery Rate. The BH procedure can be improved by estimating the proportion of anomalies in the dataset [START_REF] Robin | A cross-validation based estimation of the proportion of true null hypotheses[END_REF][START_REF] Guedj | Kerfdr: a semiparametric kernel-based approach to local false discovery rate estimation[END_REF]. Most procedure based on Benjamini-Hochberg assume that the true p-values are known. When the distribution of the scores under the null hypothesis is unknown it is generally not possible to ensure the FDR control with BH. For instance the Monte-Carlo Multiple-Testing has been suggested by [START_REF] Guo | Adaptive choice of the number of bootstrap samples in large scale multiple testing[END_REF][START_REF] Gandy | MMCTest-a safe algorithm for implementing multiple Monte Carlo tests[END_REF][START_REF] Zhang | Adaptive monte carlo multiple testing via multi-armed bandits[END_REF] to overcome this difficulty. An alternative method for controlling FDR is based on the so-called "local FDR" [START_REF] Sun | Oracle and adaptive compound decision rules for false discovery rate control[END_REF][START_REF] Wang | Online fdr controlled anomaly detection for streaming time series[END_REF]. Unfortunately this approach relies on a Gaussian assumption. In the context of ICAD, FDR can be controlled using conformal p-values with BH as shown in [START_REF] Bates | Testing for outliers with conformal p-values[END_REF][START_REF] Weinstein | A power and prediction analysis for knockoffs with lasso statistics[END_REF][START_REF] Mary | Semi-supervised multiple testing[END_REF]. Moreover the FDR control can be achieved simultaneously with upper and lower bound as suggested in [START_REF] Marandon | Machine learning meets false discovery rate[END_REF].

Online FDR Control In online multiple-testing, the decision of new observed value as an anomaly has to be done instantaneously. If the BH procedure is applied on the current time series, the time complexity will increase with the size of the data. To tackle this problem, recent papers advocate different methods for the online control of the FDR [START_REF] Wang | Online fdr controlled anomaly detection for streaming time series[END_REF][START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF][START_REF] Ramdas | SAF-FRON: an adaptive algorithm for online control of the false discovery rate[END_REF][START_REF] Xu | Dynamic algorithms for online multiple testing[END_REF]. In [START_REF] Wang | Online fdr controlled anomaly detection for streaming time series[END_REF] the author suggests using the principle of local FDR. At each observation, a decision is taken depending on the estimation of the local FDR. The [START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF][START_REF] Ramdas | SAF-FRON: an adaptive algorithm for online control of the false discovery rate[END_REF][START_REF] Xu | Dynamic algorithms for online multiple testing[END_REF] introduce a method based on alpha-investing. The p-value is compared to an adaptive threshold depending on the previous decisions. But this method is not applicable for conformal p-values because of its low detection power.

Controlling false positives for online anomaly detection remains a difficult task. In particular two challenges arise with online anomaly detection:

• The true p-values are unknown and need to be estimated.

• The decisions are made in an online context, whereas most of the multiple testing methods are done in the offline context.

The main contributions are to tackle these challenges. More precisely it is established that it is possible to design online anomaly detectors controlling the FDR of the time series.

• This paper study the relationship between the FDR and the cardinality of the calibration set used to estimate p-values. To guarantee FDR control, a calibration set cardinality tuning method is proposed.

• This paper describes an online calibration strategy for anomaly detection based on multiple testing ideas to control the False Discovery Rate (FDR).

-It explains how global control of the time series FDR can be obtained from local control of a modified version of subseries FDR. This makes it possible to control the FDR within an online context.

-A modified version of the Benjamini-Hochberg procedure is suggested to achieve local control of the modified FDR.

Description of the paper

First, the problem is explained and important objects are introduced in Section 2. Second Section 3 deals with conditions on p-values estimations to ensure local control of FDR is controlled at a desired level. Third this paper develops algorithms that allows global control the FDR time series and studies them in Section 4. Finally our solution is evaluated against one competitor from the literature in Section 5.

Statistical framework

The Anomaly Detector

Let (X , Ω, P) be a probability space and assume a realization of the random variables (X t ) t≥1 , with X t taking values in X for all t, is observed at equal time steps. T ∈ N ∪ {∞} is the size of the time series. Let P 0 be a probability distribution, called reference distribution, on the space X . For each instant t, the observation X t is called "normal" if X t ∼ P 0 . Otherwise, X t is an "anomaly". The aim of an online anomaly detector is to find all anomalies among the new observations along the time series (X t ) t≥1 : for each instant t > 1, a decision is taken about the status of X t based on past observations: (X s ) 1≤s≤t .

In this paper, the following general online anomaly detector description is suggested. It uses multiple testing ideas from [START_REF] Mary | Semi-supervised multiple testing[END_REF] and the online context from [START_REF] Laxhammar | Conformal anomaly detection: Detecting ab-normal trajectories in surveillance applications[END_REF]. It is based on the three following notions:

1. Atypicity score: A score a : X → R is a function reflecting the atypicity of an observation X t . To be more specific, the further P Xt from P 0 , the larger a(X t ).

p-value:

It is the probability of observing a(X) higher than a(X t ) when X ∼ P 0 : p t = P X (a(X) ≥ a(X t )) .

The p-value enables an interpretable criterion measuring how much unlikely X t ∼ P 0 . It is estimated using empirical p-value, given in Definition 2, and it allows to tackle the problem of low interpretability of an unknown distribution of the atypicity score.

3. Detection threshold: ε ∈ [0, 1], it discriminates observations considered as abnormal from others. The observations considered as anomalies are x t whose (estimated) p-value is lower than the threshold ε, P X (a(X) ≥ a(X t )) < ε.

This general online anomaly detector is formalized in the pseudo-code given by Algorithm 1.

Algorithm 1 Generic Online Anomaly Detector with fixed threshold ε

Require: T > 0, (Xt) T 1 time series, a an abnormality score, ε a threshold

for 1 ≤ t ≤ T do St ← a(Xt) pt ← p-value(St) if pt < ε then dt = 1 else dt = 0 end if end forOutput: (dt) T 1 detected

anomalies boolean

The threshold ε allows to control the "detection frequency": a smaller threshold will generate fewer detections. This is equivalent to defining anomalies as points above a quantile in the tail of the score distribution. Nevertheless, in practice, the calibration of the threshold is difficult. Since ε affects directly the number of false detections (false positives), it is not possible to know in advance the number of false positives due to the choice of ε.

One of the main contributions of this work consists in developing a datadriven rule allowing to choose a threshold εt at each time step t. This rule has the advantage of ensuring a global control of the false discovery rate on the complete set of observations. In this paper, it is suggested to replace the last step of the Generic Online Anomaly Detector (Algorithm 1) with an adaptive threshold that is computed in real time before any decision would be made.

Control of false positives and multiple testing

Since the present goal is to use FDR, a natural strategy is to rephrase the online anomaly detection problem as a multiple testing problem: At each step 1 ≤ t ≤ T , a statistical test is performed on the hypotheses: H 0t , "X t is not an anomaly " against H 1t "X t is an anomaly".

A natural criterion controlling the proportion of type I errors (False Positives) of the whole time series is FDR [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. For a given data-driven threshold ε and a set of estimated p-values p = (p t ) t≥1 , the FDR criterion of the sequence from 1 to T is given by

F DR T 1 (ε, p) = E[F DP T 1 (ε, p)],
with

F DP T 1 (ε, p) = t∈H0 1[p t ≤ εt ] T t=1 1[p t ≤ εt ]
,

with the convention that 0/0 = 0. In the above expression, F DP T 1 denotes the False Discovery Proportion (FDP) of the time series from 1 to T . Also

H 0 = {t ∈ N * |H 0t
is true} is called the set of null hypotheses. Let us emphasize that the anomalies (according to Algorithm 1) satisfy 1[p t ≤ εt ] = 1. The main objective of the present work is to define a data-driven sequence ε : t → εt such that, for a given control level α ∈ [0, 1], under weak assumptions on the sequence p : t → pt ,

F DR T 1 (ε, p) ≤ α (2.1)
The control is said exact when "≤" is replaced with "=". Such a control would imply that for a level α = 0.1, at most 10% of the detected anomalies along the whole time series are false positives.

The detection power of the anomaly detector is measured by means of the False Negative Rate defined, for the sequence from 1 to T , by

F N R T 1 (ε, p) = E[F N P T 1 (ε, p)], (2.2)
with

F N P T 1 (ε, p) = t∈H1 1[p t ≤ εt ] |H 1 | , ( 2.3)
where F N P T 1 denotes the False Negative Proportion (FNP) of the sequence from 1 to T and H 1 = {t ∈ N * |H 1t is true} is the set of alternative hypotheses.

However a crucial remark at this stage is that controlling FDR on the complete time series is a highly challenging task in the present online context for at least two reasons:

• The main existing approaches for controlling FDR are described in an "offline" framework where the whole series is observed first, and decisions are taken afterwards [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF][START_REF] Marandon | Machine learning meets false discovery rate[END_REF]. This makes these approaches useless in the present context.

• The already existing approaches designed in the online context [START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF][START_REF] Xu | Dynamic algorithms for online multiple testing[END_REF] are difficult to parameterize and hard to apply with estimated p-values. Let us emphasize that realistic scenarios usually exclude the knowledge of the true probability distribution of the test statistics, leading to approximating or estimating the related p-values in practice.

FDR control with Empirical p-value

A classical (offline) strategy for controlling FDR is the so-called Benjamini-Hochberg (BH) multiple testing procedure [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. Exact control relies on the knowledge of true p-values, which is usually not realistic. Actually since the true reference distribution is unknown in practical anomaly detection scenarios, there is no true p-values available.

Empirical p-value

The atypicity level of an observation is quantified by an atypicity score. The underlying scoring function assigns each observation with a real value such that the more atypical the observation, the higher the score value.

Definition 1 (Scoring function)

. The scoring function is introduced: a : X → R.

The abnormality score u at the point x is defined by u = a(x).

The interpretation is that the higher the score (value) at an observation, the more unlikely the corresponding observation has been generated from a reference distribution implicitly encoded in the scoring function.

Examples (Examples of scoring functions). Let x ∈ X and z ℓ 1 = {z 1 , . . . , z ℓ } be a training set generated from P 0 .

1. Z-score [START_REF] Zhou | Continuous anomaly detection in satellite image time series based on z-scores of season-trend model residuals[END_REF][START_REF] Buda | Deepad: A generic framework based on deep learning for time series anomaly detection[END_REF]: Let µ and σ be estimators of mean and standard deviation of

z ℓ 1 , a(x) = a Z (x, z ℓ 1 ) = |(x -µ)/σ|
2. kNN score [START_REF] Smith | Anomaly detection of trajectories with kernel density estimation by conformal prediction[END_REF][START_REF] Burnaev | Conformalized density-and distance-based anomaly detection in time-series data[END_REF]: Let d be a metric on X and k > 0 and kN N (x, z ℓ 1 ) is the k-th nearest neighbor of x in z ℓ 1 .

a(x) = a kN N (x; z ℓ 1 , k) = 1 k z∈kN N (x,z ℓ 1 )

d(x, z)

3. KDE-score [START_REF] Smith | Anomaly detection of trajectories with kernel density estimation by conformal prediction[END_REF][START_REF] Burnaev | Conformalized density-and distance-based anomaly detection in time-series data[END_REF]: let K a kernel function.

a(x) = a KDE (x, K, z ℓ 1 ) = 1 ℓ z∈z ℓ 1 (-K(x, z))
4. Auto-Encoder score [START_REF] Zhou | Continuous anomaly detection in satellite image time series based on z-scores of season-trend model residuals[END_REF][START_REF] Miele | Deep anomaly detection in horizontal axis wind turbines using graph convolutional autoencoders for multivariate time series[END_REF]. Considers having

X = R n , a compression function f θ : R n → R d with d < n, a reconstruction function g ϕ : R d → R n
and the hypothesis spaces Θ and Φ.

( θ, φ) ∈ arg min

θ∈Θ,ϕ∈Φ x∈z ℓ 1 ||g ϕ (f θ (x)) -x|| 2
Then the abnormality score is given by

∀x ∈ R n , a(x) = a AE (x, z ℓ 1 ) = ||g φ(f θ (x)) -x|| 2 .
The choice of the abnormality score depends on the structure of the time series and the type of anomalies one is looking for. Intuitively a desirable scoring function should assign a high abnormality score to any true anomaly. For example, Z-score is only able to detect anomalies that are in the tail of the distribution. By contrast it is not effective to detect abnormal point between two modes of data with a bimodal distribution [START_REF] Zhou | Continuous anomaly detection in satellite image time series based on z-scores of season-trend model residuals[END_REF][START_REF] Buda | Deepad: A generic framework based on deep learning for time series anomaly detection[END_REF]. kNN or KDE scores are more suited for multi-modal data because they raise a high score for points far from the observations of the training set. The intuition behind such a scoring function is that normal data should have a low distance from the training set. Auto-Encoders are often used for multidimensional data with complex distributions. It relies on the possibility to compress the input data in a low dimensional embedding space without losing too much information. This enables to reconstruct the input with low reconstruction error, at the computational price of training first a deep neural network [START_REF] Zhou | Continuous anomaly detection in satellite image time series based on z-scores of season-trend model residuals[END_REF][START_REF] Miele | Deep anomaly detection in horizontal axis wind turbines using graph convolutional autoencoders for multivariate time series[END_REF]. To the best of our knowledge, there does not exist any scoring function suitable for detecting all types of anomalies.

Defining a meaningful threshold from a score is the classical strategy for deciding that an observation is anomalous or not. This requires to know the true distribution of these scores, which is not realistic in general. The induced estimation step is usually made by two means. On the one hand, one can assume a parametric Gaussian distribution for the scores [START_REF] Son | Anomaly detection for big log data using a Hadoop ecosystem[END_REF][START_REF] Buda | Deepad: A generic framework based on deep learning for time series anomaly detection[END_REF][START_REF] Hu | Anomaly detection using local kernel density estimation and context-based regression[END_REF]. On the other hand, one can estimate the score distribution by use of sampling techniques [START_REF] Safin | Conformal kernel expected similarity for anomaly detection in time-series data[END_REF][START_REF] Burnaev | Conformalized density-and distance-based anomaly detection in time-series data[END_REF]. Since the Gaussian assumption can cause some troubles when it is violated, the present work rather focuses on the second strategy by considering anomaly detection relying on empirical p-values. By contrast to the Gaussian assumption, a strong asset of empirical p-value is that they can be used no matter the true score distribution or the scoring function.

Definition 2 (Empirical p-value). Let a be a scoring function (Definition 1).

Let {x 1 , ..., x n } ⊂ X be a set of data called the calibration set. The empirical p-value is a function defined by

∀x ∈ X , p-value(x; {x 1 , ..., x n }) = 1 m m i=1 1(a(x i ) ≥ a(x)).
(2.4)

Let us emphasize that Definition 2 describes an estimator of the p-value under P 0 provided the calibration set is composed of points generated from the reference distribution P 0 . However it is well known that the main difficulty with this p-value estimator is that it is not itself a p-value [START_REF] Mary | Semi-supervised multiple testing[END_REF] since the so-called super-uniformity property is violated. More precisely, super-uniformity means that, for all u ∈ [0, 1],

P X,X1,...,Xm∼P0 (p-value(X; {X 1 , . . . , X m }) ≤ u) ≤ u.
Therefore empirical p-values are usually replaced by an other p-value estimator called the conformal p-value [START_REF] Marandon | Machine learning meets false discovery rate[END_REF][START_REF] Laxhammar | Conformal anomaly detection: Detecting ab-normal trajectories in surveillance applications[END_REF], given by p-value(x; {x 1 , ..., As a consequence, an important remark is that the present work focuses on empirical p-values (and not on conformal ones). However another motivation for this choice is provided in Section 3.2.2 where it is proved that the superuniformity property also holds true with empirical p-values under some specific conditions that will be detailed later.

x n }) = 1 m + 1 1 + m i=1 1(a(x i ) ≥ a(x)) .

BH-procedure does not control FDR with empirical p-values

The present section starts by describing the behavior of the BH-procedure as well as establishing the resulting FDR control. An illustration is provided that the BH-procedure does not control FDR are the prescribed level when empirical pvalues are used. This illustration is then theoretically justified, which shows that straightforwardly using the BH-procedure in our online context is prohibited. Definition 3 [START_REF] Weinstein | A power and prediction analysis for knockoffs with lasso statistics[END_REF])). Let m be an integer and α ∈ [0, 1]. Let (p i ) 1≤i≤m ∈ [0, 1] m be a family of p-values. The Benjamini-Hochberg (BH) procedure, denoted by BH α , is given by • a data-driven threshold:

ε BHα = max{ αk m ; p (k) ≤ αk m , k ∈ 1, m },
• a set of rejected hypotheses:

BH α ((p i ; i ∈ 1, m )) = {i; p i ≤ ε BHα , i ∈ 1, m } .
The intuition behind this procedure consists in drawing the ordered statistics i → p (i) (Figure 1) with p (1) ≤ • • • ≤ p (n) and the straight line i → αi m . Then the BH-procedure amounts to rejecting all hypotheses corresponding to p-values smaller than the last crossing point between the straight line and the ordered p-values curve.

Figure 1: Illustration of the Benjamini-Hochberg procedure. p-values are sorted by increasing order. The threshold is the greatest p-value that is lower than αk/m, when k is the rank of the p-value.

The striking property of this procedure is to yield the desired control of the FDR at the prescribed level α as stated by the next result.

Theorem 1 (FDR control with BH [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]). Let m be a positive integer and (X i ) m 1 be independent random variables such that X i ∼ P 0 , 1 ≤ i ≤ m 0 , and X i ∼ P 1 , m 0 + 1 ≤ i ≤ m. Let us also define the set of true p-values, for all 1 ≤ i ≤ m by p i = P X∼P0 (X ≥ X i ) ∈ [0, 1], and assume that each p i ∼ U ([0, 1]). Then for every α ∈]0, 1], BH α applied to p = (p i ) 1≤i≤m yields the exact FDR control at the prescribed level α that is,

F DR m 1 (ε BHα , p) = m 0 α m .
The proof of the theorem is deferred to Appendix A.3. In particular, the FDR control results from the fact that under H 0 , the true p-values follow a uniform distribution. The equality could be replaced by an upper bound if the uniform distribution assumption were weakened by the super-uniform property.

By contrast with the previous framework, when performing anomaly detection, the abnormality score is computed using a scoring function (Definition 1), and the true p-value is given by

p t = P X∼P0 (a(X) ≥ a(X t )) ,
where the notation clearly emphasizes the dependence with respect to the unknown reference distribution. This justifies why empirical p-values are now substituted to true ones as earlier explained (see Eq. 2.4).

A difficulty resulting from using empirical p-values in the BH-procedure is that the FDR control does no longer hold true as illustrated by Figure 2. This figure displays the actual FDR value (plain blue curve) versus the cardinality of the calibration set used to computed the empirical p-values (see Definition 2) in the specific situation of Gaussian data. Except for some a few values of the calibration set cardinality, the FDR control is no longer achieved (red horizontal line). Furthermore the actual FDR value is higher than the desired m 0 /mα. This results from the fact that the super uniform property is violated when using empirical p-values as established by Proposition 1 below.

Proposition 1 (Distribution of empirical p-value under H 0 ). Let X ∼ P 0 where P 0 is the probability distribution under H 0 , the calibration set cardinality is denoted by n, and {X 1 , . . . , X n } ∼ P n 0 is the calibration set. If one further assumes that there are no ties among a(X 1 ), . . . , a(X n ), then the empirical p-value at X is denoted by p-value(X; {X 1 , . . . , X n }) and follows the discrete uniform distribution

U (0, 1 n , 2 n , . . . , 1).
Let us mention that under H 0 the empirical p-value has a different distribution from that of the conformal p-value [START_REF] Marandon | Machine learning meets false discovery rate[END_REF] which follows U (1/(n + 1), . . . , 1). The conformal p-value is never smaller than 1/(n + 1), which raises issues in terms of the detection power with lots of false negatives (see the right panel of Figure 11 and the discussion in Appendix A.1). With empirical p-values, it can be easily checked that

P (p-value(X; {X 1 , . . . , X n }) ≤ 0) = 1 n + 1 > 0,
which violates the super uniformity property. As a consequence, FDR is no longer controlled by the BH-procedure [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] applied to empirical p-values. Other consequences owing to the use of empirical p-values violating super uniformity are illustrated in Section 3.1.

The assumption of no ties are allowed among the scores a(X i )s is quite mild and fulfilled most of the time as supported by Example 2.3.1 as long as the reference distribution is continuous (admits a density).

Proof of Proposition 1. Since X ∼ P 0 and X 1 , . . . , X n ∼ P n 0 are independent, the empirical p-value p is a random variable computed from X and X n 1 and satisfies for any 0 ≤ ℓ ≤ n that

P(p = ℓ/n) = P(a(X (ℓ+1) ) < a(X) ≤ a(X (ℓ) )) = P({rank of a(X) among {a(X), a(X 1 ), . . . , a(X n )} is ℓ + 1}), where a(X (n) ) < a(X (n-1) ) < • • • < a(X (1)
). The conclusion comes from noticing that the probability distribution of {X, X 1 , . . . , X n } is exchangeable, and assumption of no ties in scores a(X i ) Let us also notice that Figure 2 shows that there exist particular values of the calibration set cardinality for which FDR is still controlled at the prescribed level. This perspective is further explored in Section 3.2.2, where a new multiple testing procedure yielding the desired FDR control for the whole time series is devised.

FDR control with Empirical p-values

The goal here is to describe a strategy achieving the desired FDR control for a time series of length m when using empirical p-values. A motivating example is first introduced for emphasizing the issue in Section 3.1. Then a theoretical understanding is provided along Section 3.2 which results in a new solution which applies to independent empirical p-values. An extension is then discussed to thz non-independent setup in Section 3.3.1. Finally experimental results are reported in Section 3.4 to (empirically) assess the validity of our previous theoretical conclusions.

Motivating example

The purpose here is to further explore the effect of the calibration set cardinality on the actual FDR control when using empirical p-values. This gives us more insight on how to find mathematical solutions.

Let us start by generating observations using two distributions. The reference distribution is P 0 = N (0, 1) and the alternative distribution is P 1 = N (4, 10 -4 ). The anomalies are located in the right tail of the reference distribution. The length m of the signal is m = 100. The number of observations under P 0 is m 0 = 99. The experiments have been repeated K = 10 4 times.

Figure 2 displays the actual value of FDR as a function of the cardinality n of the calibration set {x 1 , ..., x n } used to compute the empirical p-values (see Definition 2). One clearly see that FDR is not uniformly controlled at level m 0 /mα. However there exist particular values of n for which this level of control is nevertheless achieved. As long as n has become large enough (n ≥ 500), repeated picks can be observed with a decreasing height as n grows.

FDR control: main results for i.i.d. p-values

The present section aims at first explaining the shape of the curve displayed in Figure 2. This will help getting some intuition about how to design an online procedure achieving the desired FDR control for the full time series.

Proof of FDR control by BH revisited

The main focus is now given to independent p-values. In what follows, the classical proof (Proof A.3) of the FDR control by the BH-procedure is revisited then leading to the next result. Its main merit is to provide the mathematical expression of the plain blue curve observed in Figure 2.

Theorem 2. Let n be the cardinality of the calibration set and m be that of the set of tested hypotheses where {X 1 , . . . , X m } denotes a set of random variables. Let m 0 ≤ m be the cardinality of the random variables from the reference distribution P 0 . Let the empirical p-value be denoted, for any i ∈ 1, m , by pi = p-value(X i , {Z i,1 , . . . , Z i,n }), where the calibration set is {Z i,1 , . . . , Z i,n } and each Z i,j ∼ P 0 . Let the random variables R(i) be the number of detections raised by BH α when replacing X i with 0, as defined along the proof detailed in Appendix A. [START_REF] Angelopoulos | A gentle introduction to conformal prediction and distribution-free uncertainty quantification[END_REF]. Then for every α ∈]0, 1], the FDR value over the sequence from 1 to m is given by

F DR m 1 (ε BHα , p) = m 0 m k=1 ⌊ αkn m ⌋+1 n+1 k P(R(i) = k),
where εBHα denotes the BH α threshold from Definition 3 when the BH-procedure is applied to the empirical p-values p = (p i ; 1 ≤ i ≤ m).

In general it is not possible to compute the exact value of the F DR without knowing the distribution of the random variables R(i). This is in contrast with the case of true p-values where P(p i ≤ αk m ) = αk m , whereas with empirical p-

values P(p i ≤ αk m ) = ⌊ αkn m ⌋+1
n+1 , which prevents from any simplification of the final bound. Nevertheless this value still suggests a solution to circumvent this difficulty: requiring conditions on α, m, and n such that

⌊ αkn m ⌋+1 n+1
= αk m , for all k. This is precisely the purpose of next Corollary 1.

Proof of Theorem 2. When applying Proof A.3, the only modification is that P(p i ≤ αk m ) is not equal to αk m since pi now follows the discrete uniform distribution

P(p i ≤ αk m ) = ⌊nαk/m⌋ ℓ=0 P(np i = ℓ) = ⌊ αkn m ⌋ + 1 n + 1 .
Plugging this in the FDR expression, it gives

F DR m 1 (ε BHα , p) = m 0 m k=1 ⌊ αkn m ⌋+1 n+1 k P(R(i) = k).
Recall that pi follows U (0, 1/n, 2/n, . . . , 1) entails that np i follows U (0, 1, 2, . . . , n).

Tuning of the calibration set cardinality

Corollary 1. Under the same notations and assumptions as Theorem 2, the next two results hold true.

Assume that there exists an integer 1 ≤ ℓ such that ℓm α is an integer. If the cardinality n of the calibration set satisfies

n = n ℓ -1 = ℓm/α -1, then F DR m 1 (ε BHα , p) = m 0 α m .
2. For every α ∈]0, 1], assume that the cardinality of the calibration set sat-

isfies n = n ℓ -1 = ⌈ ℓm α ⌉ -1, for any integer ℓ ≥ 1. Then, n (n + 1) m 0 α m ≤ F DR m 1 (ε BHα , p) ≤ m 0 α m .
The proof is postponed to Section A.4. The first statement in Corollary 1 establishes that recovering the desired control of FDR at the exact prescribed level α is possible on condition that the calibration set cardinality is large enough and more precisely that n = ℓm/α -1. This (mild) restriction on the values of α reflects that the empirical p-values do not satisfy the super-uniformity property. By contrast, the second statement yields the desired control at the level αm 0 /m by means of lower and upper bounds. In particular, the lower bound tells us that the FDR value can be not lower than the desired level αm 0 /m up to a multiplicative factor equal to 1-1/n, which goes 1 as n grows. For instance with α = 0.1 and m = 100, n ℓ = 1000 would yield that

F DR m 1 (ε BHα , p) • m/(m 0 α) ∈ [0.999, 1]
. This small lack of control is the price to pay for allowing any value of α ∈]0, 1]. It is also important to recall that in the anomaly detection field, abnormal events are expected to be rare. As a consequence m0 m is close to 1 and the actual F DR level is close to the desired α. However in situations where m 0 /m could depart from 1 too strongly, then incorporating an estimator of m 0 /m would be helpful.

Extension to dependent p-values

In Section 3.2, Corollary 1 states that FDR is controlled at a prescribed level with empirical p-values for which the super-uniformity property is not fulfilled. A key ingredient in the proof was the independence property across empirical p-values. One purpose of the present section is to extend these results to nonindependent p-values.

Towards this extension, the concept of positive regression dependency (referred to as PRDS) [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] turns out to be useful. The PRDS property is a form of positive dependence between p-values where all pairwise p-value correlations are positive. It results that a small p-value for a given observation makes other p-values for all considered observations simultaneously small as well, and viceversa [START_REF] Bates | Testing for outliers with conformal p-values[END_REF].

Theoretical results to dependent p-values

A classical result established in [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] proves that FDR is upper bounded by αm 0 /m provided the p-value family satisfies the PRDS and super-uniformity properties. It turns out that this result can be extended to our estimator with the same choice of calibration set cardinality as the one discussed in Corollary 1. Another important achievement is the fact that FDR can be also lower bounded in the case where the calibration set is the same for all (empirical) p-values (see Definition 2). This results originally proved by [START_REF] Marandon | Machine learning meets false discovery rate[END_REF] is extended here to empirical p-values computed with a calibration set cardinality tuned as suggested in Corollary 1.

In Section 3.2, considering the size of the calibration set is correctly chosen, it has been proved that the control of the FDR can be achieved with estimated p-values for which the super-uniformity property is not fulfilled. The results obtained for i.I.d. p-values will be extended in this section for non i.i.d. p-values.

For this extension, the concept of positive regression dependency on each one from a subset called PRDS [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] is introduced. The PRDS property is a form of positive dependence of p-values where all pairwise p-value correlations are positive. Larger scores in the calibration set make the p-values for all test points simultaneously smaller, and vice-versa [START_REF] Bates | Testing for outliers with conformal p-values[END_REF].

Definition 1 (PRDS property). A family of p-values pm

1 is PRDS on a set I 0 ⊂ {1, . . . , m} if for any i ∈ I 0 and any increasing set A, the probability

P[p m 1 ∈ A|p i = u] is increasing in u.
A classical result in [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] asserts that the FDR is upper bounded by m0 m α in the case where the p-value family is PRDS and super-uniform. This result can be extended our estimator with the same choice of calibration set cardinality than in Theorem 1.

Corollary 2 (Corollary of Theorem 1.2 in [6]). Suppose the family of p-values pm

1 is PRDS on the set H 0 of true null hypotheses and suppose that pm 1 respects super-uniformity an all thresholds that can may resulting from BH

∀k ∈ 1, m P(p i < αk m ) ≤ αk m , Then, the F DR is upper-bounded by α F DR(p, εB H) ≤ m 0 α m
Unique calibration set More over, in the case where the calibration set is the same for all p-values, the FDR can also be lower bounded as shown in [START_REF] Marandon | Machine learning meets false discovery rate[END_REF]. This result can be extended to empirical p-values given in Definition 2 with a calibration set cardinality tuned as proposed in Theorem 1.

Corollary 3 (Corollary of Theorem 3.4 in [START_REF] Marandon | Machine learning meets false discovery rate[END_REF]). Assuming the following conditions: Let n be the cardinality of the calibration set, m be the cardinality of the active set and m 0 the number of normal observations. Let P 0 be the reference distribution. Let Z i for i in 1, m independents random variables, following P 0 . Let X i for i in 1, m be random independents variables and independents from (Z j ). There are exactly m 0 random variables following the P 0 distribution. Let a be a scoring function. For all i in 1, m , let pi be the empirical p-values associated with the random variables X i and computed as follows,

pi = p-value(X i , {Z 1 , ..., Z m }, a).

If the cardinality of the calibration set is a multiple of

n = n ℓ = ℓm/α -1, then the FDR using BH α on (p i ) 1≤i≤m is equal to m0α m : F DR m 1 (ε BHα , p) = m 0 α m
Overlapping calibration set In the context of online anomaly detection, moving windows are classically used to capture and process the incoming data. This is why the calibration sets of the p-value family will overlap. To have a perfect control of the FDR, an upper and lower bounds is needed.

According Proposition 2, p-values with overlapping calibration sets are PRDS.

Proposition 2 (PRDS property for overlapping calibration sets). Let X i for i in 1, m be random independents variables. There are exactly m 0 random variables following the P 0 distribution, with belong to H 0 . Let Z be the random vector that combine all calibration set, all elements of Z are generated from P 0 . The set of n indices defining the elements of calibration set related to pi in Z is noted

D i . The calibration set related to X 1 is noted Z D1 = (Z i1 , ..., Z in ). For all i in 1, m : pi = p-value(X i , Z Di ).

Under these conditions, the set of p-values is PRDS on H 0

The proof of the proposition is in delayed to Appendix A.2. Since such pvalues are are PRDS, it gives an upper bound control of the FDR using Corollary 4.

Corollary 4 (PRDS property for overlapping calibration sets). Under the same conditions as Proposition 2 and the condition on calibration set cardinality satisfy

∃ℓ ≥ 1, n = ℓ m α -1: F DR(p, εBH ) ≤ m 0 α m 3.3.

Calibration set and impact of the overlap

In the context of online anomaly detection, moving windows are usually used to capture and process the incoming observations. In this context, the calibration set coincides with the observations within this window. Since successive windows are overlapping each other depending on the size of the shift, the resulting calibration sets used for computing the successive empirical p-values are also overlapping. To have a perfect control of the FDR, an upper and lower bounds are needed. Since the Appendix A.2 proves that such p-values are PRDS, it gives an upper bound control of the FDR. No theoretical results exists to compute the lower bound, indeed the existing proof in [START_REF] Marandon | Machine learning meets false discovery rate[END_REF] did not extend to overlapping calibration sets. Therefore the next discussion suggests to establish this lower bound empirically.

The following experiments aims at drawing a comparison between the FDR values in three scenarios: independent calibration sets, partially overlapping calibration sets with an overlap size driven by the value of sn (size of the shift), and the same calibration set for all empirical p-values. To be more specific, the calibration sets (and corresponding empirical p-values) were generated according to the following scheme. Each calibration set is of cardinality n. When moving from one calibration set to the next one, the shift size is equal to sn, where s in [0, 1] is the proportion of independent data between calibration sets, resulting in an overlap of cardinality (1 -s)n. Therefore an overlap occurs as long as s < 1. All these ways to build the calibration sets are called "calibration sets strategies".

1. The independent p-values (iid Cal.) are generated according to According to these three scenarios, as s increases, the overlap cardinality decreases, which results in more and more (almost) independent calibration sets. This is illustrated by the empirical results collected in Table 1. For each calibration set strategy, presented in row, and for each calibration set cardinality in column, the estimated FDR is shown. In this experiment, the reference distribution P 0 is the Gaussian N (0, 1) and the anomalies are equal to ∆ = 4. The number of tested p-values, noted m, is equal to 100 and m 1 , the number of anomalies, is equal to 1. On each sample, BH-procedure is applied with α = 0.1 and the F DP is computed. Each FDR is estimated over 10 3 repetitions. The values of n in the columns of Table 1 are chosen such that, for each pair of columns, the FDR value is smaller for the left column and larger for the right column (see Figure 2 for a visual illustration of this phenomenon). Table 1 illustrates that, in the context of the present numerical experiments, the FDR estimation is not too strongly impacted by the value of s (proportion of the overlap). To assert that the observed differences between FDR estimations in each column are not significant, permutation tests [START_REF] Fisher | The Design of Experiments[END_REF][START_REF] Phipson | Permutation P-values should never be zero: calculating exact P-values when permutations are randomly drawn[END_REF] are performed. Under H 0n hypothesis, the F DR are the same across all calibration set strategies for the calibration set cardinality n. Under H 1n there are at least two calibration strategies leading to different F DR. The F DP samples that have been used to estimate the F DR are reused. The maximal gap between sample means is used as statistic. The test is performed using the function "permutation_test" from the Python library called Scipy. The significance level is fixed at 0.05. Since multiple tests are performed over the different cardinalities, the threshold for rejecting a hypothesis is 0.00625, according Bonferroni correction. The results are display in Table 2. All tested hypotheses have a p-values greater than the threshold 0.00625. There are no significant difference in the resulting FDR between the different proportions of overlapping in calibration sets. This would suggest that considering overlapping calibration sets should not worsen too much the control of false positives and negatives.

∀i ∈ 1, m , Z i ∼ P n 0 , p1,i = p-value(X i , Z i ). ( 3

Empirical Results: Assessing the FDR control

The purpose of the present section is to compute the actual FDR value when empirical p-value are used instead of true ones. The question raised here is to check whether the FDR of the full time series is truly controlled at a prescribed level α. The empirical results must be compared with the theoretical FDR expression that has been established in Theorem 2.

In what follows, Section 3.4.1 describes the simulation design that has been considered, Section 3.4.2 details the criteria used for the assessment, and Section 3.4.3 discusses the experimental results.

Simulation design

Two scenarios have been considered to explore how much the thickness of the distribution tails can influence the results.

Thin tails:

The reference probability distribution is P 0 = N (0, 1) for normal observations and P 1 = δ ∆ N for anomalies, where ∆ N ∈ R is a parameter encoding the strength of the shift. Here δ ∆ N denotes the Dirac measure such that δ ∆ N (z) = 1 if z = ∆ N and 0 otherwise. A Gaussian reference distribution and anomalies generated from a Dirac distribution in the right tail. ∆ N is the size of the abnormal spike in the Gaussian distribution.

Thick tails:

P 0 = T (5) is a Student probability distribution with 5 degrees of freedom and P 1 = δ ∆ T denotes the alternative distribution of anomalies, where ∆ N ∈ R is a parameter encoding the strength of the shift.

Regarding the value of the shift strength in Scenarios 1 and 2, two values of ∆ N have been considered 3.5 and 4. The values of ∆ T have been chosen such that

P X∼N (0,1) (X > ∆ N ) = P X∼T (5) (X > ∆ t )
for each choice of ∆ N . This avoids any bias in the comparison of the detection power of the considered strategy depending on the ongoing scenario.

Different cardinalities have been considered for the calibration set following the mathematical expression

n ∈ {k • 10, k ∈ 1, 200 } ∪ {ℓ • 10 -1, ℓ ∈ 1, 200 }.
In particular all integers between 10 and 2 000 are explored with a step size equal to 10 as well as all integers between 9 and 1 999 with a step size of 10. This choice is justified by the particular expression of the FDR value provided by Theorem 2.

All the n elements of the calibration set are generated from the reference distribution that is, {Z 1 , . . . , Z n } ∼ P 0 . All the m observations corresponding to the tested hypotheses {X 1 , . . . , X m } are generated according to a mixture of m 1 = 1 anomalies from P 1 and m 0 = m -m 1 normal observations from P 0 . Here m = 100 and m 0 = 99.

Each simulation condition has been repeated B = 10 4 times. For each repetition 1 ≤ b ≤ B, the observations are indexed by b such that X b,j ∼ P 1 for each j ∈ 1, m 1 , and X b,j ∼ P 0 for j ∈ m 1 + 1, m .

Criteria for the performance assessment

In the present scenarios, anomalies are all located in the right tail of the reference probability distribution. Therefore the empirical p-value are computed according to Definition 2 with the scoring function a

(x) = x. For each repetition 1 ≤ b ≤ B, ∀1 ≤ j ≤ m, pb,j = p-value(X b,j , {Z b,1 , . . . , Z b,n }).
After computing the empirical p-values, the BH α procedure (see Definition 3) is applied in such a way that, for any 1 ≤ j ≤ m, d b,j = 1 BHα( pb,1 ,..., pb,m ) (j), where 1 I denotes the indicator function of the index set I. The FDP value of the sequence from 1 to m is computed from the knowledge of the true label of the observations as "normal" or "anomaly". For each repetition 1 ≤ b ≤ B,

(F DP m 1 ) b = m j=m1+1 d b,j m j=1 d b,j
.

The results obtained after the B repetitions are averaged within the FDR estimate of the sequence from 1 to m as

F DR m 1 = 1 B B b=1 (F DP m 1 ) b .
The FNR value of the sequence from 1 to m (Equation 2.2) is estimated by According to these plots, the behavior of both FDR and FNR does not exhibit any strong dependence with respect to the reference probability distribution. The results are very close for both Gaussian and Student distributions.

(F N P m 1 ) b = 1 m -m 0 m j=m0+1 d b,j , and F N R m 1 = 1 B B b=1 (F N P m 1 ) b .

Experimental discoveries

As illustrated by Figures 3a and3c, the FDR control at the prescribed level is achieved for particular values of the calibration set cardinality. These values coincide with the ones exhibited by Theorem 1, which are multiples of α/m = 10 3 (up to a downward shift by 1).

A striking remark is that the FNR curve sharply increases from 1 to n = 999. This reflects that although the FDR value becomes (close to) optimal as n increases from 1 to n = 999, the proportion of false negatives simultaneously increases leading to a suboptimal statistical performance (because of too many false negatives). Fortunately a larger cardinality n of the calibration set, for instance n = 1999, would greatly improve the results at the price of a larger calibration set, which also increases the computational cost.

Consistently with what is established in Theorem 1, the FDR value does not depend on the strength of the distribution shift ∆ as illustrated by Figures 3a As long as FDR is concerned (which is an expectation), the shift plays no role. Let us mention that focusing of the expectation does not say anything about the probability distribution of FDP, which can be influenced by the shift strength. By contrast, the comparison of Figures 3b and3d clearly shows the impact of the shift strength on the FNR value. As the shift strength becomes lower, anomalies are more difficult to be detected which inflates the FNR value.

The best cardinality n of the calibration set depends on the number m of tested hypotheses according to Theorem 1. For instance, Figure 3a shows the value n = 999 = 100/0.1 -1 as a good candidate since it achieves the desired FDR control while reducing both the number of false negatives and the computation cost. By contrast, Figure 3e rather exhibits the value n = 1499 = 150/0.1-1 as the smallest n allowing a perfect FDR control and a small number of false negatives.

Figure 3a shows other intermediate values calibration set cardinalities yielding the FDR control. For instance n = 1499 (between n = 999 and n = 1999) is predicted by Theorem 1. However complementary experiments (summarized by Figure 12 in Appendix B.1) illustrate that these intermediate values of n allowing the FDR control actually depend on the number of anomalies m 1 . Their existence can be explained by the distribution of the number of detections. For example, Figure 12d shows a high probability of detecting 3 anomalies. Assum-

ing there exists k * ∈ 1, m such that P(R(i) = k * ) ≈ 1, Theorem 2 justifies that F DR m 1 (ε BHα , p) = n n + 1 • α m 0 m + m 0 n + 1 m k=1 (1 -q n,k ) k P(R(i) = k) ≈ n n + 1 • α m 0 m + m 0 n + 1 (1 -q n,k * ) k * .
Then the proof detailed in Appendix A.4 yields that 1 -q n,k * = α m(n+1) can be reached for all ℓ ≥ 1, such that n = ℓ m αk * -1. This allows to conclude that F DR m 1 (ε BHα , p) ≈ m0α m .

How to choose the right cardinality of the calibration set?

Intuitively an optimal choice of the cardinality n of the calibration set should enable the FDR control while minimizing the number of false negatives and avoiding any excessive computation time. To achieve this objective, the first part of Corollary 1 explains that n must be chosen from the set N = {ℓm/α -1, ℓ ≥ 1}. Using the simulation scenarios described in Section 3.4.1, the aim is to visualize the relationship between the calibration set cardinality and FNR when n ∈ N . The results are summarized by Figure 4 where the FNR value is displayed versus n. For all the considered scenarios (Fig. 4a, 4b, 4c, 4d), the FNR value converges to the value reached with true p-values (horizontal dashed line) as n grows. From Figures 4a and4b, the convergence speed depends on the "difficulty" of the problem.

In practice, generating experiments similar to the ones illustrated by Figure 3.4.1 does not require to know the true reference distribution of the time series since is uses empirical p-values. However, the lack of labeled observations prevents us from computing the abnormality score and the actual FNR value, making the choice of the optimal value of n highly challenging. To tackle this challenge our suggestion is to choose the largest possible value of n that does not exceed the computation time limit. Doing that would output a value of n minimizing the FNR criterion while meeting the computational constraints. However following this suggestion does not prevent us from computational drawbacks as illustrated by Figure 4a where the FNR optimal value is reached for n = 3999 while choosing a larger n does not bring any gain (but still increases the computational costs).

Global FDR control over the full time series

While working with streaming time series data, the anomaly detection problem requires to control the FDR value of the full time series to make sure that the global false alarm rate (FDR) remains under control at the end of the iterative process. The final criterion that is to be controlled is then the global FDR criterion given by

F DR ∞ 1 (ε, p),
where ε = (ε t ) t≥1 denotes a sequence of data-driven thresholds, and p stands for a sequence of empirical p-values (see Section 3 for further details). By contrast with this global objective, anomaly detection nevertheless requires making decision at each time step that is, for each new observation, without knowing what the next ones look like. This justifies the need for another (local) criterion that will be used to make a decision at each iteration, leading to the sequence of data-driven thresholds ε = (ε t ) t≥1 . One additional difficulty results from the connection one needs to create between this local criterion and the (global) FDR of the full time series.

To this end, Section 4.1 starts by showing that controlling the FDR criterion for subseries of the full time series does not provide the desired global FDR control. Here "global" means "on the full time series" by contrast with the local FDR control, corresponding to controlling FDR for a strict subseries of the full one. Then Section 4.2 explores the connection between FDR for the full time series and the so-called modified-FDR (mFDR) for subseries. In particular, it turns out that controlling the mFDR value for all subseries of a given length m yields the desired FDR control for the full time series. Section 4.3 then explains how the classical BH-procedure can be modified to get the mFDR control for subseries of length m, while Section 4.4 illustrates the practical behavior of the considered strategies on simulation experiments.

Local and global FDR controls are not equivalent

Let us consider a time series partitioned into 4 subseries as illustrated in Figure 5. The normal points are displayed in black and the anomalies in white. The surrounded points are those that have been detected as anomalies by the procedure. When computing the number of rejections, false positives and the False positive rate for each subseries, it comes

• Subseries 1 : 2 rejections, 1 false positive. F DP 4 1 = 0.5 • Subseries 2 : 2 rejections, 0 false positive. F DP 9 5 = 0

• Subseries 3 : 1 rejections, 1 false positive. F DP 14 10 = 1 • Subseries 4 : 1 rejections, 0 false positive. F DP 19 15 = 0 If one assumes that the same probability distribution has generated the observations within each subseries, the estimated (local) FDR can be defined as average of the successive FDP values for each subseries that is, F DR 4 1 = 0.375. Let us notice that the notation emphasizes that this FDR value is the average over subseries of respective length m = 4. If one reproduces the same reasoning for the full time series, it comes: 6 rejections, two false positives, so that F DP 16 1 = 1/3 = 0.333. This example highlights that the FDP of the full time series is not equal to that of smaller subseries. This phenomenon gives some intuition on possible reasons why applying the classical BH-procedure on local windows of length m (subseries) does control the FDR criterion for the individual subseries, but does not yield the desired global FDR control for the full time series. This intuition is confirmed by the boxplots of Figure 6, where BH α has been applied on subseries of length m = 100. The left boxplot shows that BH α Figure 6: Comparison of the calculation of the FDR computed locally on a subseries and the FDR computed globally on the whole time series with Benjamini-Hochberg procedure applied on a subseries. This result is obtained by cutting a series of cardinality 1000 into 10 subseries of cardinality 100. Then the Benjamini-Hochberg procedure is applied on each subseries.

provides the desired control at level α = 10% for each individual subseries of length m. However the right boxplot clearly departs from α, meaning that the actual FDR value for the full times series of length 1 000 is strongly larger than α (more than 20% on average) leading to more false positives at the level of the full time series. The boxplots represent the quantile of FDP over 100 repetitions.

How mFDR can help in controlling the FDR of the full time series

Mixture model and time series

In this section one assumes that the time series is generated from a mixture process between a reference distribution P 0 and an alternative distribution P 1 .

The anomaly positions are supposed to be independent and generated by a Bernoulli distribution. This is a common assumption usually in the literature [START_REF] Huber | Robust estimation of a location parameter[END_REF][START_REF] Staerman | Functional anomaly detection and robust estimation[END_REF] for simplification purposes.

Definition 4 (Time series process with anomalies). Let π ∈ [0, 1] be the anomaly proportion and P 0 and P 1 be two probability distributions on the observation domain X . P 0 is the reference distribution and P 1 denotes the alternative distribution. The generation process of a time series containing anomalies (A t ) t≥0 is given, for every t ≥ 0, by

• A t ∼ B(π) (Bernoulli distribution) • if A t = 0, then X t ∼ P 0 . • if A t = 1, then X t ∼ P 1 .
Moreover given the above scheme, (X t ) t≥0 is a random process with independent and identically distributed random variables

X t ∼ (1 -π)P 0 + πP 1 .
This definition details the way anomalies are generated. In particular it assumes that anomalies are independent from each other. Let us mention that this does not prevent us from observing a sequence of successive anomalies along the time series. However this scheme substantially differs from the case analyzed by [START_REF] Gretton | A kernel two-sample test[END_REF] where specific patterns with successive anomalies are looked for.

Preliminary discussion: Disjoint and Overlapping subseries

In the context of online anomaly detection, the main focus in what follows is put on two situations where the data-driven thresholds (ε t ) t≥0 can be defined from a set of m empirical p-values: (i) the disjoint case where disjoint subseries of length m are successively considered, and (ii) the overlapping case where the subseries (of length m) successively considered share m-1 common observations at each step.

Let us start with a subseries of length m where each observation is summarized by its corresponding empirical p-value, and let us assume that there exists a function f m : [0, 1] m → [0, 1] that is mapping a set of m empirical p-values onto a real-valued random variable. This random variable corresponds to the data-driven threshold that is applied to the subseries of length m to detect potential anomalies. This function f m is called the local threshold function since it outputs a threshold which applies to a subseries of length m.

Given the above notations, the threshold sequences εd = (ε d,t ) t and εo = (ε o,t ) t can be defined as follows.

• Disjoint subseries: εd : t → εd,t is given by ∀k ≥ 0, ∀t ∈ km + 1, (k + 1)m , εd,t = f m (p km+1 , . . . , p(k+1)m ) (4.1)

• Overlapping subseries: εo : t → εo,t is given by ∀t ≥ m, εo,t = f m (p t-m+1 , . . . , pt ). (4.2)

Figure 7 illustrates these two situations. In Figure 7a, the full time series is split into small disjoint subseries of length m. f m is applied to each such subseries and the threshold is the same for all observations within a given subseries. Figure 7b displays the situation where overlapping subseries are successively considered.

Because two successive subseries differ from each other by two observations, the thresholds are different at each time step unlike the disjoint case. Furthermore the sequences εd and εo do not enjoy the same dependence properties. Figure 7a illustrates that all thresholds εd,(k-1)m+1 , . . . , εd,km are computed by applying f m to the same subseries p(k-1)m+1 , . . . , pkm . Therefore only thresholds computed from different subseries are independent, while all thresholds from the same subseries are equal. In other words, εd,t1 and εd,t2 are independent if and only if t 1 and t 2 belong to different subseries that is, ⌊t 1 /m⌋ ̸ = ⌊t 2 /m⌋, where ⌊•⌋ denotes the integer part. By contrast Figure 7b shows that the variables εo,t and εo,t-1 are dependent because they share m -1 common observations. But all of them are still different and, for each t, εo,t is independent from εo,t-m-1 . This can be reformulated as εo,t1 , εo,t2 are independent if and only if |t 1 -t 2 | > m.

In the present online anomaly detection context, considering the overlapping case sounds more convenient since the detection threshold can be updated at each time step (as soon as a new observation has been given), which makes the anomaly detector more versatile. However for technical reasons, next Section 4.2.3 still focuses disjoint subseries as a means to introduce important notions without introducing too many technicalities, while Section 4.2.4 extends the previous results to the more realistic case of overlapping subseries.

FDR control with disjoint subseries

As illustrated in Section 4.1, controlling FDR on each subseries of length m (locally) is not equivalent to controlling FDR (globally) on the full time series. However in online anomaly detection, a decision has to be made at each time step regarding the potential anomalous status of each new observation. (This is a typical instance of a local decision since at step t, the decision making process ignores what will be observed at the next step.) This requires a criterion to be optimized locally (on subseries) in such a way that the resulting global FDR value (the one of the full time series) can be proved to be controlled at the desired level α.

This requirement for a local criterion justifies the introduction of the modified FDR criterion, denoted by mFDR [START_REF] Xu | Dynamic algorithms for online multiple testing[END_REF][START_REF] Foster | α-investing: a procedure for sequential control of expected false discoveries[END_REF], which is defined as follows.

Definition 5 (mFDR). With the previous notations, the mFDR expression of the subseries from t -m + 1 to t is given by

mF DR t t-m+1 (ε, p) = E u∈H0,t-m+1≤u≤t 1[p u ≤ εu ] E t u=t-m+1 1[p u ≤ εu ]
, where ε = (ε u ) t-m+1≤u≤t denotes a sequence of thresholds, p is a sequence of empirical p-values evaluated at each observation of the subseries from t -m + 1 to t.

Mathematically the difference between the mFDR and the FDR is that the expectation is no longer on the ratio but independently on the numerator and the denominator. The main interest for mFDR is clarified by Theorem 3, which establishes its connection to FDR. To be more specific, the control of the latter at the α level provides a global control of the FDR at the same level under simple conditions. Theorem 3 (Global FDR control with disjoint subseries). Assume that εd : t → εd,t is given by εd,t = f m (p km+1 , . . . , p(k+1)m ), for any t ∈ km+1, (k+1)m (k ≥ 0) and any integer m ≥ 1 (see Eq. (4.1)). Let us also assume that the pvalue random process p = (p t ) t≥1 follows the scheme detailed in Definition 4.

Then, the global FDR value of the full (infinite) time series is equal to the local mFDR value of the any subseries of length m from

t = km + 1, k ∈ N * . More precisely, F DR ∞ 1 (ε d , p) = mF DR (k+1)m km+1 (ε d , p) = mF DR (k+1)m km+1 (ε d , p).
Since the full time series is assumed to be infinite, Theorem 3 is an asymptotic result. It gives rise to a strategy for controlling the (asymptotic) FDR criterion at level α by means of successive local controls of mFDR on small subseries of length m. According to the asymptotic nature of Theorem 3, there is no particular constraint on the integer m. However when dealing within time series of a finite length T , the Theorem 3 proof suggests that choosing an m "not too large" would be better since then, k = T /m would take large values making the LLN applicable (see for instance Eq. (4.3)). Actually in the online anomaly detection context, practitioners only have a limited freedom regarding the choice of m. Therefore, for a given fixed m, the control of the FDR value of the full time series given by Theorem 3 will be all the more accurate as T will be large. Fortunately this is not a limitation in the online anomaly detection context. The main limitation of Theorem 3 lies in the use of disjoint subseries, which sounds somewhat restrictive (at least from a practical perspective). This limitation will be overcome by next Theorem 4.

Proof of Theorem 3. Let k ≥ 1 denote an integer and T = mk. Then, the FDP definition and the A t variables introduced in Definition 4 justify that

F DP T 1 (ε d , p) = F P T 1 (ε d , p) R T 1 (ε d , p) = T t=1 1[p t < εd,t ](1 -A t ) T t=1 1[p t < εd,t ]
,

where R T 1 (ε d , p) and F P T 1 (ε d , p) respectively denote the number of rejections (resp. false positives) at the threshold εd for the subseries p.

Using the partitioning into k subseries of length m, its first comes that

F P T 1 (ε d , p) = k i=1 F P (i+1)m im-1 (ε d , p).
It is also noticeable that the k random variables {F P (i+1)m im-1 (ε d , p)} 1≤i≤k are independent and identically distributed since the thresholds εd,i remain unchanged within each subseries, they are identically distributed from one block to another, and the empirical p-values from different blocks are independent and identically distributed as well. Therefore the random variables (F P (i+1)m im-1 (ε d , p)) 1≤i≤k are independent and identically distributed, which implies (Law of Large Numbers theorem) that, almost surely,

lim k 1 k k i=1 F P (i+1)m im-1 (ε d , p) = E[F P m 1 (ε d , p)], (4.3)
where the expectation is taken over all sources of randomness. (Here it is implicitly assumed that T can go to +∞.) Repeating the argument forR T 1 (ε d , p), it also comes that

E[R m 1 (ε d , p)] = lim k 1 k k i=1 R (i+1)m im-1 (ε d , p), a.s..
The conclusion then results from noticing that

mF DR m 1 (ε d , p) = E[F P m 1 (ε d , p)] E[R m 1 (ε d , p)] = lim k k i=1 F P (i+1)m im-1 (ε d , p) k i=1 R (i+1)m im-1 (ε d , p) = lim T F DP T 1 (ε d , p).

FDR control with overlapping windows

Unlike previous Theorem

F DR ∞ 1 (ε o , p) = mF DR t t-m+1 (ε o , p) = mF DR m 1 (ε o , p)
Theorem 4 gives a similar result to the one of Theorem 3 but in a more realistic framework corresponding to the real time anomaly detection context. In particular the main improvement lies in that a threshold can be recomputed at each time step from a (shifted) subseries of length m. An important consequence is that the desired control for the FDR of the full (infinite) time series at level α can be achieved provided one can control the successive mFDR of all (shifted) subseries of length m at level α. This point is not obvious at all and constitutes the main concern of Section 4.3 where a new multiple testing procedure is designed to yield the desired control of the mFDR criterion. The main limitation of Theorem 4 is the requirement that f m has to be permutation invariant. Let us emphasize that this property holds true with the BH-procedure for instance. Let us also mention that the empirical p-values for instance computed as pt = p-value(X t , {X t-n , . . . , X t-1 }) actually satisfy the requirements of Theorem 4 regarding the independence and the stationarity.

Proof of Theorem 4. Let us start with the FDP expression for a time series of length T .

F DP T t=1 = T t=1 1[p t < εo,t ](1 -A t ) T t=1 1[p t < εo,t ] = 1 T T t=1 1[p t < f m (p t t-m+1 )](1 -A t ) 1 T T t=1 1[p t < f m (p t t-m+1 )]
,

where f m (p j i ) = f m (p i , . . . , pj ), for i < j. Since the decision process (1[p t < εo,t ]) t and the false positives process (1[p t < εo,t ]A t ) t are not independent.

Therefore it is not possible to use the Law of Large Numbers as in the proof of Theorem 3. The alternative strategy consists first in splitting the numerator and denominator into several disjoint subseries corresponding to independent and identically distributed processes. Then partitioning the times series of length

T = T ′ (n + m) into T ′ subseries, each of length n + m, it results that 1 T T t=1 1[p t < f m (p t t-m+1 )](1 -A t ) = 1 n + m n+m k=1   1 T ′ T ′ -1 t=0 1[p t(n+m)+k < f m (p t(n+m)+k t(n+m)+k-m+1 )](1 -A t(n+m)+k )   . (4.4)
Interestingly for each k from 1 to m+n, the summands within the brackets do all belong to different subseries, which makes the sum over t a sum of independent and identically distributed random variables. It results that, for each 1 ≤ k ≤ n + m, the average within the brackets is converging to its expectation by the LLN theorem.

Since the limit of a (finite) sum is equal to the sum of the limits, the average in Eq. (4.4) is converging and lim

T →∞ 1 T T t=1 1[p t < εo,t ](1 -A t ) = m k=1 E [1[p k < εo,k ](1 -A k )] a.s. = mE [1[p m < εo,m ](1 -A m )] a.s..
Then after applying the same reasoning on the denominator, it gives:

F DP ∞ 1 (ε o , p) = E [1[p m < εo,m ](1 -A m )] E [1[p m < εo,m ]]
What remains to show is to proof that:

E [1[p m < εo,m ](1 -A m )] E [1[p m < εo,m ]] = mF DR m 1 (ε o , p)
This result comes from the stationarity of p and the permutation invariance of f m . These properties imply that all p-values inside a subseries have the same probability to be rejected.

∀t ∈ 1, m , E[1[p t < f m (p m 1 )](1 -A t )] = E[1[p m < f m (p m 1 )](1 -A m )] (4.5)
Which imply

E[ m t=1 1[p t < f m (p m 1 )](1 -A t )] = mE[1[p m < f m (p m 1 )](1 -A m )] (4.6)
Using the same argument for the denominator, is gives:

E [1[p m < εo,m ](1 -A m )] E [1[p m < εo,m ]] = mE [1[p m < εo,m ](1 -A m )] mE [1[p m < εo,m ]] (4.7) = E[ m t=1 1[p t < f m (p m 1 )](1 -A t )] E[ m t=1 1[p t < f m (p m 1 )]] (4.8) = mF DR m 1 (ε o , p) (4.9)

Modified BH-procedure and mFDR control

As shown in Section 4.2, controlling the FDR value of the full time series is possible. The strategy then consists in first controlling the mFDR criterion of all successive subseries of length m along the full time series at level α. The main challenge addressed in the present section is to design a new multiple testing procedure that controls the local mFDR criterion at a prescribed level α.

In Section 4.3.1, it is proved that applying the classical BH-procedure on a time series of length m does not yield the control of mFDR at level α. However the proof of this result gives rise to a strategy for modifying the classical BHprocedure (Section 4.3.3) in a such a way that applying the so-called modified BH-procedure provides the desired mFDR control at level α.

mFDR control with the BH-procedure

Next Proposition 3 establishes the actual mFDR level achieved by the BHprocedure.

Proposition 3. Let (X i ) m

1 satisfy the requirements detailed by Definition 4. Let (p 1 , . . . , p m ) be the true p-values corresponding to a subseries of length m that is, for any

1 ≤ i ≤ m, p i = P X∼P0 (X ≥ X i ) ∈ [0, 1]. Then for every α ∈ [0, 1], applying BH α on the p-values (p i ) 1≤i≤m leads to mF DR m 1 (p) = αm 0 m ER m 1,α (1)
ER m If the ratio ER(1) m (1)/ER(1) m were known, it could be possible to control of the mF DR criterion at level α by simply applying the BH-procedure with a preliminary level α ′ = m m0 ER1,α ERα α. Unfortunately at this stage, this ratio is not known and the latter strategy cannot be straightforwardly applied. Deriving such a modified BH-procedure is the purpose of the next sections. Let us also recall that in the anomaly detection context, m 0 is unknown but expected to be close to m since only a few anomalies are usually expected. Therefore the main challenge remains to compute ER(1) m (1)/ER(1) m .

Proof of Proposition 3. The mFDR formula is given by

mF DR m 1 (p) = E[F P m 1,α (p)] E[R m 1,α (p)]
.

Let us compute the numerator E[F P m 1 (p)] value after applying the BH α . For simplification purposes, let R = R m 1,α (p) (respectively F P m 1 (p) = F P m 1,α (p)) denote the number of rejections (resp. of false positives) resulting from BH α .

Introducing D i = 1[p i ≤ R m α], it appears that F P m 1 (p) = m i=1 (1 -A i )D i .
Now focusing on D i and using R(i), the (random) number of rejections output by BH α when the p i is replaced by 0, it comes that

D i = m k=1 1[p i ≤ k m α]1[R = k] = m k=1 1[p i ≤ k m α]1[R(i) = k],
since, on the event {R = k}, p i is already rejected that is, p i ≤ k m α. Taking the conditional expectation given p except p i on both sides, it results

E[D i |p\p

i ] = m k=1 P[p i ≤ k m α]1[R(i) = k] = α m m k=1 k1[R(i) = k]
since the true p-value p i follows a uniform distribution on [0, 1]. Now integrating over all remaining p-values yields that

E[D i ] = α m m k=1 kP[R(i) = k] = α m E[R(i)].
Since A i and D i are independent random variables, it results that

E[F P m 1 (p)] = m i=1 E[(1 -A i )E[D i | A i ]] = m i=1 E[(1 -A i )E[D i ]] = α m E[R(1)] m i=1 E[1 -A i ] = α m 0 m E[R(1)],
where the last-but-one equality stems from the fact that all random variables R(i) are identically distributed.

Evaluating the ratio of rejection numbers

Previous Section 4.3.1 raises the importance of the ratio ER m 1,α (1)/ER m 1,α of rejection numbers. The present section aims at deriving a numeric approximation to this ratio. In a first step, a first result details the value of the denominator. In a second step, an approximation to the numerator is derived based on a heuristic argument and also empirically justified on simulation experiments.

Calculating the expected number of rejections When mF DR is assumed to equal α, the expected number of rejections can be made explicit.

Proposition 4.

With the previous notation, let (X 1 , . . . , X m ) be given by Definition 4, where π denotes the unknown proportion of anomalies, and assume that mF DR m 1 = α and

F N R m 1 = β ∈ [0, 1]. Then E[R m 1,α ] = mπ(1 -β) 1 -α . (4.10)
The proof is postponed to Appendix A.5. For instance, Eq. (4.10) establishes that the expected number of rejection output by BH α increases with π, the unknown proportion of anomalies along the signal. This makes sens since the more anomalies, the more expected rejections. The expected number of rejection is also increasing with α: the larger α, the less restrictive the threshold, and the more rejections should be made. However the number of rejection decreases with the FNR value β. As β increases, the proportion of false negatives grows meaning that fewer alarms are raised, which results in a smaller number of rejections.

Heuristic arguments

In what follows, the assumption is made that anomalies are easy to detect, meaning that the F N R m 1 value β is negligible compared to 1. In this context, Proposition 4 would yield that

E[R m 1,α ] ≈ mπ 1 -α . (Power)
An another assumption is also made about the relationship between

E[R m 1,α ] and E[R m 1,α (i)].
This assumption is based on a heuristic argument supported by the results of numerical experiments as reported in Table 3. In what follows, it is assumed that

E[R m 1,α (i)] = E[R m 1,α ] + 1. (Heuristic)
No mathematical proof of this statement is given in the present paper. However, Table 3 displays numerical values which empirically support this approximation, whereas further analyzing the connection between these quantities should be necessary. Let us emphasize that Table 3 has been obtained with Gaussian data (generated similarly to those detailed in Section 3.4). For all the three considered values of α, one observes that E[R(1)] remains close to (but also slightly larger than) E[R] + 1.

BHα 0.05 0.1 0.2 E[R m 1,α ] 2.
From now on and in all what follows, (Heuristic) is assumed to hold true. In particular, (Heuristic) gives rise to a strategy for computing the ratio

E[R(1)] E[R]
. So all ingredient to build a procedure that control the mFDR are given.

Modified BH

mBH α (p 1 , . . . , p m ) = BH α ′ (p 1 , . . . , p m ).
The related mBH α threshold at level α is defined as

ε mBHα = ε BH α ′ ,
when computed with true p-values, and εmBHα = εBH α ′ when used with empirical p-values.

The above definition defines the mBH α in terms of the BH-procedure by simply changing the level of control α ′ . This new level value depends on the unknown proposition π of anomalies. Since in realistic anomaly detection scenarios observations are not labeled, [START_REF] Storey | Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach[END_REF] provides guidelines on how π could be estimated. From Proposition 4, it also appears that α ′ = α(1 + 1-α mπ(1-β) ) -1 should arise as the ideal threshold. However since β is unknown, (Power) leads to the approximation suggested within the above definition.

Corollary 5 (Control of the FDR using mBH). Under the same notations and assumptions as Theorem 4. Let m and ℓ be two integers and assuming (Power) and (Heuristic) hold, let n and α ′ be defined by

α ′ = α 1 + 1-α mπ and n = ℓm/α ′ -1.
the next two results hold true.

1. If the FDR can be controlled at level α ′ on subseries of size m, using BH α ′ :

F DR m 1 (ε BH α ′ , p) = (1 -π)α ′
, then the FDR of the whole time series can be controlled at the level α using mBH α

F DR ∞ 1 (ε mBHα , p) = (1 -π)α (4.11) 2. If F DR m 1 (ε BH α ′ , p) ≤ (1 -π)α ′ .
Then, the FDR of the whole time series can be controlled at the level α

F DR ∞ 1 (ε mBHα , p) ≤ (1 -π)α (4.12)
Proof of Corallary 5. All conditions being satisfied Theorem 4 gives that:

F DR ∞ 1 (ε mBHα , p) = mF DR m 1 (ε mBHα , p) (4.13)
By definition of mBH α in Definition 6 and with Proposition 3:

mF DR m 1 (ε mBHα , p) = mF DR m 1 (ε BH α ′ , p) = E[R(1)] E[R] F DR m 1 (ε BH α ′ , p)
Under the assumptions Heuristic and Power, it gives:

E[R(1)] E[R] = 1 + 1 -α mπ mF DR m 1 (ε mBHα , p) = (1 + 1 -α mπ )F DR m 1 (ε BH α ′ , p) From hypothesis F DR m 1 (ε BH α ′ , p) = (1 -π)α ′ .
Replacing the value of α ′ with its expression it gives:

mF DR m 1 (ε mBHα ) = (1 + 1 -α mπ )(1 -π)(1 + 1 -α mπ ) -1 α = (1 -π)α
Plugin this result into Eq. 4.13 this gives desired result.

F DR ∞ 1 (ε mBHα , p) = (1 -π)α
The next result shows how mBH α can be applied to reach a global FDR control for the full time series at the desired level. It applies Corollary 5 by specifying different ways to estimate p-values on the time series.

Theorem 5 (Global FDR control using mBH α ). Let (X t ) be a mixture process introduced in Definition 4, with π denoting the anomaly proportion. Let α ∈ [0, 1] be the desired FDR level for the full time series. For m and ℓ two integers and assuming (Power) and (Heuristic) hold, let n and α ′ be defined by

α ′ = α 1 + 1-α mπ and n = ℓm/α ′ -1.
If (p t ) are given following one of the schemes

1. pt = 1 -P X∼P0 (a(X) > a(X t )) 2. pt = p-value(X t , {Z t,1 , . . . , Z t,n }) with Z t,i ∼ P 0 then, F DR ∞ 1 (ε mBHα , p) = (1 -π)α.
Otherwise, if (p t ) are given by

4. pt = p-value(X t , {(1 -A t-n+1 )X t-n+1 + A t-n+1 Z t,1 , . . . , (1 -A t )X t + A t Z t,n }) with Z t,i ∼ P 0 , then, F DR ∞ 1 (ε mBHα , p) ≤ (1 -π)α.
The main merit of Theorem 5 is to establish the actual level of control for the global FDR of the full time series depending on the type of empirical pvalue used in the anomaly detection process. This level of control depends on the unknown proportion π of anomalies. Therefore if π is close to 0 (only a few anomalies are expected), then the FDR control is close to α. However, since alpha' is close to 0 when pi is close to 0, an estimation of π (or expert knowledge) is required to detect anomalies. Some guidelines are provided in [START_REF] Storey | Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach[END_REF]. The last type of empirical p-values is (almost) the one which is used in practice in the present work. More precisely Section 5.2.3 describes empirical p-values based on a "Sliding Calibration Set".

Proof of Theorem 5. The Corollary 5 gives the two properties that the p-values families has to verify to control the FDR of the time series:

• The p-values are stationary and independent when time distance is larger than n.

• The FDR is controlled at level α ′ on subseries of size m: F DR(εBH

α ′ 1 , p) = (1 -π)α ′
In the following, these properties are verified for the different p-values.

1. The true p-value family is i.i.d., because the time series mixture is i.i.d.

Also the Theorem 1 gives the control for subseries. Using Corollary 5 the FDR of the whole time series is controlled at level (1 -π)α.

This

Empirical results

In this section, the abilities to get local control of the mFDR and the global control of the FDR, using mBH are assessed empirically. Corollary 5 and Theorem 5 give theoretical results about the control of the mFDR for subseries under the assumptions Power and Heuristic. However, these assumptions are hard to ensure in practice. In Section 4.4.1, the assessment is done on simulated data where the level of atypicity of the anomalies varies from one sample to another. Different scenarios are tested to verify if the mFDR control hold. Theorem 3 and Theorem 4 give FDR control over the full time series. In Section 4.4.2, the abilities of thresholds computed on disjoint and overlapping subseries to control the mFDR using are compared. Theorem 3 and Theorem 4 give asymptotic FDR control over the full time series. But there is no result about the speed of convergence, which is necessary when used on finite time series. In Section 4.4.3, the FDR for the full time series is calculated across different situations, as a function of time series size. It is possible to figure out when the entire series reaches control of the FDR.

Control of the mFDR on disjoint subseries

Experiment description From Corollary 5 the mBH α controls the mF DR m 1 only if the Power assumption is satisfied. Since the power of the anomaly detector depends on how easy it is to detect anomalies, the level of atypicity δ is introduced. To quantifies the atypicity of a data point X t , the true p-value is computed as p t = P X∼P0 (X > X t ), and the atypicity level is defined as the inverse of the p-value: δ t = 1/p t . The atypicity level is preferred over the p-values because it is easier to show on the x-axis of the chart, when the p-value is small. To evaluate the effect of power, for each sample all anomalies have their level of atypicity lower bounded a given parameter δ. Therefore, it is possible to observe the effect of a variation in the level of atypicity on the mF DR m 1 , F DR m 1 and F N R m 1 . For a given scenario-meaning a proportion of anomalies π, a level of atypicity δ, and a desired level of mFDR noted α-the actual mFDR, FDR, and FNR are estimated. These quantities are estimated using K = 50 samples of m data points. To control the estimation error made when estimating on a finite number of samples, each estimation is repeated B = 100 times. The estimation proceeds as follows:

1. With 1 ≤ b ≤ B, and 1 ≤ k ≤ K, m data point are generated.

• m 0 normal data p b,k,1 , . . . , p b,k,m0 are generated according the reference law U ([0, 1]).

• m 1 abnormal data p b,k,m0+1 , . . . , p b,k,m are generated using the alternative law U ([0, 1/δ]), with δ the level of atypicity of the anomalies.

2. Then, for each sample, the thresholds are estimated with BH and mBH procedures:

• εb,k,BH = BH α (p b,k,1 , . . . , p b,k,m ),

• εb,k,mBH = mBH α (p b,k,1 , . . . , p b,k,m ).

3. The number of rejections, false positives and false negatives are computed on each sample and according each threshold. Using M ∈ {mBH, BH}:

• R b,k,M = m i=1 1[p b,k,i ≤ εb,k,M ], • F P b,k,M = m0 i=1 1[p b,k,i ≤ εb,k,M ], • F N b,k,M = m i=m0+1 1[p b,k,i > εb,k,M ].
4. The FDR, mFDR and FNR are estimated by averaging results over the K samples:

• F DR b,M = 1 K K k=1 F P b,k,m R b,k,m , • mF DR b,M = K k=1 F P b,k,m K k=1 R b,k,m , • F N R b,M = 1 K K k=1 F N b,k,m m1
.

These steps are then repeated over the different scenarios.

Results and Analysis

The results are shown in Figure 8 by varying δ, α and m 1 . In Figure 8, the level of atypicity δ in represented in the abscissa. The ordinate represents the estimated mFDR (in Figure 8a or 8c) or FNR (in Figure 8b or 8d). Different colors are used to distinguish between BH and mBH procedures.

For a low level of atypicity δ, the FNR and the mFDR are high because the anomalies are difficult to detect. By increasing δ, the FNR and the mFDR decrease. As shown in Figure 8b, with values of δ around 100, the FNR is equal to 0 which can also generate a constant mFDR as shown in Figure 8a. For the mBH-procedure, the mFDR is constant and equal to α. This is consistent with Theorem 4, which guarantees the control at level α when all anomalies are detected.

Figure 8d shows the totality of the anomalies detected for δ = 2000. The same result in figure 8b with δ = 100. This is explained by the different parameters of the experiment. The easier the anomalies are detected, faster the F N R = 0 is reached for a small δ and therefore the easier it is to guarantee mF DR = α.

Conclusion

In order to control the mFDR at the desired level α using mBH, the FNR has to be equal to 0. The capacity of mBH to control the mFDR depends of the difficulty of the problem. When abnormality proportion and level of atypicity are lower, the power of mBH decreases and the mFDR is harder to control. The results of this experiment gives an idea of the atypicity level that the detector can find.

Disjoint subseries vs overlapping subseries

Experiment Description Theorems 3 and 4 theoretically prove the control of the FDR over the full time series throw control of the mFDR over disjoint subseries or overlapping subseries. According Corollary 5, the procedure mBH α allows the control of the mFDR over subseries under assumption Heuristic and Power that are hard to verify. Empirical results from Section 4.4.1 show that control of mFDR for the disjoint subseries can be obtained for scenarios where the level of atypicity δ is high enough. It still unknown whether these results hold true in cases where the subseries overlap In this section FDR control throw disjoint and overlapping subseries are compared.

For each scenario, the quantities mF DR m 1 and F N R m 1 are estimated two times, using disjoint subseries and using overlapping subseries. All subseries are extracted from the same time series of size T = 10 4 . The distribution of these estimations is obtained by repeating the experiment across B = 100 time series. Thus, the two estimations of mF DR m 1 and F N R m 1 quantities can be compared. The experimental design is described as follows:

1. With b in 1, B and t in 1, T , the time series is generated from a mixture model:

• A b,t ∼ Ber(π)

• If A b,t = 0, p b,t ∼ U ([0, 1]) • Otherwise: p b,t ∼ U ([0, 1/δ])
2. The thresholds of mBH are estimated on each subseries p b,k,t+1 , . . . , p b,k,t+m0+m1 :

• εb,t = mBH α (p b,k,t+1 , . . . , p b,k,t+m0+m1 ).

3. The numbers of rejections, false positives and false negatives are calculated, according the different cases.

(a) In the disjoint subseries case, the quantities are computed using only thresholds on the form εb,km over disjoint subseries For 1 ≤ b ≤ L and 1

≤ k ≤ K = T /m: • R b,k,d = (k+1)m t=km+1 1[p b,k,t ≤ εb,km ], • F P b,k,d = (k+1)m t=km+1 1[p b,k,t ≤ εb,km ](1 -A t ), • F N b,k,d = (k+1)m t=km+1 1[p b,k,t > εb,km ]A t .
The mFDR and FNR are estimated:

• mF DR b,d = 1 K K k=1 F P b,k,m,d 1 K K k=1 R b,k,m,d , • F N R b,d = 1 K K k=1 F N b,k,m,d m1
.

(b) In the overlapping subseries case, the quantities are computed using the thresholds from all overlapping subseries εb,t :

For 1 ≤ b ≤ L and 1 ≤ k ≤ K = T /m: • R b,k,o = (k+1)m t=km+1 1[p b,t-m+1,t,o ≤ εb,t ], • F P b,k,o = (k+1)m t=km+1 1[p b,t-m+1,t ≤ εb,t ](1 -A t ), • F N b,k,o = (k+1)m t=km+1 1[p b,t-m+1,t > εb,t ]A t .
Notice the difference with disjoint windows case, all p-values of a subseries are compared to different thresholds and not to the same εb,km .

The mFDR and FNR are estimated:

• mF DR b,o = 1 K K k=1 F P b,k,m,o 1 K K k=1 R b,k,m,o , • F N R b,o = 1 K K k=1 F N b,k,m,o m1
.

Different scenarios are generated by varying the proportion of anomalies π and the atypicity level δ.

Results and analysis

As shown in Figure 9, disjoint and overlapping subseries control give similar results in mFDR and FNR for considered cases. Indeed, the curves are indistinguishable and decrease at the same rate.

Conclusion

The FDR control quality are similar for both strategies, overlapping windows and disjoint windows. This imply that performances of the anomaly detector to not decrease by using overlapping windows instead of disjoint windows. This is a practical result that allows to do real time detection without having to wait to complete disjoint windows.

Convergence of false discovery rate control

This section studies the convergence rate of the FDR over the full time series using mBH α .

Experiment Description

The theoretical results obtained in Theorem 4 only guarantee an asymptotic control of the FDR on the whole time series.

In practice, it is more useful to have a control of the FDR at any time, i.e. on subseries of finite size. The question is empirically studied by observing the speed of convergence of the false discovery rate towards the level α. The FDR of the full time series is calculated across different scenarios, as a function of time series size. In order to get the distribution of the FDR, the experiment is repeated on B = 100 time series. The maximal time series size explored is T = 10 4 .

1. For 1 ≤ b ≤ B and for 1 ≤ t ≤ T :

• A b,t ∼ Ber(π)

• If A b,t = 0, p b,t ∼ U ([0, 1]) • Otherwise: p b,t ∼ U ([0, 1/δ])
2. The thresholds are estimated with mBH α : εb,t,α = εmBHα (p b,t-m+1 , . . . , p b,t )

3. The proportion of false discovery (FDP) on the partial time series are calculated:

F DP b,t,α = t u=1 (1 -A b,t )1[p b,t ≤ εb,t,α ] t u=1 1[p b,t ≤ εb,t,α ]
Different scenarios are generated by varying the proportion of anomalies π and the atypicity level δ.

Results and analysis

In Figure 10, the false discovery proportion is represented in the ordinate according to the size of the time series given in the abscissa. The different levels of α used to compute mBH threshold are experimented with the results of the median FDP and its 95% band is shown in different colors. Different scenarios are represented by varying the proportion of anomalies between the sub figures.

It can be observed that the convergence is quite fast from a size of 2000 data points, since for a α of 0.05, is has 95% chance to have a false positive rate between 0.04 and 0.06, on Figure 10. Thus, the control of the false positive rate, can be ensured with a high probability, for a series of one data point per minute recorded over a few days, .

Conclusion

This ensures that the control at level is reached not only for infinite time series but also for finite time series which allows our model to be used in practice.

Empirical simulation against competitor

The control of the FDR with p-values estimated empirically has been studied at Section 3. Theorem 2 ensure the control of the F DR m 1 when the p-values are estimated on calibration set having particular cardinality value. Theorems 3 and 4 ensure the control of the FDR of the full time series throw control of the mF DR m 1 of the subseries. Corollary 5 enables to deduce that the mBH α procedure can be apply to control the FDR of the full time series under the Heuristic and Power assumptions. Even though these assumptions are hard to ensure theoretically, the experiment at Section 4.4.1 shows that the mF DR m 1 is controlled for tested scenario, provided that anomalies are sufficiently atypical. Experiment from Section 4.4.3 shows that the control of the FDR is possible even the time series is not infinite as required by Theorem 4.

These different results provide the conditions for building an anomaly detector that controls the FDR of the time series through control of the mFDR on the subseries and the p-empirical value. Our anomaly detector is evaluated under different scenarios by varying the generated anomalies and the targeted FDR. To understand the source of the difficulties that the anomaly detector may encounter, different sequences of p-values with oracle information are introduced. Our anomaly detector is compared against Levels based On Recent Discovery (LORD) which is a online multiple testing procedure, introduced in [START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF] to control the FDR.

Data

The synthetic data are generated from Gaussian distribution. With the use of the empirical p-value estimator there are no need to evaluate on other data distribution. Only anomaly proportion and the distribution shift associated to anomalies impact the performances of the anomaly detector. Data are generated accordingly with Definition 4 with Gaussian reference distribution and anomaly spike like in Section 3.4.1. The strength of the distribution shift noted by ∆ takes value in {3σ, 3.5σ, 4σ} and the abnormality proportion noted π is equal to 0.01. Each generated time series contain T = 10 4 data points. Each experiment is repeated over 100 time series.

For t in 1, T :

• A t ∼ B(π) • X t = N (0, σ 2 ) if A t = 0 ∆σ else
The value of ∆ represents the atypicity score of the anomalies. Anomalies with higher ∆ are easier to detect. In this experiment, the standard deviation σ is set to 1.

Threshold and p-value estimators description

Our proposal mBH on overlapping subseries

Using the p-value with the empirical estimator, the anomalies are detected by using mBH as the threshold estimator on overlapping subseries in the Algorithm 1. For each time t, the threshold is computed as: εmBHα,t = f m (p t-m , . . . , pt ), where f m is the mBH-procedure. To ensure FDR control according to Theorem 1, the cardinality of the calibration set to be equal to n = m α -1. In this experiment m is equal to 100 and α takes values 0.1 and 0.2 depending the tested scenario. So the calibration set takes values 999 or 1999.

LORD

LORD introduced in [START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF] is based on alpha-investing rules to define a threshold on p-values. For each time t the threshold is computed from according to the alpha-investing rules, depending on previous decision made by the algorithm. For more precision refer to the original article [START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF]. The empirical p-value specified in Definition 2.4 does not respect this property while the conformal p-value, defined in Equation 2.3.1 respects this property. Using conformal p-values to apply LORD algorithm leads to a weak power detecting anomalies. The issue is that p ≥ 1 n+1 is always verified and the threshold sequence εt decreases quickly when no rejection are made. No anomaly can be detected. For these reasons, the empirical p-value introduced in Equation 2.4 is used while applying LORD and mBH. In this experiment LORD3 from [START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF] is used with the same parameters as in the original paper.

p-value estimation

Different sequences of p-values are used to understand the limitations of our anomaly detector. The true p-values are used to evaluate the case where the only limitation comes from the multiple testing procedure. One can thus understand how the estimation of the p-values affects the detection of anomalies. One way to estimate p-values in practice is to use the same calibration set for all p-values. This is referred as the fixed calibration set. However, the p-values may be biased in that particular calibration set. In practice, the usual way to implement the estimated p-values is to use a sliding calibration set. To evaluate the p-value of a data point X t , n preceding data points are used as a calibration set. To a bias in the estimation, the points detected as abnormal cannot be part of the calibration set. However, the calibration set can be biased by undetected anomalies. To evaluate this impact, the sliding calibration set-⋆ is introduced, where the knowledge oracle of the labels is used to construct the calibration set from the previous data points.

The different p-value sequences are computed as follows:

• Oracle: The true p-value is used instead of the estimated one.

∀t ∈ 1, T , pt = Φ(X t )

• Fixed calibration set (Fixed Cal.): The p-value is estimated using the same calibration set {Z i , i ∈ [1, n]} for all observations.

∀t ∈ 1, T , pt = 1 n n i=1 1[Z i > X t ]
• Sliding Calibration set-⋆ (Sliding Cal.-⋆): The p-value is estimated using a calibration that is a sliding windows containing the n previous true normal data.

∀t ∈ 1, T , pt = 1 n n i=1 1[X h(t,i) > X t ]
With h the function that select observation that respect H 0 . For each t and i, h(t, i) gives the i-th observation lower than t and that respect H 0 hypothesis.

• Sliding Calibration set (Sliding Cal.): The calibration set is a sliding windows containing the n previous estimated normal data.

∀t ∈ 1, T , pt = 1 n n i=1 1[X ĥ(t,i) > X t ]
With ĥ the function that estimates the function h. For each t and i, ĥ(t, i) give the i-th observation lower than t and d ĥ(t,i) = 0.

Performance metrics

The anomaly detector are evaluated using their ability to control the FDR and minimize the FNR of the full time series. Therefore, the two applied metrics are the FDP and the FNP computed as:

F DP = T t=1 1[p t < εt ](1 -A t ) T t=1 1[p t < εt ]
and

F N P = T t=1 1[p t < εt ](1 -A t ) T t=1
A t where εt is estimated using mBH or LORD and pt is estimated using one of the estimator defined in Section 5.2.3. 3. and the distribution shift between the normal data and anomalies, noted ∆, varying from 4σ to 3σ. 4b needs to be taken into consideration. When ∆ = 4σ, the FNR is close to 0, while when ∆ = 3σ the FNR is equal to 0.793. Similar results are obtained with other test configurations in Table 4c and Table 4d. The FDR control at the desired level need the FNR to be close to 0.

Results

The box plots shown in

Effect of p-value estimation

To understand how the p-value estimation can prevent the control of the FDR, the first four rows in Table 4a are compared. In the column "4σ", the FDR values for the configurations "Oracle PV", "Fixed Cal." and "Sliding Cal.-⋆" are very close to the desired level α = 0.1. This control is enabled by Theorem 5, since the p-values verify all hypotheses, in particular all data in the calibration sets are generated according to the reference distribution. However, in the case of "Sliding Cal.", the FDR increases at a value of 0.335. For the same configurations, the FDR remains low, between 0.019 and 0.040 as shown at Table 4b. The increase of FDR when using "Sliding Cal." instead of "Sliding Cal.-⋆" is a consequence of calibration set contamination. Indeed, according to the procedure used to build the calibration sets, described in Section 5.2.3, all detected anomalies are removed from calibration sets used in the estimation of next p-values. When an observation is wrongly detected as an anomaly, this data point cannot be part of the calibration set at future steps of the online detection. Instead, it is replaced by an other data point having statistically a lower atypicity score. Indeed false positives have high atypicity score to be (wrongly) detected as anomalies. As a result, the calibration set contains data points with lower scores than if it had been generated under P 0 . It leads to underestimate the p-values and to increase the number of false positives. This illustrates the major drawback of mBH: it is highly sensitive to the non robustness of the pvalue estimator. Figure 15a shows that using fixed calibration instead of sliding calibration-⋆ gives a larger variance on the FDP while the FDR is the same. Using a single calibration set for the entire time series means that the FDP is highly dependent on the start of the time series. By modifying the calibration set at each time step, the statistical fluctuations in the FDP are smoothed over the course of the time series analysis.

Comparison with LORD

In this section, the results found using mBH and the ones using LORD are compared. As known from the literature, LORD controls the FDR of super-uniform p-values. In this experiment, the question is in the capacity of LORD method to control the FDR of empirical p-values that have no theoretical guaranties. It can be noticed in Figure 4a that LORD is able to ensure the control of the FDR for all calibration set definitions when anomalies are easier to detect as for ∆ = 4σ or ∆ = 3.5σ. In particular, unlike mBH, LORD is able to control the FDR in the case of the sliding calibration set. However, mBH method has a lower FNR compared to the LORD method, as shown in 4a and 4b. For example, Table 4b shows that the FNR is equal to 0.019 with mBH while it is equal to 0.781 with LORD, in the case using Sliding Calibration set-⋆ on data having ∆ = 3σ. Nevertheless, with the Sliding Calibration set case, the LORD method has quite the same FNR but with lower FDR (0.335 against 0.075). The contamination issue of mBH offsets the superior performance observed in the Sliding Calibration set-⋆.

Conclusion

In this article, an online anomaly detector that aims to have a better control of the FDR at a given level α has been proposed. The research has been developed to tackle two issues:

• the empirical p-values: it ensures conditions on the calibration cardinality to ensure FDR control when using Benjamini-Hochberg.

• and the online detection: it ensures a global control of the FDR through local control of the mFDR of subseries, using a modified version of the BH-procedure.

The results of our research is the assessment of our proposal from the theoretical point of view and from empirical experiments. Our method has been compared with a method from the state of the art. It shows the strong capability for ensuring control of the FDR even in the case of empirical p-values. The major drawback and improvement path of our method is it relies on non-robust p-value estimation.

Let us now introduce the random variables R(i) that are the number of rejections generated by BH when p i is replaced by the value 0 that is, R(i) = BH α (p 1 , .., p i-1 , 0, p i+1 , . . . , p m ). It results that

1[p i ≤ αk m ]1[R = k] = 1[p i ≤ αk m ]1[R(i) = k],
since, on the event {p i ≤ αk m }, p i is rejected and therefore R = R(i). Let us also notice that the independence between the p-values is already used at this stage since modifying the value of p i does not affect that of the others.

A.4. Proof of Corollary 1

Proof of Corollary 1. To get a deeper understanding of the FDR expression obtained in Theorem 2, q n,k the fractional part of αkn m is introduced:

q n,k = αkn m - αkn m
When plugged into the FDR expression, it gives:

F DR = m 0 m k=1 αkn m +1-q n,k n+1 k P(R(1) = k) F DR = m 0 α m n n + 1 + m 0 n + 1 m k=1 1 -q n,k k P(R(1) = k) (A.7)
In order to get lower and upper bounds of the FDR, the value of q n,k should be expressed as a function of α, k, n and m.

For the next part of the proof, it is useful to express the relation between q n,k and q n+1,k . It gives the effect of increasing the cardinality of the calibration by one. Using the definition of the fractional part:

q n+1,k -q n,k = αk(n + 1) m - αk(n + 1) m - αkn m + αkn m q n+1,k -q n,k = αk m - αk(n + 1) m + αkn m
Which can be expressed as a congruence relation:

q n+1,k -q n,k ≡ αk m (mod 1) (A.8)
Two cases are studied:

1. Particular case: there exists an integer 1 ≤ ℓ such that ℓm α is an integer. the notation n ℓ = ℓm α is introduced. Since: αkn ℓ m = αkℓm/α m = kℓ is an integer, then the fractional part is null:

q n ℓ ,k = 0
If the calibration set cardinality n is equal to n = n ℓ -1 = ℓm α -1. Then, the congruence relation in Eq. A.8 gives:

q n ℓ -1,k ≡ q n l ,k -αk/m (mod 1) q n ℓ -1,k ≡ 0 -αk/m (mod 1)

Using the fact that fractional part of a number belongs to [0, 1[, the only possible value to q n ℓ -1,k is:

q n ℓ -1,k = 1 -αk/m
Plugging the value of q n ℓ -1,k into Eq. A.7, it gives: It implies that ⌊ αk(n ℓ -1) m ⌋ < ⌊ αk(n ℓ ) m ⌋. Also, Eq. A.8 is expressed as q n ℓ ,kq n ℓ -1,k ≡ αk m (mod 1):

1 -αk m ≤ q n ℓ -1,k < 1 (A.9) Indeed, the fractional part of a number as to be larger than 1 -αk/m so that adding αk/m increase the integer part.

By plugin the bounds of q n ℓ -1,k into Eq. A.7, it can gives the bounds of the FDR. At first, to compute the upper bound of the FDR the lower bound of q n ℓ -1,k is used: With the same calculations as for the "Particular case", it gives:

F DR ≤ m 0 α m
Similarly, the lower bound of the FDR can be obtained using the q n,k upper bound from Eq. A.9 plugged into Eq. A.

(2. 5 )

 5 This definition implies the p-value property for all u in [0, 1]. But this estimator is less powerful, as illustrated by Figure11in Appendix A.1 where the FNR resulting from the use of conformal p-values is always larger than that of empirical p-values.

Figure 2

 2 Figure 2

Figure 3

 3 Figure 3 displays the FDR value (left panel) and the FNR value (right panel) as a function of the calibration set cardinality for the two scenarios (Gaussian and Student) described in Section 3.4.1. The blue (respectively orange) curve corresponds to the Gaussian (resp. Student) reference distribution. The horizontal line is the prescribed level α = 0.1 at which FDR should be controlled with true p-values (Theorem 1). Figures 3a and 3b are obtained with ∆ N = 4, while Figures 3c and 3d result from ∆ N = 3.5.

  (a) FDR depending on n with 4 sigmas anomalies (b) FNR depending on n with 4 sigmas anomalies (c) FDR depending on n with 3.5 sigmas anomalies (d) FNR depending on n with 3.5 sigmas anomalies (e) FDR depending on n with 4 sigmas anomalies and m = 150 (f) FNR depending on n with 4 sigmas anomalies and m = 150

Figure 3 :

 3 Figure 3: Effect of calibration set cardinality and abnormality score on the FDR control and the FNR

(a) ∆ = 4 ,

 4 m1 = 1 (b) ∆ = 3.5, m1 = 1 (c) ∆ = 3, m1 = 1 (d) ∆ = 3, m1 = 5

Figure 4 :

 4 Figure 4: FNR as a function of the calibration set cardinality constrained to belong to N .

Figure 5 :

 5 Figure 5: Illustration of anomaly detection in subseries.

  p1 , . . . , p(k-2)m+1 , . . . , p(k-1)m , p(k-1)m+1 , . . . , pkm fm fm εd,1 , . . . , εd,(k-2)m+1 , . . . , εd,(k-1)m , εd,(k-1)m+1 , . . . , εd,km (a) Illustration of disjoint windows. p1 , . . . , pt-2m , . . . , pt-m-1 , pt-m , pt-m+1 , . . . , pt-1 , pt fm fm fm εo,1 , . . . , εo,t-2m , . . . , εo,t-m-1 , εo,t-m , εo,t-m+1 . . . , εo,t-1 , εo,t (b) Illustration of Overlapping sliding windows.

Figure 7 :

 7 Figure 7: Comparison of disjoint window and overlapping window for the threshold function.

  (1) = |BH α (0, p 2 , .., p m )| and R m 1,α = |BH α (p 1 , p 2 , .., p m )|, with |S| denoting the cardinality of the set S.

From previous Sections 4

 4 .3.1 and 4.3.2, it is now possible to suggest and analyze the new modified BH-procedure (mBH in the sequel). Definition 6 (Modified BH-procedure (mBH)). Let m be an integer and α ∈ [0, 1]. Let us introduce the level α ′ = α(1 + 1-α mπ ) -1 , where π ∈ [0, 1] denotes the unknown proportion of anomalies (see Definition 4. Let us further assume (Power) and (Heuristic) hold true. Then the modified BH-procedure, denoted by mBH α , is given for all true p-values (p 1 , . . . , p m ) ∈ [0, 1] m by,

  (a) mFDR, α = 0.2, π = 0.07 (b) FNR, α = 0.2, π = 0.07 (c) mFDR, α = 0.1, π = 0.01 (d) FNR, α = 0.1, π = 0.01

Figure 8 :

 8 Figure 8: mFDR and FNR as a function of level of atypicity across different scenarios

  (a) mFDR, α = 0.1, π = 0.01 (b) FNR, α = 0.1, π = 0.01

Figure 9 :

 9 Figure 9: Comparison of mFDR and FNR control with disjoint and overlapping windows method.

Figure 10 :

 10 Figure 10: FDR over the full time series as a function of the time series size.

   represent the FDP and FNP distribution for 1000 repetitions. Inside each sub figure (a, b, c, d, e,..), the box plot distributions are displayed according to: 1. the multiple testing method mBH or LORD, 2. the p-value estimation model set to Oracle PV, Fixed Cal., Sliding Cal.-⋆ or Sliding Cal.

  By combining the previous argument and the independence between R(i) and the other p-values, the expectation on both sides yieldsi ≤ αk m ]1[R(i) = k] k ⇒ F DR = E[F DP ] =equality results from the fact that the true p-values follow a uniform distribution on [0, 1]. The result finally follows from noticing that for each 1 ≤ i ≤ m 0 , m k=1 P[R(i) = k], since R(i) ≥ 1 by definition.

2 .

 2 General case: With α ∈]0, 1], for each ℓ the notation n ℓ = ℓm α is introduced. Notice that this definition is consistent with the particular case. The ceiling function definition gives:Multiplying by αk on each side and the n ℓ notation:αk(n ℓ -1) m < kℓ ≤ αk(n ℓ ) m

F DR ≤ m 0 (n ℓ - 1

 1

B. 2 . 01 Figure 13 :B. 3 .

 201133 Figure 13: Effect of the atypicity level on the mFDR, FDR and FNR, according to different multiple testing procedures.

Figure 14 :

 14 Figure 14: Effect of atypicity level on mFDR and FNR, depending on whether detection is on disjoint or overlapping subseries

Figure 15 :

 15 Figure 15: Comparison of the FDPs acquired from using different multiple testing procedures, mBH or LORD, and from the way the p-values are calculated, in the case α = 0.1.

Figure 16 :

 16 Figure 16: Comparison of the FNPs acquired from using different multiple testing procedures, mBH or LORD, and from the way the p-values are calculated, in the case α = 0.1.

Figure 17 :Figure 18 :

 1718 Figure 17: Comparison of the FDPs acquired from using different multiple testing procedures, mBH or LORD, and from the way the p-values are calculated, in the case α = 0.2

Table 1

 1

		FDR results with overlapping calibration sets
	n	249	250	499	500	749	750	999	1000
	Same Cal.	0.164 0.175	0.112 0.183	0.137 0.131	0.097 0.154
	Over Cal. (s=0.1%)	0.167 0.174	0.100 0.156	0.138 0.125	0.093 0.140
	Over Cal. (s=0.2%)	0.162 0.176	0.095 0.170	0.124 0.127	0.109 0.143
	Over Cal. (s=0.5%)	0.163 0.166	0.110 0.170	0.116 0.132	0.111 0.149
	Over Cal. (s=1%)	0.151 0.180	0.094 0.177	0.127 0.128	0.099 0.143
	Over Cal. (s=2%)	0.164 0.180	0.108 0.179	0.133 0.140	0.097 0.143
	Over Cal. (s=5%)	0.168 0.172	0.108 0.169	0.125 0.130	0.096 0.144
	Over Cal. (s=10%)	0.165 0.181	0.104 0.185	0.122 0.140	0.105 0.146
	Over Cal. (s=20%)	0.173 0.207	0.109 0.171	0.136 0.149	0.101 0.140
	Over Cal. (s=50%)	0.180 0.187	0.103 0.183	0.121 0.128	0.094 0.143
	iid Cal.	0.171 0.188	0.115 0.174	0.138 0.143	0.104 0.132

Table 2

 2

		p-values resulting from permutations test
	n	249	250	499	500	749	750	999	1000
	p-value of the test	0.300 0.0326	0.572 0.313	0.588 0.435	0.735 0.690

  3, following Theorem 4 establishes a similar control of the global FDR criterion on the full time series by means of successive local controls of the mFDR criterion on subseries that are allowed to overlap each other. This is closer to the practical situation arising in online anomaly detection where one new observations is collected at each time step, inducing a shift by one of the set of observations for which a decision has to be made.

Theorem 4 (Global FDR control with overlapping subseries). Assume that εo : t → ε o,t is given by εo,t = f m (p t-m+1 , . . . , pt ), for any t ≥ 1, with f m : [0, 1] m → [0, 1] permutation invariant (see Eq. (4.1)). Let us also assume that the p-value random process p = (p t ) t≥1 follows the scheme detailed in Definition 4 and there exists n such that |t 1 -t 2 | > n implies that pt1 and pt2 are independent. Then, the global FDR value of the full (infinite) time series is equal to the local mFDR value of the any subseries of length m computed at time t ∈ N * . More precisely

Table 3

 3 Numerical evaluations for different values of α (10 3 repetitions)

		14	2.32	2.78
	E[R m 1,α (i)]	3.18	3.44	3.99

  p-value family is i.i.d., because the time series mixture is i.i.d. and the calibration sets are also i.i.d. Using the Theorem 1, it gives the FDR control on subseries. With Corollary 5 the FDR of the whole time series is controlled at level (1 -π)α. 3. This p-value family is not i.i.d. However, because the calibration are build using a sliding window of size n, two p-values p t1 and p t2 are independent when |t 1 -t 2 | > m. The calibration sets of the p-values overlapped, Then Corollary 4 ensures the upper bound of the FDR for subseries of size m.

With Corollary 5 the FDR of the whole time series is upper bounded controlled at level (1 -π)α.

Table 4

 4

	gives a summary using FDR and FNR estimations. It enables eas-
	ily the comparison between these values coming from the different strategies
	combining:
	1. the multiple testing method mBH or LORD,
	2. the choice of the level α varying from 0.1 to 0.2,
	3. the p-value estimation model set to Oracle PV, Fixed Cal., Sliding Cal.-⋆
	or Sliding Cal.
	4. and the distribution shift between the normal data and anomalies, varying
	from 4σ to 3σ.

Table 4

 4 Comparison of mBH versus LORD for online anomaly detection in Gaussian white noisewith different abnormality levels. Power from Theorem 5, mBH α enables control of the FDR at level α if all anomalies are detected. To test this assertion, the different columns of the Table4a, are compared. In the row "mBH with Oracle PV", with ∆ = 4σ the FDR is estimated at 0.101 which is close to the desired level α = 0.1. While, when ∆ = 3σ the FDR level is estimated at 0.281 which is almost three times the desired level α. The FNR results in Table

	5.5. Analysis
	5.5.1. Effect of the strength of the distribution shift ∆
	According the assumption

Acknowledgments

The authors would like to thank Cristian Preda, director of the MODAL team at Inria, for valuable discussions.

Appendix A: Proofs

A.1. Comparison of p-values estimators

The control of the FDR is not achievable using classical multiple testing [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF][START_REF] Ramdas | SAF-FRON: an adaptive algorithm for online control of the false discovery rate[END_REF] since the p-value estimator, shown in Definition 2, is not super-uniform. Conformal p-value estimator p, shown in Equation 2.3.1, verifies the superuniform property. However, this estimator p ≥ 1 m+1 has lower power because zero anomalies are detected with thresholds below 1 m+1 . Figure 11 To conclude, the choice between conformal p-values and estimated p-values depends on the calibration set size. Indeed, for calibration set n = 1000 the performances are similar. But for other calibration set sizes as n = 1500 the FDR control are similar but the FNR is better for estimated p-values.

A.2. PRDS property for p-values having overlapping calibration set

The following construction is used to describe a family of p-values with overlapping calibration set. Let Z the vector that combine all calibration set, the Z i are i.i.d. with marginal probability P 0 . The set of the n indices defining the elements of the calibration set related to pi in Z is noted D i . The calibration related to

To proof that p-values with overlapping calibration set are PRDS as described in Definition 1, the methodology used in [START_REF] Bates | Testing for outliers with conformal p-values[END_REF] to be extended in the case of overlapping calibration set. For i in 1, m the calibration set associated to X i is noted Z Di . The law of total probabilities gives:

If these two lemma are suppose to be true, the PRDS property is verified. Lemma 5.1. For non-decreasing set A and vectors z, z ′ such that z ⪰ z ′ , then

Indeed, take i ∈ 1, m and u ≥ u ′ and define Z Di,1 and Z Di,2 as in the statement of Lemma 5.2.

The PRDS property is satisfied. To complete the proof, the introduced lemmas are proven.

Proof of Lemma 5.1. Let be i in 1, m and vectors z, z ′ and z vectors such that z ⪰ z ′ . The vectors z, z ′ are used to define the calibration set related to the p-values pi and z is used to define elements of calibrations sets that are not in the calibration set of pi . By conditioning on the calibration sets defined by (z, z) and (z ′ , z) it gives:

This result comes from the decomposition the following decomposition, for all j in 1, m

2 can be integrated over Z Di to give: n+1) the order statistics of (a(Z Di,1 ), . . . , a(Z Di,n ), a(X i )). And R i the rank of a(X i ) among these. [START_REF] Ahmed | A survey of network anomaly detection techniques[END_REF] , . . . , S ′ i,(n+1) :

(S (1) , . . . , S (n)

The right-hand side is not increasing with k and pi = Ri-1 n

A.3. Proof of Theorem 1

Proof of Theorem 1. Let R be a random variable describing the number of rejections made by BH α that is, R = m i=1 D i , where D i = 1 if hypothesis H 0,i is rejected. Let also F P be the number of false positives made by BH α . Then,

where A i is a random variable equal to 1 if hypothesis H 0,i is true and 0 otherwise. Furthermore

A.5. Proof of Proposition 4

Proof of Proposition 4. By definition mF DR m 1 =

ER m 1 , and R m 1 = F P m 1 + T P m 1 . With hypothesis the mF DR is equal to α, this gives:

Then, the expectation of true positives is expressed using the proportion of false negatives β, the proportion of anomaly π in the m observations, A i the random variable equal to 1 if the observation X i is an anomaly and D i the random variable equal to 1 if the observation X i is detected as anomaly :