1 / 111

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

Méthodes numériques pour l’astrophysique M2 « Astrophysique et Milieux Dilués ». Hervé Beust. Laboratoire d’Astrophysique de Grenoble. Méthodes numériques pour l’astrophysique. Techniques de base Estimateurs et statistique Modélisation de données

Download Presentation

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

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  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

More Related