Détection de ruptures

En analyse statistique, le problème de détection de ruptures (ou détection de points de changement) est un problème de régression ayant pour but d'estimer les instants où un signal présente des changements dans la distribution. Ces instants sont matérialisés sur les deux figures par des lignes verticales bleues. Classiquement, on réalise de la détection de rupture pour un signal ayant des changements dans la moyenne. De manière plus générale, on réalise de la détection de ruptures pour un signal ayant des changements dans la distribution (par exemple, dans la moyenne et la variance).

La détection de ruptures peut s'appliquer à un signal sonore d'une émission dont on souhaite estimer les instants où l'on change de situations, à la détection d'attaques informatiques (variations de flux réseaux)^[1], ou encore au contrôle qualité.

Cet article traite du problème de détection de ruptures rétrospective (dite offline) où l'on dispose de l'ensemble des données du signal. Ce contexte est différent d'une détection temps réel (dite online) où les données du signal arrivent à la volée, mais moins à même de détecter précisément l'instant de rupture.

Formalisation du problème

Soit ${\displaystyle \{X_{t}\}_{t\in [\![1,T]\!]}}$ un signal réel, provenant d'observations recueillies au cours des instants ${\displaystyle 1,\dots ,T}$, présentant des changements dans la distribution. En notant ${\displaystyle P_{X_{i}}}$ la loi de probabilité de ${\displaystyle X_{i}}$, la distribution de ${\displaystyle \{X_{t}\}_{t\in [\![1,T]\!]}}$ vérifie :

${\displaystyle {\begin{aligned}P_{X_{\tau _{0}^{*}}}=\dots =P_{X_{\tau _{1}^{*}-1}}\neq P_{X_{\tau _{1}^{*}}}=\dots =P_{X_{\tau _{2}^{*}-1}}\neq \dots \neq P_{X_{\tau _{D^{*}-1}^{*}}}=\dots =P_{X_{\tau _{D^{*}}^{*}-1}},\end{aligned}}}$

avec ${\displaystyle \tau _{0}^{*},\dots ,\tau _{D^{*}}^{*}}$ étant les vrais instants de ruptures (on note ${\displaystyle m^{*}=\{\tau _{0}^{*},\dots ,\tau _{D^{*}}^{*}\}}$ (${\displaystyle D^{*}}$ est le vrai nombre de segments) avec la convention ${\displaystyle \tau _{0}^{*}=1}$ et ${\displaystyle \tau _{D^{*}}^{*}=T+1}$). On cherche à estimer ces instants de ruptures à l'aide d'un algorithme.

Dans le cas de la détection de rupture dans la moyenne, le modèle est :

${\displaystyle {\begin{aligned}\forall t\in [\![1,T]\!],~X_{t}=s_{t}+\epsilon _{t},\end{aligned}}}$

avec ${\displaystyle \{s_{t}\}_{t\in [\![1,T]\!]}}$ est la fonction de régression et ${\displaystyle ~\{\epsilon _{t}\}_{t\in [\![1,T]\!]}}$ est un bruit d'espérance nulle et de variance ${\displaystyle \sigma ^{2}}$. La fonction de régression ${\displaystyle \{s_{t}\}_{t\in [\![1,T]\!]}}$ est supposée constante par morceaux avec des discontinuités à chaque vrai instant de ruptures ${\displaystyle \{\tau _{i}^{*}\}_{i\in [\![0,D^{*}]\!]}}$.

Dans le cas de la détection de ruptures dans la distribution, on recode les observations initiales ${\displaystyle \{X_{t}\}_{t\in [\![1,T]\!]}}$ par de nouvelles observations ${\displaystyle \{Y_{t}\}_{t\in [\![1,T]\!]}}$ définies par ${\displaystyle Y_{t}=k_{X_{t}}=k(X_{t},.)}$ où ${\displaystyle k}$ est un noyau symétrique et semi-défini positif (en) (par exemple ${\displaystyle \forall x,y\in \mathbb {R} ,~k^{1}(x,y)=x*y}$ : ${\displaystyle k^{1}}$ est le noyau linéaire ; autre exemple ${\displaystyle \forall x,y\in \mathbb {R} ,k_{\delta }^{2}(x,y)=\exp \left(-{\frac {(x-y)^{2}}{2\delta ^{2}}}\right)}$ : ${\displaystyle k_{\delta }^{2}}$ est le noyau Gaussien de paramètre ${\displaystyle \delta >0}$ ). Pour un noyau symétrique et semi-défini positif ${\displaystyle k}$, le théorème de Moore-Aronszahn assure l'existence d'un espace de Hilbert à noyau reproduisant ${\displaystyle {\mathcal {H}}}$ de noyau reproduisant ${\displaystyle k}$.

Le modèle est :

${\displaystyle {\begin{aligned}\forall t\in [\![1,T]\!],~Y_{t}=\mu _{t}^{*}+\epsilon _{t},\end{aligned}}}$

avec ${\displaystyle \{\mu _{t}^{*}\}_{t\in [\![1,T]\!]}}$ est la fonction de régression et ${\displaystyle ~\{\epsilon _{t}\}_{t\in [\![1,T]\!]}}$ est un bruit d'espérance nulle et de variance ${\displaystyle \sigma ^{2}}$. De plus, ${\displaystyle \forall t\in [\![1,T]\!],~Y_{t},~\mu _{t}^{*},~\epsilon _{t}}$ appartiennent à ${\displaystyle {\mathcal {H}}}$. La fonction de régression ${\displaystyle \{\mu _{t}^{*}\}_{t\in [\![1,T]\!]}}$ est supposée constante par morceaux avec des discontinuités à chaque vrai instant de ruptures ${\displaystyle \{\tau _{i}^{*}\}_{i\in [\![0,D^{*}]\!]}}$.

Méthodes existantes

Le problème de détection de ruptures peut être vu comme un problème de sélection de modèle^[2]: chaque segmentation candidate (i.e. liste d'instants de ruptures candidats) est reliée à un modèle qu'il faut choisir. Nous présentons deux approches : l'une utilisant le principe de la programmation dynamique et l'autre l'heuristique de segmentation binaire dans le cadre classique de la détection de ruptures dans la moyenne puis dans le cas général de la détection de ruptures dans la distribution. Pour chacune de ces deux méthodes, nous présentons les algorithmes permettant de calculer une segmentation en ${\displaystyle D}$ segments ${\displaystyle {\hat {m}}_{D}}$ optimisant un critère (ici, le risque empirique).

Notations

Nous présentons les notations communes aux deux méthodes puis celles qui leur sont spécifiques :

Notation commune

${\displaystyle {\mathcal {M}}(D)}$ l'ensemble des segmentions ${\displaystyle m=\{\{\tau _{0},\dots ,\tau _{1}-1\},\{\tau _{1},\dots ,\tau _{2}-1\},\dots ,\{\tau _{D-1},\dots ,\tau _{D}-1\}\}}$ en ${\displaystyle D}$ segments ${\displaystyle \left({\text{avec la convention}}~\tau _{0}=1,~\tau _{D}=T+1\right)}$.

Notations pour le cas de la détection de ruptures dans la moyenne

• ${\displaystyle X=(X_{1},\dots ,X_{T})}$ et ${\displaystyle s=(s_{1},\dots ,s_{T})}$.
• ${\displaystyle S_{m}^{\mathbb {R} }=\{u\in \mathbb {R} ^{T}:u_{\tau _{0}}=\dots =u_{\tau _{1}-1}\neq u_{\tau _{1}}=\dots =u_{\tau _{2}-1}\neq \dots \neq u_{\tau _{D-1}}=\dots u_{\tau _{D}-1}\}}$ : ${\displaystyle S_{m}^{\mathbb {R} }}$ est l'ensemble des vecteurs ${\displaystyle u\in \mathbb {R} ^{T}}$ constants par morceaux sur les segments de ${\displaystyle m}$.
• Estimateur du risque empirique : ${\displaystyle {\hat {s}}_{m}={\underset {u\in S_{m}^{\mathbb {R} }}{\operatorname {arg\,min} }}\|X-u\|_{2}^{2}}$.

Notations pour le cas de la détection de ruptures dans la distribution

• ${\displaystyle Y=(Y_{1},\dots ,Y_{T})}$ et ${\displaystyle \mu ^{*}=(\mu _{1}^{*},\dots ,\mu _{T}^{*})}$.
• ${\displaystyle Y_{s}^{e}=\{Y_{s},\dots ,Y_{e}\}}$ pour ${\displaystyle s<e}$ avec ${\displaystyle (s,e)\in [\![1,T]\!]^{2}}$.
• ${\displaystyle S_{m}^{\mathcal {H}}=\{u\in {\mathcal {H}}^{T}:u_{\tau _{0}}=\dots =u_{\tau _{1}-1}\neq u_{\tau _{1}}=\dots =u_{\tau _{2}-1}\neq \dots \neq u_{\tau _{D-1}}=\dots u_{\tau _{D}-1}\}}$ : ${\displaystyle S_{m}^{\mathcal {H}}}$ est l'ensemble des vecteurs ${\displaystyle u\in {\mathcal {H}}^{T}}$ constants par morceaux sur les segments de ${\displaystyle m}$.
• Norme dans ${\displaystyle {\mathcal {H}}^{T}}$ : ${\displaystyle \forall h\in {\mathcal {H}}^{T},~\|h\|_{{\mathcal {H}}^{T}}^{2}=\sum _{i=1}^{T}\|h_{i}\|_{\mathcal {H}}^{2}}$.
• Estimateur du risque empirique : ${\displaystyle {\hat {\mu }}_{m}={\underset {u\in S_{m}^{\mathcal {H}}}{\operatorname {arg\,min} }}\|Y-u\|_{{\mathcal {H}}^{T}}^{2}}$.

Les méthodes proposées ci-dessous utilisent le risque empirique comme critère à minimiser (${\displaystyle \|X-{\hat {s}}_{m}\|_{2}^{2}}$ pour la détection de ruptures dans la moyenne ; ${\displaystyle \|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}}$ pour la détection de ruptures dans la distribution). Pour le noyau linéaire ${\displaystyle k^{1}}$, ${\displaystyle \|X-{\hat {s}}_{m}\|_{2}^{2}=\|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}}$, les méthodes utilisées dans le cadre classique se déduisent de celles du cas des noyaux par le biais de ${\displaystyle k^{1}}$ : donc on ne présentera les algorithmes que dans le cas des noyaux.

Programmation dynamique

La méthode de la programmation dynamique utilise le principe d'optimalité de Bellman : toute solution optimale s'appuie elle-même sur des sous-problèmes résolus localement de façon optimale. On utilise cette méthode exacte pour récupérer pour ${\displaystyle D\in [\![1,D_{max}]\!]}$ la meilleure segmentation en ${\displaystyle D}$ segments minimisant le risque empirique i.e. :

${\displaystyle {\begin{aligned}{\hat {m}}_{D}={\underset {m\in {\mathcal {M}}(D)}{\operatorname {arg\,min} }}\|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}.\end{aligned}}}$

Nous présentons dans cette section l'algorithme de programmation dynamique appliquée au problème de détection de ruptures^[3],[4]. Dans un premier temps, nous exprimons le risque empirique ${\displaystyle \|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}}$ en fonction du noyau ${\displaystyle k}$ et de ${\displaystyle \{X_{i}\}_{i\in [\![1,T]\!]}}$ à l'aide des deux résultats suivants ci-dessous.

On montre tout d'abord que, pour ${\displaystyle m\in {\mathcal {M}}(D),~d\in [\![1,D-1]\!]}$, pour ${\displaystyle i\in [\![\tau _{d},\tau _{d+1}-1]\!]}$, ${\displaystyle \left({\hat {\mu }}_{m}\right)_{i}={\frac {1}{\tau _{d+1}-\tau _{d}}}\sum _{j=\tau _{d}}^{\tau _{d+1}-1}Y_{j}}$.

On montre que, pour ${\displaystyle m\in {\mathcal {M}}(D)}$,

${\displaystyle {\begin{aligned}\|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}=\sum _{d\in [\![1,D-1]\!]}C_{\tau _{d},\tau _{d+1}},\end{aligned}}}$

avec, pour ${\displaystyle ~(\tau ,\tau ^{'})\in [\![1,T]\!]^{2},~\tau <\tau ^{'}}$, ${\displaystyle C_{\tau ,\tau ^{'}}=\sum _{i=\tau }^{\tau ^{'}-1}k(X_{i},X_{i})-{\frac {1}{\tau ^{'}-\tau }}\sum _{i=\tau }^{\tau ^{'}-1}\sum _{j=\tau }^{\tau ^{'}-1}k(X_{i},X_{j}).}$

KPGD est l'implémentation du principe de la programmation dynamique à noyau appliquée au problème de détection de ruptures. Elle prend en paramètre la matrice de coût ${\displaystyle C=(C_{i,j})_{(i,j)\in [\![1,T]\!]}(i<j)}$ et elle renvoie ${\displaystyle \{{\hat {m}}_{D}\}_{D\in [\![1,D_{max}]\!]}}$.

 algorithme KPGD (${\displaystyle C}$) 
   for ${\displaystyle d\in [\![2,D_{max}]\!]}$ do
     for ${\displaystyle \tau ^{'}\in [\![d,T]\!]}$ do
       ${\displaystyle L_{d,\tau ^{'}}=\min _{\tau \leq \tau ^{'}}\{L_{d-1,\tau }+C_{\tau ,\tau ^{'}+1}\}}$
       ${\displaystyle m_{d,\tau ^{'}+1}={\underset {\tau \leq \tau ^{'}}{\operatorname {arg\,min} }}\{L_{d-1,\tau }+C_{\tau ,\tau ^{'}+1}\}}$
     end for        
   end for
 Fin algorithme

avec ${\displaystyle \forall \tau >1,L_{1,\tau }=C_{1,\tau }}$.

La sélection de modèle^[4] permet de récupérer ${\displaystyle {\hat {D}}}$ un estimateur du nombre de segments ${\displaystyle D^{*}}$. ${\displaystyle {\hat {D}}}$ est défini par

${\displaystyle {\begin{aligned}{\hat {D}}={\underset {D\in [\![1,D_{max}]\!]}{\operatorname {arg\,min} }}\{\|Y-{\hat {\mu }}_{{\hat {m}}_{D}}\|_{{\mathcal {H}}^{T}}^{2}+pen({\hat {m}}_{D})\},\end{aligned}}}$

avec ${\displaystyle \mathrm {pen} ({\hat {m}}_{D})={\frac {CD}{n}}\left(1+\log \left({\frac {n}{D}}\right)\right)}$. La méthode utilisée pour calibrer la constante ${\displaystyle C}$ est l'heuristique de pente^[5]. On obtient ainsi un estimateur ${\displaystyle {\hat {m}}_{\hat {D}}}$.

Segmentation binaire

L'heuristique de segmentation binaire^[6] est une méthode, fonctionnant par dichotomie, permettant de récupérer un minimiseur local ${\displaystyle {\hat {m}}_{D}}$ du risque empirique ${\displaystyle \|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}}$. Plus précisément, la segmentation binaire cherche à la première itération l'indice ${\displaystyle j_{1}\in [\![2,T-1]\!]}$ de l'instant de ruptures candidat qui minimise le risque empirique ${\displaystyle \|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}}$ : cet indice est l'indice de notre premier instant de ruptures estimé. Puis, elle détermine récursivement, à la deuxième itération, deux instants de ruptures candidats ${\displaystyle j_{2},j_{3}}$ sur chacun des segments délimités par les instants de ruptures estimés. Elle retient comme second instant de ruptures estimé celui (parmi ces deux instants de ruptures candidats) qui minimisent le risque empirique. Puis, on procède de la même manière pour les itérations suivantes. Nous illustrons sur un exemple le fonctionnement de l'algorithme utilisant le principe de la segmentation binaire :

• Étape 1 : A l'itération ${\displaystyle k=1}$, on cherche ${\displaystyle j_{1}}$ qui minimise le risque empirique ${\displaystyle \|Y-{\hat {\mu }}_{m_{1}}\|_{{\mathcal {H}}^{T}}^{2}}$ avec ${\displaystyle m_{1}=\{1,j,T\}}$ pour ${\displaystyle j\in [\![2,T-1]\!]}$ avec :

${\displaystyle {\begin{aligned}\|Y-{\hat {\mu }}_{m_{1}}\|_{{\mathcal {H}}^{T}}^{2}=C_{1,j}+C_{j,T}.\end{aligned}}}$

${\displaystyle j_{1}}$ est notre premier instant de ruptures estimé noté ${\displaystyle {\hat {\tau }}_{1}}$.

• Étape 2  : A l'itération on ${\displaystyle k=2}$, on cherche ${\displaystyle j_{2},~j_{3}}$ minimisant le risque empirique sur ${\displaystyle [\![2,{\hat {\tau }}_{1}-1]\!]}$ et ${\displaystyle [\![{\hat {\tau }}_{1}+1,T-1]\!]}$ respectivement. Par exemple,

${\displaystyle j_{2}={\underset {j\in [\![2,{\hat {\tau }}_{1}-1]\!]}{\operatorname {arg\,min} }}\|Y_{1}^{{\hat {\tau }}_{1}}-{\hat {\mu }}_{m_{2}}\|_{{\mathcal {H}}^{{\hat {\tau }}_{1}}}^{2},}$

avec ${\displaystyle m_{2}=(1,j,{\hat {\tau }}_{1},T)}$. Puis on choisit parmi les instants de ruptures candidats ${\displaystyle j_{2},~j_{3}}$ celui qui minimise le risque empirique ${\displaystyle \|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2}}$
et on le note ${\displaystyle {\hat {\tau }}_{2}}$ . Puis, on continue récursivement.

Ainsi, au bout de ${\displaystyle k}$ itérations, on récupère une segmentation en ${\displaystyle D=k+1}$ segments vérifiant :

${\displaystyle {\begin{aligned}{\hat {m}}_{D}={\underset {m\in {\mathcal {M}}_{D}({\hat {\tau }}_{1},\dots ,{\hat {\tau }}_{D-2})}{\operatorname {arg\,min} }}\|Y-{\hat {\mu }}_{m}\|_{{\mathcal {H}}^{T}}^{2},\end{aligned}}}$

avec ${\displaystyle {\mathcal {M}}_{D}({\hat {\tau }}_{1},\dots ,{\hat {\tau }}_{D-2})}$ l'espace des segmentations en ${\displaystyle D}$ segments où les instants de ruptures estimés ${\displaystyle {\hat {\tau }}_{1},\dots ,{\hat {\tau }}_{D-2}}$ ont été calculés aux itérations précédentes.

Une méthode alternative de segmentation binaire avec temps d'arrêt^[7] permet d'estimer le nombre de segments et donc de récupérer un estimateur ${\displaystyle {\hat {m}}_{\hat {D}}}$ de ${\displaystyle m^{*}}$ .

Notes et références

Articles

• Portail des probabilités et de la statistique