Réduction - Problème inverse : Page pour l'impression

La modélisation physique du phénomène nous a permis d'obtenir un modèle $f(\theta,T)$ où $\theta \in \mathbb{R}^n$ sont les paramètres du modèle et $T \in \mathbb{R}^m$ sont les instants d'observation. De nombreuses sources d'erreurs ont aussi été listées. On se pose maintenant la question suivante: comment estimer les paramètres du modèle compte tenu des observations ?

Après avoir listé diverses sources de signal et de bruits de l'émission de la lumière à la valeur donnée par le déteteur, nous allons maintenant donner une expression finale au modèle. Cette expression n'est pas universelle, et selon l'étoile observée, la précision recherchée, d'autres formulations peuvent être préférables.

Après avoir établi ce modèle, on donnera quelques principes d'analyse statistique et des moyens algorithmiques pour estimer les paramètres du modèle. On verra en particulier que l'analyse "dans le domaine fréquenciel" est particulièrement importante.

Modèle final

Le modèle comportera une partie déterministe et une partie aléatoire. Pour l'astrométrie, la position $x(t) = \alpha(t) \cos \delta(t)$ $y(t)) = \delta(t)$ sur la sphère céleste est

$\begin{array}{ccc} x(t) = x_0 + \mu_x t + \eta_x t^2 + \sum\limits_{k=1}^{n_p} B_j X_j(t) + G_j Y_j(t) + \varphi \Pi_x(t) + S_x(t) + \Delta x_{atm}+ \epsilon_{S_x} + \epsilon_{x,atm} + \epsilon_x\\ y(t) = y_0 + \mu_y t + \eta_y t^2 + \sum\limits_{k=1}^{n_p} A_j X_j(t) + F_j Y_j(t) + \varphi \Pi_y(t) + S_y(t) + \Delta y_{atm} + \epsilon_{S_y}+ \epsilon_{y,atm} + \epsilon_y \end{array}$

Où (x_0,y_0) est la position initiale de l'étoile, $\mu_x, \mu_y$ sont les composantes du mouvement propre $\eta_x, \eta_y$ sont des termes d'accélération de perspective, B_j, G_j, A_j, F_j sont les constantes de Thiele-Innes de la planète , X_j et Y_j sont ses coordonnées sur son plan orbital, $\varphi$ est la parallaxe, $\Pi_x, \Pi_y$ sont les coefficients parametrant le mouvement de la Terre, Les $\epsilon$ sont les bruits résiduels modélisés par des bruits gaussiens.: $\epsilon_{S_x}, \epsilon_{S_y}$ sont les bruits stellaires, $\epsilon_{x,atm}, \epsilon_{y,atm}$ sont les bruits atmosphériques et $\epsilon_x, \epsilon_y$ représentent des bruits instrumentaux.

$V(t) = \dot{\overrightarrow{OT}} + \mu_z + K\left( \cos(v+\omega) + \cos \omega \right) + S_z(t) + z_{atm}(t)+ \epsilon_{S} + \epsilon_{mes}$

Où $\dot{\overrightarrow{OT}$ est la vitesse de l'observateur dans le référentiel barycentrique du système solaire, $\mu_z$ est la composante du mouvement propre dans la direction radiale, S_z(t) est le signal stellaire dû à la granulation, aux oscillations et à l'activité, $\epsilon_S$ est le bruit stellaire résiduel et $\epsilon_{mes}$ le bruit associé à la mesure.

Dans les deux cas, les techniques de réduction de données visent à trouver des paramètres qui sont "plausibles", en l'occurrence, qui reproduisent les observations.

Traitement statistique

Rappelons qu'une variable aléatoire peut être vue comme un programme informatique qui délivre des valeurs suivant une certaine distribution de probabilité lorsqu'on lui demande. Le problème que nous posons maintenant est équivalent au suivant. Supposons qu'un ordinateur ait en mémoire des paramètres $\theta \in \mathbb{R}^p$ (nombre de planètes, leurs caractéristiques orbitales, la période de rotation de l'étoile etc.). S'il n'y avait pas de bruit, une mesure à l'instant t_k reviendrait à demander à l'ordinateur d'évaluer une fonction $f(\theta,t_k)$ . Nous connaissons t_k et (c'est l'un des modèles de la page précédente), mais nous ne connaissons pas $\theta$ . Notre but est de le déterminer à partir des mesures $y_k = f (\theta,t_k)$ . Ce principe est similaire à la résolution d'une énigme: quelqu'un connaît une information $\theta$ et nous essayons de la deviner en posant une question. Dans le cas sans bruit, celui qui pose l'énigme ne nous induit pas en erreur, mais cela ne veut pas dire que la résolution est facile !

Exemple: on veut trouver les paramètres d'une fonction affine du temps f(a,b,t) = at+c (ici $\theta = \left( \begin{array}{c} a \\c \end{array} \right)$ ). On évalue la valeur de en t_1 , on obtient y_1 = at_1+c : on a deux inconnues pour une équation, on ne peut pas résoudre. Si on a y_2 = at_2 + c , avec $t_2 \neq t_1$ , alors on a $a = \frac{y_2-y_1}{t_2-t_1}$ et c = y_1 - at_1 = y_2 - at_2 .

Notre cas est plus compliqué. Les valeurs que nous obtenons sont $y_k = f(\theta, t_k) + b_k$ où b_k est la réalisation d'une variable aléatoire $\epsilon_k$ . A l'appel numéro du programme, l'ordinateur fait appel à un autre programme qui délivre une variable alétoire selon une certaine loi. En l'occurrence, nous supposons cette loi gaussienne. Si nous faisons mesures, tout se passe comme si un programme principal évaluait la fonction $f(\theta, t_k)$ et programmes secondaires $\epsilon_k$ , k=1..m retournent chacun une valeur b_k . Si nous pouvions remonter le temps et faire les mesures plusieurs fois aux mêmes instants, on aurait des des vecteurs de mesures $y_1 = \left( \begin{array}{c} y_1 \\y_2 \\ \\y_m \end{array} \right) = f(\theta,T) + b_1$ , puis $y_2 = f(\theta,T) + b_2$ etc. (notez qu'ici b_1 et b_2 sont des vecteurs.

Comme nous ne connaissons que la loi suivie par les variables $\epsilon_k$ , ils nous est impossible de connaître les paramètres avec certitude. On leur attache une "erreur", qui quantifie l'incertitude que l'on a sur eux. En reprenant l'exemple précédent on mesure y_1 = at_1 + c +b_1 et y_2 = at_2 + c +b_2 . Si on estime et avec les mêmes formules, on fera une erreur $\frac{\sqrt{\sigma_1^2 + \sigma_2^2}}{t_1-t_2}$ sur et une erreur $\sqrt{\left(\frac{t_1}{t_2-t_1} \sigma_2\left)^2 + \right((1-\frac{t_1}{t_2-t_1}) \sigma_1 \left)^2}$ sur (admis).

Bruits

Dans le modèle $f(\theta,T) + \epsilon$ , le symbole $\epsilon$ désigne un bruit gaussien. Comme ils apparaissent constamment en détection de planètes extrasolaires et ailleurs, nous allons en donner quelques propriétés.

A une expérience donnée, $\epsilon$ prendra une valeur imprévisible. La probabilité que la valeur de soit comprise entre b_1 et b_2 est $\int_{b_1}^{b_2} f(x) dx$ où f(x) est la densité de probabilité de $\epsilon$ . Dire que $\epsilon$ est un bruit gaussien veut dire que sa densité est de la forme $f(x) = \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{(x-\mu)^2}{\sigma^2}}$ où $\mu$ et $\sigma$ sont des réels, qui sont égaux respectivement à la moyenne et à l'écart-type de $\epsilon$ . On note souvent $\epsilon \sim g(\mu, \sigma^2)$ , qui signifie " $\epsilon$ suit une loi gaussienne de moyenne $\mu$ et de variance $\sigma^2$ . Dans la plupart des cas, le bruit est de moyenne nulle (c'est le cas ici).

Dans le modèle, des bruits d'origines différentes s'additionnent. Sachant que le résidu de l'activité stellaire que nous n'avons pas ajusté $\epsilon_S$ et le bruit de mesure $\epsilon_{mes}$ suiven une certaine loi, quelle loi suivra $\epsilon = \epsilon_{S} + \epsilon_{mes}$ ? Nous pouvons déjà dire que la moyenne de $\epsilon$ sera égale à la somme des moyennes de $\epsilon_S$ et $\epsilon_{mes}$ car l'espérance est un opérateur linéaire. Peut-on dire plus ? Si ces bruits dépendaient l'un de l'autre, la réponse pourrait être complexe. En l'occurrence, la physique de l'étoile cible et les erreurs instrumentales sont totalement indépendantes. On peut montrer que dans ces conditions, la variance de $\epsilon$ est égale à la somme des variances de $\epsilon_S$ et $\epsilon_{mes}$ . Nous pouvons même aller plus loin car la somme de deux variables gaussiennes indépendante est une variable gaussienne. En résumé, $\epsilon \sim g\left(\mu_{S} +\mu_{mes}, \sigma_{S}^2+ \sigma_{mes}^2\right)$ en l'occurrence $\mu_S$ et $\mu_{mes}$ sont nulles.

Lorsqu'on dispose de plusieurs mesures, à l'expérience numéro on a un certain bruit b_k réalisation d'une variable $\epsilon_k$ de densité f_k . La plupart du temps, on fait l'hypothèse que les brutis $\epsilon_k$ sont indépendants, c'est à dire que la probabilité d'obtenir le bruit b_k à l'expérience ne dépend pas des valeurs prises aux expériences précédentes et suivantes. Lorsque ce n'est pas le cas on parle de bruits corrélés. Pour les caractériser, on utilise souvent leur densité spectrale de puissance. Un certain profil de densité spectrale correspond à une "couleur" du bruit.

A retenir: la somme de variables gaussienne indépendantes $\epsilon = \sum\limits_{k=1}^n \epsilon_k$ où $\epsilon_k \sim g\left( \mu_k, \sigma_k^2)$ est une variable gaussienne suivant la loi $g \left( \sum\limits_{k=1}^n \mu_k , \frac{1}{n} \sum\limits_{k=1}^n \sigma_k^2 \right)$ .

Signification statistique

Il esiste plusieurs outils pour s'assurer qu'une détection possible n'est pas due au bruit. L'un des plus utilisé est le test de signification (significance en anglais), qui consiste à calculer la probabilité d'avoir le signal observé "au moins aussi grand" s'il n'y avait en réalité que du bruit. Par exemple, qupposons que l'on veuille mesurer une quantité qui est perturbée par un bruit gaussien additif de moyenne nulle et d'écart-type $\sigma=1$ , donnant une mesure y = a+b .

Si il n'y avait en réalité pas de signal ( a=0 ), les mesures seraient uniquement dues au bruit. On imagine deux cas de figures:

1er cas: on mesure . La probabilité d'avoir obtenue une valeur au moins aussi grande par hasard, c'est à dire la probabilité que $|y| \geq 0.1$ sachant que est $\int_{-\infty}^{-0.1} e^{-x^2}dx + \int_{0.1}^{\infty} e^{-x^2}dx = 2 \int_{0.1}^{\infty} e^{-x^2}dx \approx 0.818$
2ème cas: on mesure . La probabilité pour que alors que est $\int_{-\infty}^{-10} e^{-x^2}dx + \int_{10}^{\infty} e^{-x^2}dx = 2 \int_{10}^{\infty} e^{-x^2}dx \approx 2.06 \times 10^{-9}$

Ces valeurs ont une interprétation: si on réalisait exactement le même type de mesure alors qu'il n'y a pas de signal, on observerait $|y| \geq 0.1$ dans 81.8 % des cas et $|y| \geq 10$ dans $2.06 \times 10^{-7}$ % des cas.Dans le premier cas, la probabilité d'avoir un signal aussi grand que celui que l'on a mesuré est grande. On ne peut pas assurer qu'un signal a été détecté. Par contre, dans le deuxième cas on serait dans un des deux cas sur un milliard où le signal serait dû au bruit. On peut alors dire qu'on a détecté un signal à $10 \;\sigma$ , car la valeur de l'écart-type du bruit est , sa moyenne est , donc on a $\frac{y}{\sigma-0} = \frac{10}{1} = 10$ .

Détecter un signal "à 10 sigmas" est un luxe que l'on peut rarement se payer. Les détections sont annoncées plutôt pour des valeurs de sigmas.

Remarque importante: on calcule la probabilité d'avoir les observations sachant qu'il n'y a pas de signal et non la probabilité d'avoir un signal sachant les observations qui est une quantité qui a davantage de sens. Le calcul de cette dernière quantité se fait dans le cadre du calcul bayésien, outil très puissant qui ne sera pas développé dans ce cours.

Exercices

Auteur: Nathan Hara

Espérance et variance de la moyenne empirique

Difficulté : ☆☆☆

Question 1)

On rappelle que l'espérance et la variance d'une variable alétoire de densité de probabilité sont données par $\mathbb{E}(X) = \int_{-\infty}^{+\infty} xf(x)dx$ et $\text{Var}(X) = \mathbb{E}\left((X-\mathbb{E}(X))^2\right)= \int_{-\infty}^{+\infty} (x-\mathbb{E}(X))^2f(x)dx$

Soit une variable aléatoire et un réel. Montrer que $\mathbb{E}(aX) = a \mathbb{E}(X)$ et $\text{Var}(aX) = a^2 \text{Var}(X)$
Soient et deux variables aléatoires indépendantes. Montrer que $\mathbb{E}(X+Y) = \mathbb{E}(X) + \mathbb{E}(Y)$ .
On définit la covariance par $\text{Cov}(X,Y) = \mathbb{E}\{ (X-\mathbb{E}(X)) (Y-\mathbb{E}(Y))\}$ . Montrer que $\text{Var}(X+Y) = \text{Var}(X) + \text{Var}(Y) + 2 \text{Cov}(X,Y)$
Nous allons maintenant voir un cas où la précision d'un estimateur est facilement calculable. On considère mesures entachées de bruits gaussiens d'une quantité fixe . Plus précisément, la mesure numéro est modélisée par une variable aléatoire $X_j = K + \epsilon_j$ où $\epsilon_j$ est un bruit gaussien de moyenne nulle et de variance $\sigma^2$ . On suppose que les mesures sont indépendantes, ce qui implique en particulier que $\text{Cov}(\epsilon_j,\epsilon_i) = 0$ pour $i \neq j$ . Pour obtenir une estimation de , on fait la moyenne empirique des expériences, c'est à dire $M = \frac{X_1 + X_2 +... X_m}{m}$ . Montrer que $\mathbb{E}(M) = K$ et $\text{Var}(M) = \frac{\sigma^2}{m}$
La précision de l'estimateur augmente-t-elle avec le nombre de mesures ?

Loi du chi 2

Nous nous sommes toujours ramenés à des modèles du type: modèle déterministe + bruit gaussien. Afin de vérifier que les observations sont compatibles avec le modèle, on étudie les résidus, définis comme "les observation - le modèle ajusté".

La loi du $\chi^2$ est un outil commode pour étudier le comportement de plusieurs variables gaussiennes. Considérons d'abord une famille de variables aléatoires gaussiennes indépendantes $(\epsilon_k)_{k=1..m}$ , de moyenne nulle et de variance unité. On forme la quantité $\chi^2_m= \epsilon_1^2 + \epsilon_2^2 + ... + \epsilon_m^2$ . Comme les $\epsilon_k$ sont des variables aléatoires, les $\epsilon_k^2$ le sont aussi. La somme de variables alétoires étant toujours une variable alétoire, $\chi^2_m$ suit une certaine loi de probabilité. Dans l'analogie avec un programme informatique, la variable alétoire $\chi^2_m$ se comporte comme un programme qui appelle programmes générant une variables gaussiennes, puis additionne leurs carrés. Elle est appelée loi du $\chi^2$ à m degrés de liberté. On peut montrer qu'en moyenne une variable gaussienne au carré $\epsilon^2$ a une moyenne de 1. En conséquence, $\chi^2_m$ vaudra typiquement .

Pourquoi cette loi serait utile pour notre cas ? Si le modèle est bien ajusté, les résidus doivent se comporter comme un bruit gaussien. En supposant que les erreurs sont toutes indépendantes, de moyenne nulle et de variance unité, les résidus $r_k = y-f(\theta,t_k})$ en sont une réalisation. Donc $R^2=\sum\limits_{k=1}^m r_k^2$ est une réalisation d'une loi du $\chi^2$ à degrés de liberté. Si R^2 est de l'ordre de , le modèle est cohérent. Sinon, le modèle ou les paramètres ajustés sont à revoir.

En pratique, les erreurs $\epsilon__k$ ne sont évidemment pas de variance unité et parfois pas indépendantes. Par contre on peut à bon droit supposer qu'elles sont de moyenne nulle. Pour se ramener au cas précédent, on calcule non pas une réalisation de $\epsilon_1^2+...+\epsilon_m^2$ mais de $\chi^2_m'=\epsilon^T \Sigma^{-1} \epsilon$ où $\epsilon = \left( \begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \\ \epsilon_m \end{array} \right)$ , $\epsilon^T$ est sa transposée et $\Sigma$ est la matrice des variances-covariances de $\epsilon$ . Dans le cas où les $\epsilon_k$ sont indépendantes, la matrice des variances-covariances est diagonale, son -ème terme diagonal étant $\frac{1}{\sigma_k^2}$ , soit l'inverse de la variance de $\epsilon_k$ .

La vraisemblance

Etant donné des paramètres $\theta$ , le modèle global $f(\theta,T)+\epsilon$ est une variable aléatoire: il est somme d'une variable aléatoire valant $f(\theta,T)$ avec une probabilité 1 et d'un vecteur de variables aléatoires gaussienne $\epsilon$ . A ce titre, il a une certaine densité de probabilité que l'on note $L(y|\theta)$ . Le symbole | se lisant "sachant". La lettre L vient de Likelihood, qui veut dire vraisemblance en anglais. Il s'agit dans l'idée de la probabilité d'obtenir $y = f(\theta,T)+\epsilon$ pour une valeur de $\theta$ donnée..

La fonction $\theta \rightarrow L(y|\theta)$ est souvent appelée "fonction de vraisemblance". La valeur de $\theta$ maximisant $L(y|\theta)$ est appelé l'estimateur du maximum de vraisemblance. Il a de bonnes propritétés statistiques. En effet, on peut montrer que c'est un estimateur:

Non biaisé: lorsque le nombre de mesure tend vers l'infini, l'estimateur du maximum de vraisemblance tend vers la bonne valeur $\theta$
De variance minimale: parmi les estimateurs non biaisés, il a la variance la plus faible que l'on puisse espérer.

Dans notre cas, si les $\epsilon_k$ sont des variables indépendantes, leur densité de probabilité jointe est égale au produit de leurs densité de probabilité. $L(y|\theta) = \prod\limits_{k=1}^m g_k(y-f(\theta,t_k))$ où g_k est la densité de probabilité de la variable $\epsilon_k$ . De plus, si ces varibles sont gaussiennes et indépendantes, on a: $L(y|\theta) = \prod\limits_{k=1}^m \frac{1}{\sqrt{2\pi}\sigma_k} e^{-\frac{(y-f(\theta,t_k))^2}{\sigma_k^2}} = \frac{1}{(\sqrt{2\pi})^m \prod\limits_{k=1}^m \sigma_k} e^{-\sum\limits_{k=1}^m \frac{(y-f(\theta,t_k))^2}{\sigma_k^2}}}$

Méthode des moindres carrés

Pour des bruits gaussiens indépendants, maximiser la vraisemblance, revient à minimiser $\theta \rightarrow F(y,\theta)=-ln(L(y|\theta)) =\sum\limits_{k=1}^m \frac{(y-f(\theta,t_k))^2}{\sigma_k^2}}$ puisque $x \rightarrow e^{-x}$ est une fonction décroissante. La méthode consistant à minimiser $F(y,\theta)$ s'appelle la méthode des moindres carrés. C'est la méthode d'estimation de loin la plus utilisée dans tous les domaines. Elle est parfois utilisée quand les bruits ne sont pas gaussiens, mais il faut garder à l'esprit qu'elle n'a alors plus de propriétés statistiques sympathiques (sauf quand le modèle est linéaire en $\theta)$ .

Dans notre cas, les paramètres $\theta$ sont les éléments des orbites, les paramètres du bruit stellaire, du mouvement propre, etc. La fonction $F(y,\theta)$ a donc de nombreux paramètres, et trouver son minimum global est une tâche ardue qui fait l'objet d'une littérature très vaste.

Lorsque le modèle est linéaire en $\theta$ i. e. $f(T,\theta) = A \theta$ où est une certaine matrice dépendant des instants d'observation $T = (t_k)_{k=1..m}$ , l'ajustement est beaucoup plus simple car il a une solution explicite (voir mini-projet). On essaye de se ramener autant que possible à des ajustements linéaires. La plupart du temps, on estime les paramètre les uns après les autres, puis un ajustement global est réalisé. Une démarche classique consiste à:

Ajuster le mouvement propre et la parallaxe, puis soustraire ces ajustements du signal initial
Chercher les périodes éventuelles dans le signal à l'aide du périogogramme (voir page suivante)
Si on trouve une période, on ajuste une orbite, que l'on soustrait du signal
On recommence les étapes 2 et 3 jusqu'à ce que le périodogramme ne présente plus de pic significatif

Sur la figure, on représente les étapes d'un ajustement d'un signal astrométrique simulé (2 planètes, 45 observations). De gauche à droite et de haut en bas:

Signal brut en rouge En vert signal après correction des bruits instrumentaux et atmosphériques . Le mouvement propre (tendance linéaire) et la parallaxe (oscillation) dominent très nettement le signal.
En vert: le signal planétaire réel, sans aucun bruit $(s^{\star})$ . En rouge: signal estimé en soustrayant le mouvement propre, l'accélération de perspective et la parallaxe au signal .
En rouge: signal . En vert : orbite d'une planète ajustée sur .
En rouge : signal dont on a soustrait . En vert :orbite d'une planète ajustée sur .
En rouge: signal . En vert: , ajustement de deux planètes en partant des estimations et .
En vert, signal $(s^{\star})$ . On ajuste tous les paramètres simultanément sur en partant des estimations précédentes. Le signal en rouge représente le signal planétaire de cet ajustement.

Exemple d'ajustement en astrométrie (données simulées)

Six étapes d'ajustement aux données. Pour représenter la chronologie des mesures, on trace des segments reliant deux points de mesure consécutifs.

Périodogramme

Evaluer les résidus sur une grille d'un modèle à paramètres, où chacun d'eux peut prendre valeurs recquiert n^p évaluations, ce qui devient rapidement ingérable numériquement. Les planètes ont un mouvement périodique, donc il est raisonnable de checher des signaux périodiques dans le signal en ne faisant varier que la période du signal recherché. Pour des signaux échantillonnés à intervalles réguliers, on utilise la transformée de Fourier. Le périodogramme est un moyen de checher des signaux périodiques dans des données échantillonnées irrégulièrement. On les notera . Le périodogramme de Lomb-Scargle d'un signal $(y_k)_{k=1..m}$ échantillonné aux instants $(t_j)_{j=1..m}$ est défini comme suit pour une fréquence quelconque $\omega$ :

$P_y(\omega) = \frac{1}{2} \left( \frac{ \left( \sum\limits_{j=1}^m y_j \cos \omega(t_j-\tau) \right)^2}{ \sum\limits_{j=1}^m \cos^2 \omega(t_j-\tau)}} + \frac{ \left( \sum\limits_{j=1}^m y_j \sin \omega(t_j-\tau) \right)^2}{ \sum\limits_{j=1}^m \sin^2 \omega(t_j-\tau)}} \right)$ où $\tau$ vérifie:

$\tan 2 \omega \tau = \frac{\sum\limits_{j=1}^m \sin 2\omega t_j}{\sum\limits_{j=1}^m \cos 2\omega t_j}$

Cette expression est équivalente à $P_y(\omega) = A^2 +B^2$ où et sont les paramètres minimisant $\sum\limits_{j=1}^m (y_j-A\cos(\omega t_j) - B\sin(\omega t_j))^2$ . Le modèle $A,B, t \rightarrow A \cos \omega t + B \sin \omega t$ est linéaire en A, B , on a donc une solution explicite à la minimisation.

Le périodogramme a une propriété très intéressante: si le signal d'entrée est un bruit gaussien $\epsilon$ de variance unité, p_0 une valeur réelle fixée et $\omega$ une fréquence quelconque, la probabilité que $P_{\epsilon}(\omega)$ dépasse p_0 est $Pr\{P_\epsilon(\omega)>p_0\} = e^{-p_0}$ . En d'autres termes, la probabilité qu'une valeur du périodogramme à $\omega$ fixée soit "au moins aussi grand que p_0 " par hasard décroît exponentiellement. Supposons que l'on ait un signal y(t) = x(t) + b(t) où b(t) est un bruit gaussien de variance unité et nous trouvons un pic de taille p_0 , on calcule la probabilité de trouver un pic au moins aussi grand si le signal n'est composé que de bruit: $e^{-p_0}$ . Si cette valeur est petite, on pourra confirmer la détection d'un signal avec une erreur de fausse alarme de $e^{-p_0}$ . Ce procédé n'est autre qu'un test de signification statistique.

La figure montre un exemple de périodogramme. Il s'agit d'un périodogramme d'une des coordonnées d'un signal astrométrique simulé dont on a soustrait le mouvement propre et la parallaxe. En bleu, on représente un périodogramme idéal, sans bruit, avec 10000 observations. Le périodogramme représenté en rouge est lui calculé pour 45 observations. Le pic le plus haut correspond bien à une fréquence réelle. Par contre, le deuxième pic le plus important (à 0.37 rad/s) ne correspond pas à une sinusoïde. C'est ce qu'on appelle un alias de la fréquence principale.

En pratique, la variance du bruit n'est pas unitaire et dépend de l'instant de mesure. On peut corriger ce problème en minimisant un critère pondéré $\sum\limits_{j=1}^m \frac{(y_j-A\cos(\omega t_j) - B\sin(\omega t_j))^2}{\sigma_j^2}$ où $\sigma_j$ est la variance du bruit à la mesure .

Exemple de périodogramme

En rouge: périodogramme d'un système simulé à deux planètes après soustraction du mouvement propre et de la parallaxe (signal très peu bruité), pour 45 mesures. En bleu: périodogramme du même système sans bruit et avec 10000 observations.

Réduction - Problème inverse

Objectif