Méthodes numériques pour l’astrophysique M2 « Astrophysique et Milieux Dilués »

1. Méthodes numériques pour l’astrophysiqueM2 « Astrophysique et Milieux Dilués » Hervé Beust Laboratoire d’Astrophysique de Grenoble Master 2 AMD - Méthodes numériques pour l'astrophysique

2. Méthodes numériques pour l’astrophysique • Techniques de base • Estimateurs et statistique • Modélisation de données • Résolution numérique d’équations différentielles • Equations aux dérivées partielles Master 2 AMD - Méthodes numériques pour l'astrophysique

3. Méthodes numériques pour l’astrophysique • Techniques de base • Résolution numérique d’équations • Intégration numérique • Minimisation de fonctions • Estimateurs et statistique • Modélisation de données • Résolution numérique d’équations différentielles • Equations aux dérivées partielles Master 2 AMD - Méthodes numériques pour l'astrophysique

4. Résolution numérique d’équations • Problème : résoudre numériquement une équation de type f(x) = 0, en général par une méthode itérative. • Le problème est très différent suivant que f est une fonction scalaire ou vectorielle. • En général, les méthodes classiques fonctionnent assez bien pourvu que • On soit sûr qu’il y ait une racine; • On parte du voisinage de la racine. • Ca marche mieux si on encadre la racine entre deux réels a et b et si la fonction est monotone dans l’intervalle. Master 2 AMD - Méthodes numériques pour l'astrophysique

5. Différents cas de figure... Cas standard  Cas plus délicats  Master 2 AMD - Méthodes numériques pour l'astrophysique

6. Méthode de dichotomie • On prend c=(a+b)/2 et on calcule f(c). Si f(c) est du même signe que f(a), la racine est entre c et b; sinon entre a et c. On recommence avec le nouvel intervalle.Do while (abs(a-b)>…) c = (a+b)/2 if f(c)*f(a)>0 then a = c else b = c end if end do • Cette méthode a l’avantage de marcher tout le temps pourvu qu’il y ait une racine. Elle ne nécessite pas le calcul de la dérivée. • C’est une méthode d’ordre 1 a convergence moyennement rapide. Master 2 AMD - Méthodes numériques pour l'astrophysique

7. Autres méthodes d’ordre 1 • Méthodes de la sécante et de la fausse position • Ces méthodes convergent en général plus vite que la dichotomie, mais pas toujours… Master 2 AMD - Méthodes numériques pour l'astrophysique

8. Méthode de Newton-Raphson • On part d’un point, on approxime la fonction par sa tangente, on recommence avec la racine trouvée, etc… f(x+h) ≈ f(x)+h f ’(x)+(h2/2) f ’’(x)+..A xn+1 = xn - f(xn)/f ’(xn) • Cette méthode ne nécessite pas d’encadrer la racine au préalable, mais juste de partir d’un point. Elle fonctionne dans le complexe. • C’est une méthode d’ordre 2 très efficace. • MAIS, il faut être sûr d’être au voisinage de la racine. La méthode peut échouer sinon. Master 2 AMD - Méthodes numériques pour l'astrophysique

9. Méthode de Newton-Raphson : Problèmes • Il faut être sûr d’être au voisinage d’une racine… • Il faut aussi pouvoir calculer la dérivée. Sinon:f’(x) ≈ (f(x+h)-f(x-h)) / 2h Master 2 AMD - Méthodes numériques pour l'astrophysique

10. Amélioration : Newton-Raphson quartique SUBROUTINE KEPLER0(M,E,U) IMPLICIT NONE REAL*8 M, ! Mean anomaly & E, ! Eccentricity & U, ! Eccentric anomaly & F,F1,F2,F3, ! Intermediaires & DELA ! Correction REAL*8, PARAMETER :: EPSI=1d-12 LOGICAL OK INTEGER I U = M OK = .FALSE. I=0DO WHILE (.NOT.OK)F = U-E*SIN(U)-MF1 = 1.d0-E*COS(U) F2 = U-M-FF3 = 1.d0-F1 DELA = -F/F1 DELA = -F/(F1+0.5d0*DELA*F2) DELA = -F/(F1+0.5d0*DELA*F2+DELA*DELA*F3/6.d0)U = U+DELA I = I+1OK = ((ABS(DELA).LT.EPSI).OR.(I.GE.500))END DOEND • On peut améliorer la convergence si on sait calculer des dérivées d’ordre supérieur • Exemple : équation de KéplerM = u - e sin u Master 2 AMD - Méthodes numériques pour l'astrophysique

11. En multidimensions … • En dimensions plus que 1, c’est beaucoup plus difficile. • La seule méthode générale, c’est Newton-Raphson, sans garantie de succès… • Jxn = Matrice jacobienne en x • A chaque étape il faut résoudre un système linéaire • La convergence est quadratique pouvu qu’on soit au voisinage de la racine… Master 2 AMD - Méthodes numériques pour l'astrophysique

12. Intégration numérique • Problème : calculer numériquement une intégrale de la forme • w(x)est une fonction de poids • La technique consiste toujours à faire une combinaison linéaire d’évaluations de la fonction f en un certain nombre de pointsxi , i=1..n • Tout repose dans le choix des coefficientsaiet des pointsxi. • Les formules convergent vers l’intégrale lorsque n ∞ Master 2 AMD - Méthodes numériques pour l'astrophysique

13. Formules à points régulièrement espacés • Formule des trapèzes : • Formule de Simpson : Master 2 AMD - Méthodes numériques pour l'astrophysique

14. Quadratures de Gauss – polynômes orthogonaux • Idée : ne plus imposer que les xi soient régulièrement espacés. • Choisir les xiet les ai de manière à ce que la formule d’approximation soit exacte pour des polynômes de degré leplus élevé possible. • C’est possible jusqu’à des polynômes de degré 2n-1 : Pour un choix donné de a, b, w(x), et n il existe une unique jeu de xiet les ai telle que l’approximation soit exacte pour toute fonction f polynôme de degré ≤ 2n-1. • La théorie est liée à celle des polynômes orthogonaux. On définit • La suite des polynômes Pn lorsque n varie est orthogonale au sens du produit scalaire suivant Master 2 AMD - Méthodes numériques pour l'astrophysique

15. Quadratures de Gauss – polynômes orthogonaux • Les Pnobéissent à la récurrence suivante : • Les xise calculent comme les racines de Pn et les aivérifient : Master 2 AMD - Méthodes numériques pour l'astrophysique

16. Quadratures de Gauss – polynômes orthogonaux Cas particuliers : Dans le cas de Gauss-Legendre, on a en outre Et pour tout intervalle (a,b) on a Master 2 AMD - Méthodes numériques pour l'astrophysique

17. Exemple résolu : Gauss-Legendre SUBROUTINE GAULEG(X1,X2,X,W,N) IMPLICIT NONE INTEGER*4 N,M, ! Ordres & I,J ! Indices REAL*8 X(N), ! Tableau de racines & W(N), ! Tableau de poids & Z, ! Racine & X1,X2,XM,XL, ! Bornes & Z1,P1,P2,P3, ! Intermediaires & PP ! Derivee LOGICAL OK ! Test d'arret REAL*8, PARAMETER :: EPS=3.d-14, & PI=3.14159265358979323846d0 M = (N+1)/2 XM = 0.5d0*(X2+X1) XL = 0.5d0*(X2-X1) ! Conversion d’intervalle DO I = 1,M Z = COS(PI*(I-0.25d0)/(N+0.5d0)) ! Guess initial OK = .FALSE. DO WHILE(.NOT.OK) P1 = 1.0d0 P2 = 0.0d0 DO J = 1,N P3 = P2 P2 = P1 P1 = ((2.0d0*J-1.0d0)*Z*P2-(J-1.0d0)*P3)/J END DO ! Calcul de Pn(Z) par récurrence PP = N*(Z*P1-P2)/(Z*Z-1.0d0) Z1 = Z Z = Z1-P1/PP ! Newton-Raphson OK = (ABS(Z-Z1).LT.EPS) END DO X(I) = XM-XL*Z X(N+1-I) = XM+XL*Z W(I) = 2.d0*XL/((1.0d0-Z*Z)*PP*PP) W(N+1-I) = W(I) END DO END Master 2 AMD - Méthodes numériques pour l'astrophysique

18. Intégration multidimensionnelle • C’est en général beaucoup plus difficile. • Technique de base : intégrales emboîtées • 2 difficultés : • Nombre de points nd • Ca ne marche que si le domaine à intégrer n’est pas trop compliqué • Dans les autres cas : Intégration Monte-Carlo Master 2 AMD - Méthodes numériques pour l'astrophysique

19. Intégration Monte-Carlo • On veut intégrer une fonction f sur un domaine W. La moyenne de f sur W vaut par définition • On tire au hasard un nombre n de points dans l’ensemble et on approxime la moyenne par • L’erreur moyenne sur l’intégrale esta convergence en n-1/2, lente • Si W n’est pas calculable directement, on inclut le volume W dans un volume V plus grand, facilement calculable, et on prolonge f par 0 en dehors de W. On intègre ensuite sur V. Mais il faut choisir un volume V proche de W sinon, la convergence est encore plus lente. Master 2 AMD - Méthodes numériques pour l'astrophysique

20. Minimisation de fonctions • Problème : Trouver numériquement un point minimum (ou maximum) d’une fonction f donnée • Un minimum peut être local ou global : les méthodes itératives garantissent en général la convergence vers un minimum local. • Le problème est très différent suivant que la fonction est mono- ou multidimensionnelle • Il y a deux familles de méthodes : les méthodes avec gradient et celles sans gradient. Master 2 AMD - Méthodes numériques pour l'astrophysique

21. Méthodes en dimension 1 • On veut minimiser une fonction f d’une variable réelle. • Avant de converger, il faut encadrer un minimum : trouver 3 points a<b<c, tels que f(b)<f(a) et f(b)<f(c). • Pour encadrer, on déplace le triplet dans le sens de la descente jusqu’à trouver un encadrement. • Première méthode simple une fois qu’on a un encadrement: On tire un réel d entre a et c, on évalue f(d), et on cherche un nouvel intervalle d’encadrement… puis on recommence… • Ca marche toujours, mais ce n’est pas très rapide. Master 2 AMD - Méthodes numériques pour l'astrophysique

22. Méthode : interpolation par paraboles inverses • Idée : près d’un minimum, la fonction s’approxime bien par une parabole. • Quand on a un encadrement (a,b,c), on calcule la fonction de degré 2 (parabole) qui passe par ces trois points et on en cherche le minimum d. On a • Le nouvel encadrement est (a,d,b) ou (b,d,c). • La méthode de Brent combine les deux méthodes précédentes. Master 2 AMD - Méthodes numériques pour l'astrophysique

23. Méthodes multidimensionnelles • On cherche à minimiser une fonction f(x1,…,xn) = f(x) • La plupart des méthodes procèdent à des minimisations successives le long de plusieurs directions. • Minimisation le long d’une direction: étant donné deux vecteurs x et n, trouver le réel l qui minimise f(x+ln).aOn utilise une méthode de dimension 1 • A partir de x+ln, on recommence à minimiser dans une autre direction n’. • Problème : En minimisant le long de n’, il ne faut pas détruire la minimisation le long de n.aLes directions n et n’ doivent être conjuguées. Master 2 AMD - Méthodes numériques pour l'astrophysique

24. Directions conguguées • On se place en un point P • Le changement de gradient dans un déplacement infinitésimal vaut • Une fois qu’on a minimisé le long de u, si on recommence le long de v, il faut que le gradient reste perpendiculaire à u. • Les directions u et v sont alors dites conjuguées. Master 2 AMD - Méthodes numériques pour l'astrophysique

25. Méthode de Powell • Idée : faire des minimisations successives le long de directions conjuguées. • Problème : Comment trouver des directions conjuguées en chaque point ? • La méthode de Powell minimise le long de directions qui tendent à être conjuguées au fil des itérations : On a à chaque instant n directions ui, i=1..n. Au départ, on prend ui=ei(vecteurs de base). On répète ensuite la séquence suivante: • On part d’une position de départ x0. • On fait n minimisation successives le long des n directions. On appelle x1,… xn les points trouvés. • On remplace u0 par u1, puis u1 par u2, …, un-1 par un. • On remplace un par xn-x0 • On minimise dans la direction de un à partir de xn et on remplace x0 par le point trouvé. • La méthode a un défaut : elle tend à produire des directions ui linéairement dépendantes • Solution : Réinitialiser les directions (ui=ei) toutes les n itérations. • Remarque :Cette méthode est une méthode sans gradient. Master 2 AMD - Méthodes numériques pour l'astrophysique

26. Méthodes avec gradient • Quand on sait calculer le gradient f en tout point… • Méthode de descente maximale : On part de x, on minimise dans la direction de -f , et on recommence. • Cette méthode ne donne pas de bons résultats car les directions successives ne sont pas conjuguées. Master 2 AMD - Méthodes numériques pour l'astrophysique

27. Méthode de Fletcher-Reeves-Polak-Ribiere • On part d’un point x0, on fixe une direction u0, et on calcule • A chaque étape on dispose de • On minimise le long de ui. On appelle xi+1 le nouveau point trouvé. On calcule • On calcule ui+1 comme • On recommence. • Cette méthode garantit que les directions ui sont conjuguées. Elle est recommandée quand on n’a aucune idée de la matrice Hessienne. Master 2 AMD - Méthodes numériques pour l'astrophysique

28. Méthodes numériques pour l’astrophysique • Techniques de base • Estimateurs et statistique • Propriétés et exemples • Convergence • Recherche d’estimateurs • Estimation par ajustement • Modélisation de données • Résolution numérique d’équations différentielles • Equations aux dérivées partielles Master 2 AMD - Méthodes numériques pour l'astrophysique

29. Estimateurs et statistiques • Une série de mesures doit être considérée comme une séries de réalisations d’une ou plusieurs variables aléatoires • Le phénomène aléatoire dépend d’un paramètre physique θque l’on cherche à estimer. • Un estimateur est une fonction T arbitraire d’un échantillon de données(X1,…,Xn) censée représenter (estimer) le paramètre inconnu θ • En général, on a plutôt une suite d’estimateurs(Tn) de θen fonction de la taille n de l’échantillon • Exemple : On cherche à estimer la moyenne (espérance) d’une variable aléatoire. Si on fait n tirages successifs, un estimateur de la moyenne est • On dit que l’estimateur est convergent si pour tout ε>0, on a • La loi faible des grands nombre montre que Xn est un estimateur convergent de l’espérance. Master 2 AMD - Méthodes numériques pour l'astrophysique

30. Propriétés des estimateurs • Si (Tn) est un estimateur convergent du paramètre θ, si φ est une fonction de IR dans IR, continue en θ, alors φ(Tn) est un estimateur convergent de φ(θ). • Exemple :X suit la loi uniforme sur [0,θ].Xn est un estimateur convergent de la moyenne = θ/2 Tn=2Xn est un estimateur convergent de θ. • On pose Y=ln X. On a alors E(Y)=ln θ – 1 (calcul). Doncest un estimateur convergent de ln θ – 1 . Et ensuiteest un autre estimateur convergent de θ • On peut en créer beaucoup d’autres. Comment sélectionner le meilleur ? Master 2 AMD - Méthodes numériques pour l'astrophysique

31. Propriétés des estimateurs • Variance de Tn : V[Tn] = E[(Tn-E(Tn))2] • Erreur quadratique de Tn par rapport à θ: EQ(Tn,θ) = E[(Tn-θ)2] • Biais de Tn par rapport à θ: B(Tn,θ) = E[Tn-θ]. • Propriété 1 : si EQ(Tn,θ))0 quand n∞, alors l’estimateur Tnest convergent. • Un estimateur est meilleur qu’un autre si son erreur quadratique est inférieure. • Avec un estimateur biaisé, E[Tn] ≠ θ (décalage). On préfère en général des estimateurs sans biais (B(Tn,θ) = 0). • Un estimateur est asymptotiquement sans biais si B(Tn,θ)0 quand n∞. Master 2 AMD - Méthodes numériques pour l'astrophysique

32. Exemples • On reprend X, loi uniforme entre 0 et θ. On a Tn et Tn’, estimateurs de θ. • E(Tn) = θ B(Tn,θ) = 0; EQ(Tn,θ) = θ2/(3n). • B(Tn’) ~ θ/(2n) ; EQ(Tn’,θ) ~ θ2/naTn est meilleur que Tn’ (non biaisé et plus petite erreur quadratique). • Autre estimateur de θ : Tn″ = max {X1,…,Xn}, convergent et biaisé. • E(Tn″) = nθ/(n+1)  B(Tn″,θ) = -θ/(n+1); EQ(Tn,θ) = 2θ2/(n+1)(n+2) ~ 2θ2/n2. • On construit alors Tn″′=(n+1)Tn″ /n, non biaisé… EQ(Tn,θ) = θ2/n(n+2) ~ θ2/n2.C’est le meilleur des 4 (non biaisé et convergence en n2) Master 2 AMD - Méthodes numériques pour l'astrophysique

33. Convergence des estimateurs • EQ(T,θ) = E[(T-θ)2]= E[(T-E(T)+E(T)-θ)2] = E[(T-E(T))2] + (E(T)-θ)2 + 2(E(T)-θ)×E[T-E(T)] = Var(T) + B(T,θ)2 • Conséquence : Si un estimateur est sans biais (même asymptotiquement), et si sa variance tend vers 0, il est convergent. Master 2 AMD - Méthodes numériques pour l'astrophysique

34. Estimateurs standards • On a (X1,...,Xn), un échantillon de loi inconnue. • Estimateur de l’espérance convergent et non biaiséVar(Xn)=Var(X)/n • Estimateur de la variance Convergent, mais biaisé (variance empirique)E[Sn2]=(n-1)/n×Var(X) • Estimateur non biaisé(variance empirique non biaisée) Master 2 AMD - Méthodes numériques pour l'astrophysique

35. Recherche d’estimateurs • Je cherche à estimer un paramètre physique à partir de statistiques : Comment trouver un bon estimateur ? • Ceci nécessite de faire une hypothèse sur le processus physique, par exemple une loi de probabilité (modélisation de données) • Exemple : méthode des moments. Pour k>0 E[Xk] et E[(X-E(X))k] sont les moments de la variable aléatoire X. • Je suppose une loi de probabilité d’une forme particulière, dépendant de k paramètres a1,...,ak que je cherche à estimer. • Les k premiers moments peuvent s’exprimer en fonction des a1,...,ak Les paramètres s’expriment en fonction des moments. • Il suffit d’avoir des estimateurs des moments (p.ex la moyenne et la variance) pour en déduire des estimateurs des paramètres. • Inconvénient de la méthode : les estimateurs qu’on tire sont souvent peu précis. Master 2 AMD - Méthodes numériques pour l'astrophysique

36. Estimation par ajustement • Idée : définir une métrique sur les lois de probabilités, et trouver la loi de probabilité P qui minimise une « distance » par rapport à la loi empirique Pe. • Si (x1,…,xn) est un échantillon observé, (c1,…,ck) l’ensemble des valeurs prises par les xi’s, on définit Pe(ci)=n(ci)/n, où n(ci) est le nombre de fois où ci est réalisé (fréquence empirique). • Distance du χ2 : • Distance de Kolmogorov-Smirnov = comparaison entre les fonctions de répartition théoriques (F) et empiriques (Fe) Master 2 AMD - Méthodes numériques pour l'astrophysique

37. Méthodes numériques pour l’astrophysique • Techniques de base • Estimateurs et statistique • Modélisation de données • Régression linéaire • Modèles linéaires de moindres carrés • Modèles non-linéaires : méthode de Levenberg-Marquardt • Recuit simulé et algorithmes génétiques • Résolution numérique d’équations différentielles • Equations aux dérivées partielles Master 2 AMD - Méthodes numériques pour l'astrophysique

38. Modélisation de données par régression • On mesure une série de données (x1,y1)… (xn,yn). On cherche à trouver une relation physique y=F(x). • Dans la pratique on se donne un modèle dépendant d’un nombre k de paramètres α1,…,αk: y=f(α1,…,αk; x). • On va chercher les paramètres α1,…,αk qui collent au mieux avec les observations (xi,yi) . • Chaque mesure yi vient avec son erreur standard σi2. Yi=f(α1,…,αk; xi)+Ei , où Ei est une variable aléatoire de variance σi2. • On va chercher les paramètres qui minimisent • Toute la difficulté consiste à • Trouver les paramètres réalisant le minimum • Estimer l’incertitude sur les paramètres. Master 2 AMD - Méthodes numériques pour l'astrophysique

39. Régression linéaire • But : modéliser des données par une relation linéairey = ax+b • Méthode : Minimiser le χ2pour trouver a et b Master 2 AMD - Méthodes numériques pour l'astrophysique

40. Régression linéaire (2) • Erreur sur a et b : • Var(yi)=σi2 Var(aiyi)=ai2 σi2(Somme de variables aléatoires indépendantes) Master 2 AMD - Méthodes numériques pour l'astrophysique

41. Modèles linéaires de moindres carrées • On cherche un modèle qui dépend linéairement des paramètresy = f(α1,…,αk; x) = α1X1(x)+∙∙∙+αkXk(x)où les X1,…,Xk sont des fonctions données à l’avance • On va chercher à minimiseren disant ∂χ2/∂αi=0 pour i=1…k • On introduit la matrice A, n×k, et le vecteur b de longueur n tels que Master 2 AMD - Méthodes numériques pour l'astrophysique

42. Modèles linéaires de moindres carrés • ∂χ2/∂αj=0 pour j=1…k revient à …autrement dit résoudre un système linéaire ! (équations normales) • On appelle fij = [c]-1ij. On a Master 2 AMD - Méthodes numériques pour l'astrophysique

43. Modèles non linéaires • y = f(α1,…,αk; x) où la dépendance est quelconque. • Le jeu de paramètres (α1,…,αk) minimisant le χ2 doit être trouvé de manière itérative. • C’est conceptuellement identique à un problème de minimisation dans un espace de dimension k. • La non-linéarité ne garantit pas l’existence d’un minimum unique problèmes de minima locaux  nécessité de faire de nombreux essais ! • Quand on peut calculer le gradient de la fonction f par rapport aux paramètres, la méthode recommandée est celle de Levenberg-Marquardt. Master 2 AMD - Méthodes numériques pour l'astrophysique

44. Méthode de Levenberg-Marquardt • On appelle α le vecteur (α1,…,αk) courant. • Si le modèle est linéaire, la fonction χ2(α) est une forme quadratique v est un vecteur et H une matrice k×k • Dans ce cas on sait où est le minimum : αm = H-1∙v. De manière équivalente, si on part d’un vecteur α0 , on a • Si la fonction n’est pas trop éloignée d’une forme quadratique, cette formule peut-être une bonne formule d’itération. • Si la fonction est plus compliquée, on se contentera d’une simple descente de gradient • Idée de la méthode : Alterner entre les deux formules • Problème : Il faut la matrice H = matrice Hessienne Master 2 AMD - Méthodes numériques pour l'astrophysique

45. Méthode de Levenberg-Marquardt (2) • Problème : Les termes cjldépendent des dérivées secondes de f qui ne sont pas forcément accessibles • Souvent on laisse tomber les termes en question… Master 2 AMD - Méthodes numériques pour l'astrophysique

46. Méthode de Levenberg-Marquardt (3) • Deux questions : • Comment décider de l’alternance entre les deux formules (1) et (2) ? • Comment fixer la constante dans la formule (2) ? • Levenberg-Marquardt : • cte ~ 1/cjj. En fait cte =1/(λcjj)  (2) devient λcjjδαj = dj • (1) et (2) peuvent se combiner en une seule équation. On définit C’= C + λ×diag(c11,…,ckk) et on remplace (1) et (2) par • Quand λ>>1, (3) revient à (2); quand λ<<1, (3) équivaut à (1) • Recette : • On part d’un choix initial α, et on calcule χ2(α). On prend un petit λ =0.001 • On résout (3) et on calcule l’incrément da. On calcule χ2(α+da). • Si χ2(α+da) > χ2(α), on est loin du minimum  on multiplie λ par 10 et on revient à 2. • Si χ2(α+da) < χ2(α), on approche  On divise λ par 10, on remplace αparα+da, et on revient à 2. Master 2 AMD - Méthodes numériques pour l'astrophysique

47. Intervalles de confiance des paramètres • Une fois qu’on a obtenu les paramètres (α1,…,αk) qui minimisent, on souhaite obtenir des intervalles de dispersion (δα1,…,δαk)des paramètres à un niveau de confiance donné. • La méthode de Levenberg-Marquardt fournit la matrice des covariances F = C-1. Mais celle-ci n’est utilisable que si les erreurs de mesure sont distribuées de manière gaussienne ! Si ce n’est pas le cas, on peut en déduire n’importe quoi !! • Si les erreurs ne sont pas gaussiennes (ou si on ne sait pas), il faut avoir recours à d’autres méthodes statistiques beaucoup plus robustes… Master 2 AMD - Méthodes numériques pour l'astrophysique

48. Avec la matrice des covariances… • Si les erreurs sont gaussiennes, alors on a avec un degré de confiance de 68% (1σ dans la loi normale). • Attention : le jeu d’intervalles de confiances δαj, j=1..n ne constitue pas une région de confiance pour les k paramètres conjointement… Master 2 AMD - Méthodes numériques pour l'astrophysique

49. Méthodes statistiquesou quand les erreurs ne sont pas gaussiennes… • Simulation de jeu de données Monte-Carlo • Idée : tirer plusieurs (N) jeux de données fictives dans les intervalles de confiance ±σi (de manière gaussienne…), et effectuer N minimisations. • On observe la distribution des jeux de paramètres αobtenus, on en déduit des intervalles de confiance. • Méthode Bootstrap • Au lieu de tirer des jeux de données fictives, on tire au hasard n données (avec éventuelle répétition) parmi les n de base, et on refait une minimisation. • On répète l’opération N fois, on obtient la distribution des α. • Avantage de ces méthodes : Elles sont simples, robustes et ne présupposent rien sur la distribution des erreurs de mesure. • Désavantage : Elles sont coûteuses en temps de calcul (il faut refaire N minimisations). Master 2 AMD - Méthodes numériques pour l'astrophysique

50. Méthode MCMC(Monte-Carlo + Chaînes de Markov) • Idée : ne plus chercher d’intervalle de confiance pour les paramètres, mais échantillonner leur distribution statistique. • C’est une amélioration de la méthode de simulation de jeu de données. • On a un vecteur de données y. Pour un jeu de paramètres , on connaît , probabilité d’observer y sachant . On cherche • Par le théorème de Bayes • On suppose connaître (prior). • Une chaîne de Markov = une suite de modèles i, où chaque i+1se déduit de ivia une probabilité de transition . A la fin, la distribution des i échantillonne Master 2 AMD - Méthodes numériques pour l'astrophysique