Cours Chapitre 3 Régression Linéaire PDF A3S6 2023-2024

16 avril 2024 Statistiques Année 3 Semestre 6 Chapitre 3 : Régression Linéaire : Simple (RLS) & Multiple (RLM) au sens des Moindres Carrés (MC) Laetitia DELLA MAESTRA Enseignant-chercheur en Mathématiques [email protected] - Bureau L405 Références & Bibliographie Références données aux CMO précédents ; Equipe de Statistique de l’IRMAR (Institut de Recherche Mathématique de Rennes) ; attention : avec le logiciel R, mais présentation des méthodes & explications mathématiques de très grande qualité Site de François Husson : https://husson.github.io Livre Régression avec R, 2ème édition, P.-A. Cornillon, N. Hengartner, E. Matzner-Lober, L. Rouvière, ed. EDP Sciences Machine-Learning avec Scikit-Learn, 3ème édition, Aurélien Géron, ed. Dunod Exploration de données & méthodes stat., L. Bellanger et R. Tomassone, ed. Ellipses Cours de Statistique de la donnée d’A. Jaghdam (ESILV, années précédentes, A2S4 & A3S6). Attention : les seules notations, définitions, appellations, etc... qui font foi dans ce module sont celles du présent cours. 2/97 Plan du cours A) Présentation & Résolution (algébrique/géométrique & analytique) Objectif 1 : approximer, avec les meilleurs paramètres possibles, le nuage de points représentant nos données quantitatives par une forme affine ds R2 par une droite affine, ds R3 par un plan affine, ds Rp+1 (p > 2) par un hyperplan affine Objectif 2 : effectuer des prévisions à l’aide du modèle approximé Objectif 3 : évaluer la qualité du modèle choisi Objectif 4 : quantifier l’erreur commise ds le choix de ce modèle approx B) Modèle probabiliste des erreurs (on impose que l’erreur entre notre modèle théorique et nos données est aléatoire, les paramètres optimaux sont alors eux-mêmes aléatoires et sont des estimateurs des paramètres théoriques) : Estimation des paramètres, calcul de leurs Espérance & Variance ↝ contrôle local d’ordre 1 et 2 de l’erreur commise avec le modèle choisi C) Modèle gaussien des erreurs (on impose que l’erreur entre notre modèle théorique et nos données est gaussienne, les paramètres optimaux sont alors eux-mêmes gaussiens) : Loi des paramètres, IC & Tests ↝ contrôle global de l’erreur commise avec le modèle choisi D) Limites du modèle : Retour critique sur les conditions imposées Objectif 5 : Sélection de variables ↝ parcimonie : peut-on faire quasiment aussi bien avec moins de facteurs explicatifs, c-à-d en dim + petite ? 3/97 A) Présentation Problématique de la Régression On dispose de données quantitatives, discrètes ou continues, concernant n individus I1 ,... , In et p + 1 caractéristiques appelées variables x1 ,... , xp , y y joue le rôle de la variable à expliquer x1 ,... , xp jouent le rôle de variables explicatives autrement dit on dispose de n séries de p + 1 valeurs : (x1. = (x11 ,... , x1p ) , y1 ),... , (xn. = (xn1 ,... , xnp ) , yn ) et l’on cherche à expliquer et à prédire y par x1 ,... , xp ⇒ on cherche quel type de relation il existe entre les xi. et les yi c-à-d qu’on cherche une fonction f telle que pour tout i ∈ J1, nK, f (xi. ) ≃ yi prev ⇒ on cherche à prévoir le + précisément possible l’étiquette yn+1 associée prev à un nouveau vecteur non étiqueté x(n+1). en posant yn+1 ∶= f (x(n+1). ). 4/97 Problématique de la Régression Linéaire On cherche une relation de type linéaire entre xi. et yi , càd que l’on se restreint aux fcts affines f ∶ x. = (x1 ,... , xp ) ∈ Rp ↦ β0 + β1 x1 +... + βp xp où β0 , β1 ,... , βp ∈ R autrement dit, l’on suppose que les points (xi. = (xi1 ,... , xip ) , yi ) de l’espace affine Rp+1 sont regroupés autour d’un hyperplan affine théorique d’équation y = β0 + β1 x1 +... + βp xp Remarque : yi peut être vue comme l’étiquette associée au i-ème vecteur xi. : la Régression est une méthode d’Apprentissage Supervisé 5/97 Exemple pour p = 1 : Un agent immobilier reçoit un nouveau client lui demandant de prendre en charge la vente de son appartement. Pour fixer au mieux le prix de vente de ce bien en fonction de sa surface, l’agent étudie d’abord l’historique de ses ventes dans le quartier où est situé cet appartement : il s’agit de n = 20 appartements pour lesquels il dispose du prix ↝ variable y de la surface ↝ variable x1 , que l’on notera x (il n’y a pas d’autre variable explicative !) Objectif : expliquer y à partir de x Numéro i du bien Surface xi (en m2 ) Prix de vente yi (en K euros) Numéro i du bien Surface xi (en m2 ) Prix de vente yi (en K euros) 6/97 1 28 130 11 32 155 2 50 280 12 52 245 3 55 268 13 40 200 4 110 500 14 70 325 5 60 320 15 28 85 6 48 250 16 30 78 7 90 378 17 105 375 8 35 250 18 52 200 9 86 350 19 80 270 10 65 300 20 20 85 Autres exemples pour p = 1 : Exemple 2 pour p = 1 : Une étude a été menée auprès de 12 étudiants afin d’expliquer le score à un examen d’anglais à partir du temps consacré à la préparation de cet examen. Pour chaque étudiant, on dispose du temps de révision en heures (variable x), du score obtenu sur 800 points (variable y) étudiant n○ i xi yi 1 4 390 2 9 580 3 10 650 4 14 730 5 4 410 6 7 530 7 12 600 8 1 350 9 3 400 10 8 590 11 11 640 12 5 450 Ainsi avec une préparation de 4 heures, l’étudiant n○ 1 a obtenu le score de 390 à l’examen, avec une préparation de 9 heures, l’étudiant n○ 2 a obtenu le score de 580 à l’examen. Exemple 3 pour p = 1 : On étudie l’évolution du nombre d’inscriptions à un jeu en ligne au cours du temps. Pour chaque mois de l’année 2022, on dispose du rang du mois (variable x ; janvier est de rang 1, février est de rang 2,...) et du nombre d’inscriptions en milliers (variable y). Les résultats sont xi = mois n○ i yi 1 37 2 43 3 41 4 40 5 51 6 47 7 48 8 54 9 56 10 64 11 66 12 73 Ainsi, au mois de janvier 2022, il y a eu 37000 inscriptions au jeu, en février 2022 il y a eu 43000 inscriptions au jeu. Exemple 4 pour p = 1 : Avant la commercialisation d’un produit, une entreprise effectue une étude de marché afin de déterminer la quantité demandée en milliers (variable y) en fonction du prix de vente en euros (variable x). Pour 6 prix de vente différents, les résultats sont : n○ i xi yi 7/97 1 15 44.4 2 20 27.0 3 25 16.3 4 30 10.0 5 35 6.2 6 40 3.5 Exemple pour p = 2 : On souhaite expliquer le poids d’un basketteur professionnel de la NBA à partir de sa taille et de son âge. Ainsi, pour n = 505 basketteurs de la NBA (données d’il y a dix ans) , on dispose : de leur poids (en livres) de leur taille (en inches) ↝ variable y ↝ variable x1 de leur âge ↝ variable x2 Objectif : expliquer y à partir de x1 et x2 8/97 Exemple pour p = 5 : Nous disposons d’un tableau de notes obtenues par n = 8 élèves A(riane), B(asile), C(alliope), D(enis), E(ugénie), F(ilippo), G(eorgie), H(ector) dans p = 5 matières Statistique, Informatique, Marketing, Finance, Comptabilité ↝ les variables x1 ,... , x5 et de leur moyenne pondérée arrondie calculée à partir de ces notes ↝ variable y Si l’on ne sait pas quels coefficients ont été appliqués pour chaque matière, comment les retrouver (approximativement) ? 9/97 Formalisation du problème de Régression Linéaire au sens des Moindres Carrés (MC) Cette méthode est dite Régression Linéaire au sens des Moindres Carrés ( = least squares method) lorsqu’en plus d’imposer que la fonction f soit affine : f ∶ x. = (x1 ,... , xp ) ∈ Rp ↦ β0 + β1 x1 +... + βp xp où β0 , β1 ,... , βp ∈ R on impose le critère n 2 (∀i ∈ J1, nK, f (xi. ) ≃ yi ) ⇔ ( ∑ (yi − f (xi. )) i=1 le + petit possible ) c-à-d que l’on impose de choisir (β0 , β1 ,... , βp ) ∈ Rp+1 tel que n 2 ∑ ( yi − (β0 + β1 xi1 +... + βp xip ) ) soit le plus petit possible i=1 10/97 Lorsque p > 1, on parle de Régression Linéaire Multiple au sens des Moindres Carrés p = 1, on parle de Régression Linéaire Simple au sens des Moindres Carrés, et cette méthode revient à imposer, en plus du fait que la fonction f soit affine f ∶ x ∈ R ↦ β0 + β1 x où β0 , β1 ∈ R le critère n 2 (∀i ∈ J1, nK, f (xi ) ≃ yi ) ⇔ ( ∑ (yi − f (xi )) le + petit possible ) i=1 c-à-d que l’on impose de choisir (β0 , β1 ) ∈ R2 tel que n 2 ∑ ( yi − (β0 + β1 xi )) ) soit le plus petit possible i=1 A partir de maintenant nous dirons juste RLS (resp. RLM) à la place de Régression Linéaire Simple (resp. Multiple) au sens des Moindres Carrés 11/97 Réécriture du problème d’optimisation On cherche donc le (p + 1)-uplet (β̂0 , β̂1 ,... , β̂p ) qui minimise la somme des erreurs quadratiques ϵ2i , où ϵi ∶= yi − (β0 + β1 xi1 +... + βp xip ) : (β̂0 , β̂1 ,... , β̂p ) = = n Argmin ∑ ( yi − (β0 + β1 xi1 +... + βp xip ) ) (β0 ,β1 ,...,βp )∈Rp+1 i=1 n Argmin 2 2 ∑ ϵi (β0 ,β1 ,...,βp )∈Rp+1 i=1 Pr i ∈ J1, nK, ŷi ∶= β̂0 + β̂1 xi1 +... + β̂p xip ↝ valeur prédite pour yi Pr i ∈ J1, nK, ϵ̂i ∶= yi − ŷi ↝ résidu dans la prédiction de yi par ŷi n min (β0 ,β1 ,...,βp )∈R n = ∑ ( yi i=1 i=1 2 n − (β̂0 + β̂1 xi1 +... + β̂p xip ) ) = ∑ ϵ̂2i ↝ somme des carrés des résidus (scr) 12/97 2 ∑ ( yi − (β0 + β1 xi1 +... + βp xip ) ) p+1 i=1 Réécriture du problème d’optimisation : cas p = 1 On cherche ici le couple (β̂0 , β̂1 ) qui minimise la somme des erreurs quadratiques ϵ2i , où ϵi ∶= yi − (β0 + β1 xi ) : n (β̂0 , β̂1 ) = Argmin ∑ ( yi − (β0 + β1 xi ) ) 2 (β0 ,β1 )∈R2 i=1 n = Argmin ∑ ϵ2i (β0 ,β1 )∈R2 i=1 Pr i ∈ J1, nK, ŷi ∶= β̂0 + β̂1 xi ↝ valeur prédite pour yi Pr i ∈ J1, nK, ϵ̂i ∶= yi − ŷi ↝ résidu dans la prédiction de yi par ŷi n min ∑ ( yi − (β0 + β1 xi ) ) 2 (β0 ,β1 )∈R i=1 n 2 2 n = ∑ ( yi − (β̂0 + β̂1 xi ) ) = ∑ ϵ̂2i i=1 i=1 ↝ somme des carrés des résidus (scr) 13/97 Interprétation graphique de la RL (au sens des MC) cas p = 1 (RLS) : β̂0 , β̂1 sont construits pour minimiser les "distances verticales", au sens portées par des droites parallèles à (Oy ), entre les observations (xi , yi ) et la droite de régression théorique y = β0 + β1 x. Exemple pour p = 1 : pour les données (50, 280), (55, 268), (60, 320) extraites de nos données appartements = (surfaces,prix), et deux droites quelconques du type y = b0 + b1 x 14/97 cas p = 2 : β̂0 , β̂1 , β̂2 sont construits pour minimiser les "distances verticales", au sens portées par des droites parallèles à (Oz), entre les observations (xi1 , xi2 , yi ) et le plan affine de régression théorique y = β0 + β1 x1 + β2 x2 Exemple pour nos données nba avec un plan quelconque : cas p ≥ 3 : représentation graphique impossible ! 15/97 Récapitulatif des objectifs On veut donc : estimer de manière optimale β0 , β1 ,... , βp à l’aide des données (x1. = (x11 ,... , x1p ) , y1 ),... , (xn. = (xn1 ,... , xnp ) , yn ) mesurer l’importance des variables x1 ,... , xp dans l’explication de y prev la "valeur moyenne" de y pour une prédire avec précision ŷn+1 nouvelle valeur (xn+1,1 ,... , xn+1,p ) de x1 ,... , xp , Lorsque p = 1 cela devient estimer de manière optimale β0 , β1 à l’aide des données (x1 , y1 ),... , (xn , yn ) mesurer l’importance de la variable x dans l’explication de y prev prédire avec précision ŷn+1 , la "valeur moyenne" de y pour une nouvelle valeur xn+1 de x 16/97 A) Résolution du problème de RL au sens des MC Conditions de dimension & de rang Ns ns placerons ds tte la suite sous la double condition (A0) suivante : n > p + 1 ↝ dans nos données il y a plus d’individus que de variables ⎛1⎞ 1n,1 = ⎜ ⋮ ⎟ , x.1 ,... , x.p forment une famille libre de p + 1 vecteurs de ⎝1⎠ n R (c-à-d ces p + 1 vecteurs sont linéairement indépendants) Lorsque p = 1, la condition est que les vecteurs 1n,1 et x. soient linéairement indépendants, c-à-d ∃ 1 ≤ i ≠ k ≤ n tq xi ≠ xk 1 Lorsque p > 1, la condition équivaut à x.. ∶= ( ⋮ 1 x11 ⋮ xn1......... x1p ⋮ xnp ) de rang p + 1 ou encore à t x.. x.. inversible (on dit alors que x.. est de rang plein puisque x.. ∈ Mn,p+1 (R) et n > p + 1) (puisque t x.. x.. ∈ Mp+1 (R), Ker(x.. ) = Ker( t x.. x.. ) et, d’après le Théorème du rang, p + 1 = rg(x.. ) + dim(Ker(x.. )) = rg( t x.. x.. ) + dim(Ker( t x.. x.. )) 17/97 Résolution du problème d’optimisation : cas p = 1 On veut résoudre le problème d’optimisation de la RLS au sens des MC : (β̂0 , β̂1 ) = Argmin C (β0 , β1 ) (β0 ,β1 )∈R2 n 2 où C ∶ (β0 , β1 ) ∈ R × R ↦ ∑ (yi − (β0 + β1 xi )) i=1 ↝ C peut être vue comme une fonction modélisant le coût d’utilisation de la droite de régression théorique y = β0 + β1 x C est une fonction deux fois dérivable sur R × R Recherche des points critiques de C c-à-d des points (β0 , β1 ) ∈ R2 tel 0 que ∇C (β0 , β1 ) = ( ) où ∇C (β0 , β1 ) est le gradient de C en (β0 , β1 ) : 0 n ∇C (β0 , β1 ) = ∂C (β0 , β1 )⎞ ⎛ ∂β 0 ∂C ⎝ ∂β1 (β0 , β1 )⎠ ⎛ −2 ∑ (yi − (β0 + β1 xi )) ⎞ i=1 ⎟ =⎜ ⎜ n ⎟ ⎝−2 ∑ xi (yi − (β0 + β1 xi ))⎠ i=1 18/97 Cela revient à résoudre le syst. de 2 équations à 2 inconnues (β0 , β1 ) : n ⎧ ⎪ (yi − (β0 + β1 xi )) = 0 −2 ⎪ ∑ ⎪ ⎪ i=1 ⎨ n ⎪ ⎪ −2 ∑ xi (yi − (β0 + β1 xi )) = 0 ⎪ ⎪ ⎩ i=1 n n ⎧ ⎪ y = nβ + β xi ⎪ ∑ ∑ i 0 1 ⎪ ⎪ i=1 ⇔ ⎨ i=1 n n n ⎪ ⎪ x y = β x + β xi2 ⎪ ∑ ∑ ∑ i i 0 i 1 ⎪ ⎩ i=1 i=1 i=1 On trouve un unique point critique, que l’on note (β̂0 , β̂1 ) : ⎧ β̂0 = y − β̂1 x ⎪ ⎪ ⎪ n ⎪ 1 ⎪ (x −x )(yi −y ) n∑ i ⎨ i=1 ⎪ β̂ = n 1 ⎪ ⎪ 1 ⎪ (x −x )2 ⎪ n∑ i ⎩ i=1 Rq : 1 n n ∑ (xi − x )2 = σ̂x2 ≥ 0 et est σ̂x2 = 0 ssi x1 = x2 =... = xn , ce qui i=1 est impossible puisqu’on a supposé qu’il existe 1 ≤ i ≠ k ≤ n tq xi ≠ xk , donc σ̂x2 > 0 19/97 On détermine la Hessienne de C en (β̂0 , β̂1 ) : 2 ⎛ ∂ C (β̂0 , β̂1 ) ⎜ ∂β02 HessC (β̂0 , β̂1 ) = ⎜ ⎜ ∂2 C (β̂ , β̂ ) ⎝ ∂β0 ∂β1 0 1 ↝ elle est définie positive ∂ 2 C (β̂ , β̂ )⎞ 0 1 ∂β0 ∂β1 ⎟ ⎟ ∂ 2 C (β̂ , β̂ ) ⎟ 0 1 ⎠ ∂β 2 1 n ⎛ 2n ⎜ =⎜ n ⎜ ⎝2 ∑ xi i=1 local n ⎛ 2n ⎜ HessC (β0 , β1 ) = ⎜ n ⎜ ⎝2 ∑ xi i=1 i=1 (son déterminant est égal à 4n2 σ̂x2 qui est strictement positif car les xi ne sont pas tous égaux, et sa trace est égale à 2n(1 + σ̂x2 + (x )2 ) qui est strictement positive) bien un minimum 2 ∑ xi ⎞ i=1 ⎟ ⎟ n ⎟ 2 ∑ xi2 ⎠ en (β̂0 , β̂1 ) donc C atteint (en fait, pour tout (β0 , β1 ) ∈ R2 , 2 ∑ xi ⎞ i=1 ⎟ ⎟ donc C est convexe sur R2 , et C atteint bien son minimum global en (β̂0 , β̂1 ) ) n ⎟ 2 ∑ xi2 ⎠ i=1 Conclusion pour la RLS : La droite de régression linéaire au sens des ⎧ β̂0 = y − β̂1 x ⎪ ⎪ ⎪ n ⎪ 1 ⎪ (x −x )(yi −y ) n∑ i moindres carrés est y = β̂0 + β̂1 x avec ⎨ ĉ σ̂ i=1 ⎪ β̂ = = σ̂x ,y2 = σ̂yx ρ̂x ,y n 1 ⎪ ⎪ 1 x 2 ⎪ (x −x ) ∑ i ⎪ n ⎩ i=1 où ĉx ,y ∶= 1 n n ∑ (xi − x )(yi − y ) est la covariance de x et y et ρ̂x ,y ∶= i=1 ĉx ,y σ̂x σ̂y le coefficient de corrélation linéaire de x et y (aussi dit de Pearson) 20/97 est ĉx ,y ↝ x.cov(y, ddof=0) ou np.cov(x, y, ddof=0)[0,1] ρ̂x ,y ↝ x.corr(y) ou np.corrcoef(x, y)[0,1] Comme y = β̂0 + β̂1 x , cette droite passe par le point de coordonnées (x , y ) càd le centre de gravité du nuage de points (aussi appelé point moyen, centre d’inertie, isobarycentre du nuage de points) Comme σ̂x , σ̂y > 0, le coefficient directeur β̂1 de cette droite et ρ̂x ,y sont de même signe : la fonction associée est croissante (resp. décroissante) ssi ρ̂x ,y est positif (resp. négatif) Cela permet de deviner le signe de ρ̂x ,y avec la silhouette du nuage de points (xi , yi )i∈J1,nK ! ρ̂x ,y mesure à quel point x et y sont liées linéairement ; ρ̂x ,y ∈ [−1, 1] (conséquence de l’inégalité de Cauchy-Schwarz) et, plus ∣ρ̂x ,y ∣ est proche de 1, plus la liaison linéaire entre x et y est forte. Rq : ĉx ,x = σ̂x2 et ρ̂x ,x = 1 Pour α, α′ , β, β ′ , γ, γ ′ ∈ R, u, v ∈ Rn , ĉαx +βu+γ,α′ y +β ′ v +γ ′ = αα′ ĉx ,y + αβ ′ ĉx ,v + βα′ ĉu,y + ββ ′ ĉu,v 21/97 Les graphiques ci-dessous illustrent le lien existant entre la pertinence de l’ajustement d’un nuage de points par une droite, caractérisée par la corrélation linéaire entre x et y, et son coefficient ρ̂x ,y : Source : wikipedia 22/97 Résolution du problème d’optimisation : cas p > 1 ⎛β0 ⎞ ⎛1 x11... x1p ⎞ ⎛y1 ⎞ ⎜β ⎟ ⋮... ⋮ ⎟, β. ∶= ⎜ 1 ⎟, nous avons En notant : y. = ⎜ ⋮ ⎟, x.. ∶= ⎜ ⋮ ⎜ ⋮ ⎟ ⎝1 xn1... xnp ⎠ ⎝yn ⎠ ⎝βp ⎠ n 2 C (β. ) = C (β0 , β1 ,... , βp ) ∶= ∑ ( yi − (β0 +β1 xi1 +...+βp xip ) ) = ∣∣y. −x.. β. ∣∣2Rn i=1 β̂. = argmin C (β0 , β1 ,... , βp ) = argmin ∣∣y. − x.. β. ∣∣2Rn On cherche donc : β. ∈ Rp+1 β. ∈ Rp+1 ce qui équivaut à chercher argmin ∣∣y. − z. ∣∣2Rn z. ∈ Im(x.. ) Si β̂. = argmin ∣∣y. − x.. β. ∣∣2 n R β. ∈ Rp+1 , on a : ∀β. ∈ Rp+1 , ∣∣y. − x.. β̂. ∣∣2 n ≤ ∣∣y. − x.. β. ∣∣2 n ce qui est équivaut (puisque Im(x.. ) R R est par déf. l’ensemble des x.. β. avec β. variant dans Rp+1 ), en notant ŷ. = x.. β̂. (qui par déf. ∈ Im(x.. )) à ∀z. ∈ Im(x.. ), ∣∣y. − ŷ. ∣∣2 n ≤ ∣∣y. − z. ∣∣2 n , ce qui équivaut aussi à ŷ. = argmin ∣∣y. − z. ∣∣2 n R R Inversement, si z.⋆ ∈ Im(x.. ) est tel que z.⋆ = z. ∈ Im(x.. ) R argmin ∣∣y. − z. ∣∣2 n , on a ∀z. ∈ Im(x.. ), ∣∣y. − z.⋆ ∣∣2 n ≤ ∣∣y. − z. ∣∣2 n , ce qui R R R z. ∈ Im(x.. ) équivaut à ∀β. ∈ Rp+1 , ∣∣y. − z.⋆ ∣∣2 n ≤ ∣∣y. − x.. β. ∣∣2 n Comme z.⋆ ∈ Im(x.. ), ∃β.⋆ ∈ Rp+1 tq z.⋆ = x.. β.⋆ , et l’on a R R ∀β. ∈ Rp+1 , ∣∣y. − x.. β.⋆ ∣∣2 n ≤ ∣∣y. − x.. β. ∣∣2 n donc β.⋆ = argmin ∣∣y. − x.. β. ∣∣2 n R R R β. ∈ Rp+1 23/97 Propriété fondamentale de la RLM (A CONNAITRE PAR COEUR) Sous la condition double (A0) : n > p + 1 et t x.. x.. inversible le problème de RLM argmin ∣∣y. − x.. β. ∣∣2Rn admet une unique solution : β. ∈ Rp+1 β̂. = ( t x.. x.. )−1 t x.. y ↝ ∀j ∈ {0, 1,. , p}, l’estimateur des moindres carrés de βj est β̂j ↝ ŷ. = x.. β̂. = β̂0 1n,1 + β̂1 x.1 +... + β̂p x.p = x.. (t x.. x.. )−1 t x.. y. = ΠIm(x.. ) (y. ) c-à-d que le vecteur des valeurs prédites, ŷ. , est le projeté orthogonal du vecteur y. des valeurs de la variable à expliquer sur l’espace engendré par les colonnes de x.. , c-à-d par 1n,1 (l’intercept) et x.1 ,... , x.p les vecteurs contenant les valeurs des variables explicatives ↝ l’hyperplan affine de régression linéaire au sens des moindres carrés est y = β̂0 + β̂1 x1 +... + β̂p xp ⊥ (y ) où Π ⊥ est ↝ ϵ̂. = y. − ŷ. = (In − x.. (t x.. x.. )−1 t x.. )y. = Π (Im(x.. )). (Im(x.. )) ⊥ la matrice de projection orthogonale sur (Im(x.. )) 24/97 −1 Remarque 1 : Cette matrice x.. ( t x.. x.. ) t x.. de projection orthogonale sur Im(x.. ) est dite hat matrix (ou encore projection matrix ou influence matrix ) car elle met des −1 "hat" sur y. puique ŷ. = x.. ( t x.. x.. ) t x.. y.. ↝ elle envoie le vecteur "valeurs réponses" (= response values ) sur le vecteur "valeurs prédites" (fitted values). En tant que matrice de projection orthogonale, elle est −1 t symétrique ↝ t (x.. ( t x.. x.. ) −1 t idempotente ↝ (x.. ( t x.. x.. ) 25/97 −1 t x.. ) = x.. ( t x.. x.. ) 2 x.. −1 t x.. ) = (x.. ( t x.. x.. ) x.. ) ⎛ 1... 1 ⎞ 1 x11... xn1 ⎟ ⎛ ⎜x ⎟⎜⋮ ⋮ Remarque 2 : on a t x.. x.. = ⎜ 11 ⎜ ⋮... ⋮ ⎟⎝ ⎝x1p... xnp ⎠ 1 xn1 n x.1 nx.2... ⎛ n 2 ⎜n x.1 ∣∣x ∣∣ ⟨x , x ⟩.1 Rn.1.2 Rn... ⎜ t ⎜ ∣∣x.2 ∣∣2Rn... donc x.. x.. = ⎜n x.2 ⟨x.2 , x.1 ⟩Rn ⎜ ⋮ ⋮ ⋮ ⋱ ⎜ ⎝n x.p ⟨x.p , x.1 ⟩Rn ⟨x.p , x.2 ⟩Rn... n n x. Cela donne lorsque p = 1 : t x.. x.. = ( ) n x. ∣∣x. ∣∣2Rn 2 n n i=1 2 i=1... x1p ⎞... ⋮ ⎟... xnp ⎠ nx.p ⎞ ⟨x.1 , x.p ⟩Rn ⎟ ⎟ ⟨x.2 , x.p ⟩Rn ⎟ ⎟ ⎟ ⋮ ⎟ 2 ∣∣x.p ∣∣Rn ⎠ d’où 2 det( t x.. x.. ) = n∣∣x. ∣∣2Rn − n2 ( x. ) = n2 ( n1 ∑ xi2 − ( n1 ∑ xi ) ) = n2 σ̂x2. (qui est > 0 si les ∣∣x. ∣∣Rn −n x. ) ce qui permet, avec −n x. n ⎧ β̂0 = y. − β̂1 x. ⎪ ⎪ ⎪ n ⎪ 1 ⎪ (x −x. )(yi −y. ) n∑ i β̂. = ( t x.. x.. )−1 t x.. y. , de retrouver ⎨ i=1 ⎪ β̂ = n 1 ⎪ ⎪ 1 ⎪ (x −x. )2 ⎪ n∑ i ⎩ i=1 xi ne sont pas ts égaux) 26/97 −1 ( t x.. x.. ) = 1 n2 σ̂x2. ( Preuve n○ 1 : Projection orthogonale (1/4) argmin ∣∣y. − z. ∣∣2Rn est atteint en z. = ΠIm(x.. ) (y. ) = ΠIm(x.. ) y. z. ∈ Im(x.. ) où ΠIm(x.. ) ∈ Mn (R) est la matrice de projection orthogonale sur Im(x.. ) (et également, par abus de notation, l’application linéaire qui lui est canoniquement associée) ↝ donc ŷ. = x.. β̂ = ΠIm(x.. ) (y. ) (Rq : on a bien ŷ. ∈ Rn , puisque, d’une part, x.. ∈ Mn,p+1 (R) et β̂ ∈ Rp+1 , et d’autre part ΠIm(x.. ) ∈ Mn (R) et y. ∈ Rn ) En effet, d’après le Théorème de Pythagore, ∀z. ∈ Im(x.. ), ∣∣y. − z. ∣∣2Rn = ∣∣y. − ΠIm(x.. ) (y. ) + ΠIm(x.. ) (y. ) − z. ∣∣2Rn = ∣∣y. − ΠIm(x.. ) (y. )∣∣2Rn + ∣∣ΠIm(x.. ) (y. ) − z. ∣∣2Rn car y. − ΠIm(x.. ) (y. ) = Π ⊥ (y. ) ∈ (Im(x.. )) (Im(x.. )) ΠIm(x.. ) (y. ) − z. ∈ Im(x.. ) puisque Im(x.. ) est un e.v., ΠIm(x.. ) (y. ) ⊥ ∈ Im(x.. ) et z. ∈ Im(x.. ) Donc le minimum est atteint en un unique point : lorsque ∣∣ΠIm(x.. ) (y. ) − z. ∣∣2Rn = 0 , c-à-d lorsque z. = ΠIm(x.. ) (y. ) 27/97 Preuve n○ 1 : Projection orthogonale (2/4) Notons z.⋆ = ΠIm(x.. ) (y. ). On a alors z.⋆ ∈ Im(x.. ), il existe donc β.⋆ ∈ Rp+1 tel que z.⋆ = x.. β.⋆ , et, comme x.. est de rang plein, ce β.⋆ est unique, (en effet, s’il existait β̃. ∈ Rp+1 tel que x.. β.⋆ = x.. β̃. , alors on aurait x.. ( β.⋆ − β̃. ) = 0n,1 ; or comme x.. est supposé de rang plein (c’est-à-dire de rang p + 1 puisque x.. ∈ Mn,p+1 (R) et n > p + 1), on a d’après le Théorème du rang que dim(Ker(x.. )) = 0, d’où Ker(x.. ) = (0p+1,1 ), et donc β.⋆ = β̃. ) Avec les notations ŷ. = est l’unique solution de et l’on a β.⋆ = argmin ∣∣y. − x.. β. ∣∣2Rn β. ∈ Rp+1 = ΠIm(x.. ) (y. ) et β̂. = β.⋆ argmin ∣∣y. − x.. β. ∣∣2Rn β. ∈ Rp+1 z.⋆ on a que ŷ. = x.. β̂. et β̂. Il reste à déterminer l’expression explicite : ↝ de la solution β̂. du pb et de ŷ. ↝ de la matrice ΠIm(x.. ) de projection orthogonale sur l’espace engendré par les colonnes de la matrice x.. 28/97 Preuve n○ 1 : Projection orthogonale (3/4) Méthode 1 : ∀v ∈ Im(x.. ) , ⟨v , (In − ΠIm(x.. ) ) y. ⟩ = 0 ⇔ ∀u ∈ Rp+1 , ⟨x.. u, (In − ΠIm(x.. ) ) y. ⟩ = 0 ⇔ ∀u ∈ Rp+1 , t u. t x.. (In − ΠIm(x.. ) ) y. = 0 ⇔ t x.. (In − ΠIm(x.. ) ) y. = 0p+1,1 ⇔ t x.. y. = t x.. ΠIm(x.. ) y. ⇔ t x.. y. = t x.. x.. β̂. t ⇔ ( x.. x.. ) −1 t x.. y. = β̂. car t x.. x.. inversible Méthode 2 : ⊥ (In − ΠIm(x.. ) ) y. ∈ (Im(x.. )) ⇔ ∀j ∈ J0, pK, x.j ⊥ (In − ΠIm(x.. ) ) y. ⇔ ∀j ∈ J0, pK, x.j ⊥ (y. − x.. β̂. ) t −1 t ⇔ t x.. y. = t x.. x.. β̂. ⇔ β̂. = ( x.. x.. ) 29/97 x.. y. car t x.. x.. inversible Preuve n○ 1 : Projection orthogonale (4/4) Conclusion ↝ Détermination de l’expression explicite de ΠIm(x.. ) : t Comme β̂. = ( x.. x.. ) −1 t x.. y. et ŷ = x.. β̂. = ΠIm(x.. ) y. , on a t −1 t x.. ( x.. x.. ) x.. y. = ΠIm(x.. ) y. Et comme ceci est vrai pour tout y. ∈ Rn , on a finalement t −1 t ΠIm(x.. ) = x.. ( x.. x.. ) x.. Rq : x.. est de taille n × (p + 1), donc t x.. est de tailles (p + 1) × n, et t x.. x.. est de taille (p + 1) × (p + 1), de même que (t x.. x.. )−1 d’où finalement (t x.. x.. )−1 t x.. est de taille (p + 1) × n et x.. (t x.. x.. )−1 t x.. est de taille n × n ; enfin comme y. est de taille n × 1, β̂. est de taille (p + 1) × 1 30/97 Preuve n○ 2 : Optimisation On cherche β̂. = argmin C (β0 , β1 ,. , βp ) où C (β0 , β1 ,. , βp ) = ∣∣y. − x.. β. ∣∣2Rn β. ∈ Rp+1 Condition d’ordre 1 c-à-d recherche des points critiques de C càd des points β. = (β0 , β1 ,... , βp ) ∈ Rp+1 annulant le gradient de C ∂C ⎛ ∂β0 (β0 , β1 ,... , βp )⎞ 0 ⎜ ∂C ⎟ ⎛ ⎞ ⎜ ∂β1 (β0 , β1 ,... , βp )⎟ ⎜0⎟ ⎟ = ⎜ ⎟ = 0p+1,1 c-à-d tel que ∇C (β. ) = ∇C (β0 , β1 ,... , βp ) = ⎜ ⎟ ⎜⋮⎟ ⎜ ⋮ ⎜ ⎟ ⎝ ⎠ 0 ⎝ ∂C (β , β ,... , β )⎠ ∂βp 0 1 p Méthode 1 ↝ différentiabilité : comme C différentiable sur Rp+1 et ∀h. ∈ Rp+1 , ∣∣y. − x.. (β. + h. )∣∣2Rn = ∣∣y. − x.. β. ∣∣2Rn − 2⟨y. − x.. β. , x.. h. ⟩Rn + ∣∣x.. h. ∣∣2Rn h. ∈ Rp+1 ↦ −⟨2(y. − x.. β. ), x.. h. ⟩Rn = −⟨2 t x.. (y. − x.. β. ), h. ⟩Rp+1 linéaire ∣∣x.. h. ∣∣2Rn = o( ∣∣h. ∣∣Rp+1 ) on a par identification que ∇C (β. ) = −2 t x.. (y. − x.. β. ) , et donc les points critiques β. de C sont caractérisés par t x.. y. = t x.. x.. β. , c-à-d, comme t x.. x.. t −1 inversible, par β. = ( x.. x.. ) t x.. y. 31/97 Preuve n○ 2 : Optimisation Méthode 2 ↝ calcul composante par composante de ∇C (β. ) ∣∣y. − x.. β. ∣∣2Rn = t (y. − x.. β. ) (y. − x.. β. ) = ( t y. − t β. t x.. ) (y. − x.. β. ) = t y. y. − t y. x.. β. − t β. t x.. y. + t β t x.. x.. β. = t y. y. − 2 t β. t x.. y. + t β. t x.. x.. β. car t (x.. β. ) = t β. t x.. et, comme t y. x.. β. est un réel, il est égal à sa transposée c-à-d t y. x.. β. = t β. t x.. y. Notons (e.0 , e.1 ,... , e.p ) la base canonique de Rp+1 (c-à-d ∀j ∈ {0, 1,... , p}, e.j est le vecteur colonne à p + 1 composantes avec p composantes nulles sauf la j + 1-ème qui vaut 1 ; rq : on a donc pr tt u = t (u0 , u1 ,... , up ) ∈ Rp+1 , t e.j u = uj est la (j + 1)-ème composante de u) Soit j ∈ {0, 1,. , p}, déterminons ∂C ∂βj (β. ) : par linéarité de la dérivation, on a ∂C ∂(β. ↦ t y. y. ) ∂(β. ↦ −2 t β. t x.. y. ) ∂(β. ↦ t β. t x.. x.. β (β. ) = (β. ) + (β. ) + ∂βj ∂βj ∂βj ∂βj = 0 − 2 t e.j t x.. y. + t e.j t x.. x.. β. +t β. t x.. x.. e.j = −2 t e.j t x.. y. + 2 t e.j t x.. x.. β car t e.j t x.. x.. β. est un réel, donc il est égal à sa transposée, d’où t e.j t x.. x.. β. = t β.t x.. x.. e.j ∂C Donc ∇C (β. ) = 0 ⇔ ∀j ∈ {0, 1,... , p} , ∂β (β. ) = 0 ⇔ −2 t e.j t x.. y. + j 2 t e.j t x.. x.. β. = 0 ⇔ t e.j t x.. x.. β. = t e.j t x.. y. ⇔ t x.. x.. β. = t x.. y. 32/97 Condition d’ordre deux : montrons que la fct C est convexe c-à-d que sa Hessienne est définie positive, ce qui impliquera que le point t −1 critique (unique) β̂. = ( x.. x.. ) t x.. y. est le minimum (unique) de C Notons HC (β. ) la matrice Hessienne de C au point β. : ∀β. ∈ Rp+1 , HC (β. ) = 2 t x.. x.. (et donc HC est constante) : ∀j, k ∈ {0, 1,. , p} ∂C ) ∂( ∂β ∂(β. ↦ −2 t e.j t x.. y. + 2 t e.j t x.. x.. β. ) ∂2C j (β. ) = (β. ) = (β. ) ∂βj ∂βk ∂βk ∂βk = −2 ∂(β. ↦ t e.j t x.. y. ) ∂(β. ↦ 2 t e.j t x.. x.. β. ) +2 ∂βk ∂βk t t t t = 0 + 2 e.j x.. x.. e.k = 2 e.j x.. x.. e.k HC est définie positive : pour tout z. = (z0 , z1 ,... , zp ) ∈ Rp+1 non nul, t z HC (β. ) z. = t z. (2 t x.. x.. ) z. = 2∣∣x.. z. ∣∣2p+1 > 0. En effet, R 2 ∗ ∣∣x.. z. ∣∣ p+1 ≥ 0 en tant que norme au carré R ∗ ∣∣x.. z. ∣∣2p+1 = 0 si et seulement si x.. z. = 0n,1 ce qui est impossible puisque, R comme x.. est supposé de rang p + 1, on a d’après le Théorème du rang que dim(Ker(x.. )) = 0, d’où Ker(x.. ) = (0p+1,1 ), or z. ≠ 0p+1,1 par hypothèse donc ∣∣x.. z. ∣∣2 Rp+1 33/97 ≠0 Remarque : comme ∂β0 C (β̂0 , β̂1 ,... , β̂p ) = 0, on a n n i=1 i=1 ⇔ ∑ yi = ∑( β̂0 + β̂1 xi1 +... + β̂p xip ) ⇔ y = β̂0 + β̂1 x.1 +... + β̂p x.p ⇔ y = β̂0 + β̂1 x.1 +... + β̂p x.p autrement dit le centre de gravité (x.1 ,... , x.p , y. ) appartient à l’hyperplan affine solution au sens des moindres carrés. 34/97 Application numérique & graphique : cas p = 1 ∗ avec les formules explicites : hatbeta1 = np.cov(x, y, ddof=0)[0, 1] / np.var(x) hatbeta0 = np.mean(y) - hatbeta1 * np.mean(x) ∗ avec lstsq (LeaST SQuares) de scipy.linalg : n = len(x) xaugm = np.ones((n,2)) xaugm[:,1] = x hatbeta = linalg.lstsq(xaugm, y) ∗ avec LinearRegression de sklearn.linearmodel : regmod = linearmodel.LinearRegression().fit(xaugm, y) hatbeta0 = regmod.intercept hatbeta1 = regmod.coef ∗ avec OLS (Ordinary Least Squares) de statsmodels : xaugm = statsmodels.api.add− constant(x) modelelin = statsmodels.api.OLS(y, xaugm) reg = modelelin.fit() hatbeta = reg.params 35/97 On trouve pour l’exemple des appartements : β̂0 ≃ 33.64381565 ≃ 33.644 β̂1 ≃ 3.84782015 ≃ 3.848 c-à-d : prix ≃ 33.644 + 3.848 × surface (en noir : la droite y = β̂0 + β̂1 x ; en bleu : les appartements ; en rouge : un point de repère, le 14 ème appartement (70, 325) ) Application numérique & graphique : cas p = 2 x− tab = nba− tab[["taille", "age"]] n = nba.index.size x− tab["intercept"] = np.ones((n,1)) y− vect = nba− tab["poids"] ∗ avec les formules matricielles explicites : hat− beta = np.dot(linalg.inv(np.dot(np.transpose(x− tab), x− tab)), np.dot(np.transpose(x− tab), y− vect)) ∗ avec lstsq de scipy.linalg : hat− beta = linalg.lstsq(x− tab, y− vect) ∗ avec sklearn.linear− model.LinearRegression : reg− mod = LinearRegression().fit(x− tab, y− vect) hat− beta− 0 = reg− mod.intercept− hat− beta− 1, hat− beta− 2 = reg− mod.coef− [0:2] ∗ avec OLS (Ordinary Least Squares de statsmodels) : modelelin = statsmodels.api.OLS(y− vect, x− tab) reg = modelelin.fit() hat− beta = reg.params 36/97 y = β̂0 + β̂1 x1 + β̂2 x2 où β̂0 ≃ −296.599, β̂1 ≃ 6.327, β̂2 ≃ 0.651 c-à-d : poids ≃ −296.599 + 6.327 ∗ taille + 0.651 ∗ âge Application numérique & graphique : cas p = 5 matieres = tab− notes.columns[0:5] x− tab = tab− notes[matieres] x− tab[’intercept’] = np.ones((n,1)) y− vect = tab− notes[’Moyenne− ponderee’] hat− beta = np.dot(linalg.inv(np.dot(np.transpose(x− tab), x− tab)) , np.dot(np.transpose(x− tab), y− vect) ) ce qui donne : [0.3938768 0.30034548 0.1089102 0.10380356 0.08649371 1.5899103 ] ou avec sklearn.linear− model.LinearRegression matieres = tab− notes.columns[0:5] x− tab = tab− notes[matieres] reg = LinearRegression().fit(x− tab, y− vect) print(reg.coef− ) ce qui donne : [0.3938768 0.30034548 0.1089102 0.10380356 0.08649371 0. ] print(reg.intercept− ) ce qui donne : 1.5899102986314961 Rq : les vrais coefficients étaient [0.4, 0.3, 0.1, 0.1, 0.1] et 1.5 37/97 Deux questions supplémentaires : Comment utiliser ces résultats pour faire des prévisions ? Comment évaluer l’adéquation du modèle aux données ? 38/97 Prévision de l’étiquette d’une nouvelle donnée Supposons que l’on dispose d’une nouvelle donnée, non étiquetée, (xn+1,1 ,... , xn+1,p ) , mais qu’a priori l’on peut supposer son comportement "proche" de celui des données déjà connues. p Autrement dit, l’étiquette yn+1 de xn+1 doit vérifier yn+1 ≃ β0 + ∑ βj xn+1,j j=1 On utilise alors comme prévision de yn+1 la valeur p prev ∶= β̂0 + ∑β̂j xn+1, j = xn+1,. β̂. avec xn+1,. ∶= (1, xn+1,1 ,... , xn+1,p ). ŷn+1 j=1 Ex pr p = 2 (basketteurs) : comme β̂0 ≃ −296.599, β̂1 ≃ 6.327, β̂2 ≃ 0.651, prev prev pour un nouveau basketteur, de taille x(n+1),1 = 88 et d’âge x(n+1),2 = 20, on peut supposer qu’il pèsera prev prev prev ŷn+1 ∶= β̂0 + β̂1 x(n+1),1 + β̂2 x(n+1),2 ≃ −296.599 + 6.327 ∗ 88 + 0.651 ∗ 20 ≃ 273 avec OLS taille− new, age− new = 88, 20 x− new = np.array([1, taille− new, age− new]) y− prev = reg.predict(x− new) 39/97 Cas de la RLS : supposons que l’on dispose d’une nouvelle donnée, non étiquetée, xn+1 , mais qu’a priori l’on peut supposer son comportement "proche" de celui des données déjà connues. Autrement dit, l’étiquette yn+1 de xn+1 doit vérifier yn+1 ≃ β0 + β1 xn+1 On utilise alors comme prévision de yn+1 la valeur prev ŷn+1 ∶= β̂0 + β̂1 xn+1 Ex pr p = 1 (appartements) : en postulant que le nouvel appartement, de surface xn+1 = 77, que notre agent prend en charge suit la même tendance que les appartements du même quartier qu’il a déjà vendus, autrement dit, que son étiquette yn+1 vérifie yn+1 ≃ β0 + β1 xn+1 , notre agent peut fixer le prix de vente de ce nouvel appartement à prev ŷn+1 ∶= β̂0 + β̂1 xn+1 ≃ 33.644 + 3.848 ∗ 77 ≃ 330 (point orange) surface21 = 77 yprev21 = hatbeta0 + hatbeta1 * surface21 ou bien, en utilisant OLS surface21augm = np.array([1, surface21]) yprev21 = reg.predict(surface21augm) ou encore, en utilisant les notations matricielles yprev21 = np.dot(surface21augm, hatbetam) 40/97 Evaluation de l’adéquation du modèle ANOVA (Analyse de la Variabilité) & R 2 n Somme des Carrés Totale sct ∶= ∣∣y. − y. 1n,1 ∣∣2Rn = ∑ (yi − y. )2 = nσ̂y2 i=1 n Somme des Carrés Expliquée sce ∶= ∣∣ŷ. − y. 1n,1 ∣∣2Rn = ∑ (ŷi − y. )2 i=1 n n i=1 i=1 Somme des Carrés Résiduelle scr ∶= ∣∣y. − ŷ. ∣∣2Rn = ∑ (yi − ŷi )2 = ∑ ϵ̂2i Rq : d’après le Théorème de Pythagore, sct = sce + scr (exercice !) scr Coefficient de détermination r 2 ∶= sce sct = 1 − sct ↝ r 2 ∈ [0, 1] correspond à la proportion de variabilité des yi expliquée p par le modèle linéaire y = β0 + ∑ βj xj j=1 Plus r 2 est proche de 1, plus le modèle est adapté aux données Remarque : lorsque p = 1, on a r 2 = ρ̂2x ,y 41/97 Exemple pour p = 1 (appartements) : en noir : la droite y = β̂0 + β̂1 x ; hatprix = hatbeta0 + hatbeta1 * np.array(surface) hatepsilon = np.array(prix) - hatprix scr = np.sum(hatepsilon**2) sct = n * np.var(prix) sce = np.sum( ( hatprix - np.mean(prix) )**2) r2 = sce/sct points bleus : les données appartements initiales (xi , yi ) ; points oranges : les prédictions i.e. les couples (xi , ŷi ) en orange pointillé : les longueurs des résidus point rouge : le 14 ème appartement (70, 325). ou bien, avec la fonction OLS : hatprix = reg.fittedvalues hatepsilon = reg.resid scr = reg.ssr sce = reg.ess sct = reg.centered− tss r2 = reg.rsquared ou encore, avec les notations matricielles hatprixm = np.dot(surfaceaugm, hatbetam) hatepsilonm = np.array(prix) - hatprixm scr = np.dot(np.transpose(hatepsilonm), hatepsilonm) 42/97 ce qui donne pour la somme des carrés des résidus : scr ≃ 36476.88 pour la somme des carrés expliquée : sce ≃ 195068.321 pour la somme des carrés totale : sct ≃ 231545.200 pour le coefficient de détermination : r 2 ≃ 0.842 Exemple pour p = 2 (basketteurs) : avec la fonction OLS : hatprix = reg.fittedvalues hatepsilon = reg.resid scr = reg.ssr sce = reg.ess sct = reg.centered− tss r2 = reg.rsquared en rouge : le plan y = β̂0 + β̂1 x1 + β̂2 x2 ; points bleus : les données "joueurs de NBA" initiales (xi1 , xi2 , yi ) ; points oranges : les prédictions i.e. les triplets (xi1 , xi2 , ŷi ) en orange pointillé : les longueurs des résidus ou encore, avec les notations matricielles : hat− y = np.dot(x− tab, hat− beta) hat− epsilon = y− vect - hat− y scr = np.dot(np.transpose(hat− epsilon), hat− epsilon) ce qui donne pour la somme des carrés des résidus : scr ≃ 112784.883 pour la somme des carrés expliquée : sce ≃ 244982.206 pour la somme des carrés totale : sct ≃ 357767.089 pour le coefficient de détermination : r 2 ≃ 0.685 43/97 Vers la modélisation probabiliste des erreurs Etant donné que nos points ne sont PAS sur un hyperplan affine, nous avons un terme d’erreur incompressible. Nous allons modéliser ce terme d’erreur par un vecteur aléatoire. Pourquoi modéliser de manière probabiliste ce terme d’erreur ? Parce que ces erreurs peuvent provenir de causes inconnues, qui peuvent être en très grand nombre, et dont on n’a aucune idée, d’une part intrinsèque d’aléatoire dans la situation considérée 44/97 B) Modélisation probabiliste des erreurs Nous allons donc à présent supposer que la fluctuation de y autour de p l’hyperplan affine de régression linéaire théorique y = β0 + ∑ βj xj (resp. j=1 est de nature aléatoire : autrement dit, nous supposons que y est une réalisation d’une v.a.r. Y où y = β0 + β1 x pour le cas de la RLS p = 1) p Y = β0 + ∑βj xj + j=1 (resp. Y = β0 + β1 x + pour la RLS) Y est la variable à expliquer, à valeurs dans R x1 ,... , xp sont les variables explicatives, chacune appartenant à R (resp. pr la RLS, x est la variable explicative, appartenant à R) est le terme d’erreur aléatoire du modèle, à valeurs dans R ; β0 , β1 ,... , βp sont les (p + 1) paramètres inconnus à estimer, chacun appartenant à R (resp. pr la RLS, β0 et β1 sont les deux paramètres inconnus à estimer, tous deux appartenant à R) 45/97 Nous supposons donc que les données ((x1. , y1 ),... , (xn. , yn )) ((x1 , y1 ),... , (xn , yn )) pr la RLS) sont une réalisation de l’observation ((X1. , Y1 ),... , (Xn. , Yn )) (resp. ((X1 , Y1 ),... , (Xn , Yn )) pr la RLS) où : X1. ,. , Xn. vecteurs aléat. constants : ∀i ∈ J1, nK, Xi. = xi. = (resp. ⎛xi1 ⎞ ⎜ ⋮ ⎟ ⎝xip ⎠ (resp. pr la RLS : X1 ,... , Xn v.a.r. constantes : ∀i ∈ J1, nK, Xi = xi ) ∀i ∈ J1, nK, Yi = β0 + ⟨ (resp. pr la RLS : ∀i ∈ J1, nK, Yi = β0 + β1 xi + p , xi. ⟩Rp + i = β0 + ∑ βj xij + i ⎛β1 ⎞ ⎜ ⋮ ⎟ ⎝βp ⎠ i ) j=1 1 ,... , n v.a.r. représentant les termes d’erreurs et vérifiant (A1) "les erreurs sont centrées" : ∀i ∈ J1, nK, E[i ] = 0 (A2) "condition d’homoscédasticité" : ∀i ∈ J1, nK, Var[i ] = σ 2 (A3) "les termes d’erreurs sont non corrélés" : ∀i, i ′ ∈ J1, nK, i ≠ i ′ ⇒ Cov(i , i ′ ) = 0 β0 , β1 ,. , βp ∈ R, σ 2 ∈ R⋆+ (resp. β0 , β1 ∈ R, σ 2 ∈ R⋆ + pr la RLS) paramètres inconnus L’observation est donc ((x1. , Y1 ),... , (xn. , Yn )) (resp. ((x1 , Y1 ),... , (xn , Yn )) pr la RLS), avec les xi. déterministes ds Rp (resp. xi déterministes ds R pr la RLS) , et les Yi aléatoires à valeurs ds R. 46/97 On peut également réécrire matriciellement le modèle p ∀i ∈ J1, nK, Yi = β0 + ∑ βj xij + i j=1 ⎛Y1 ⎞ Y. = ⎜ ⋮ ⎟ ∈ Mn,1 (R) ↝ ⎝Yn ⎠ sous la forme Y. = x.. β. + . où vecteur à expliquer , (ALÉATOIRE & OBSERVÉ) ⎛1 x11... x1p ⎞ ⋮ ⋮ ⎟ ∈ Mn,(p+1) (R) ↝ x.. = ⎜ ⋮ ⎝1 xn1... xnp ⎠ β. = ⎛ β0 ⎞ ⎜ β1 ⎟ ⎜ ⎟ ⎜ ⋮ ⎟ ⎝βp ⎠ ∈ M(p+1),1 (R) ↝ matrice explicative , (DÉTERMINISTE & OBSERVÉE) , vecteur des paramètres (DÉTERMINISTE & INCONNU (NON OBSERVÉ)) ⎛1 ⎞ . = ⎜ ⋮ ⎟ ∈ Mn,1 (R) ↝ ⎝n ⎠ vecteur d’erreurs. (ALÉATOIRE & NON OBSERVÉ) Les conditions (A1)-(A2)-(A3) se réécrivent matriciellement : (A1) : E[. ] = 0n,1 et (A2)-(A3) : Var[. ] = σ 2 In (Ici on ne suppose PAS que les 47/97 i sont indépendants, NI qu’ils sont de même loi, mais on le supposera ds le C)) Rappels : E[. ] = ⎜ ⎛E[ 1 ]⎞ ⋮ ⎟ ⎝E[ n ]⎠ . est le vecteur d’espérance du vecteur aléat. Var[. ] est la matrice de variance-covariance du vecteur aléat. et est définie par Var[. ] = [Cov(i , k )]1≤i,k≤n = ⎛ Var[ 1 ] ⎜Cov( 2 , 1 ) ⎜ ⋮ ⎜ ⎜ ⎜ ⎝Cov( n , 1 ) Cov( 1 , 2 ) Var[ 2 ] n )...... Cov( 1 , Var[ n−1 ] Cov( n , n−1 ) ⎞ ⎟ ⎟ ⎟ ⎟ Cov( n−1 , n )⎟ ⎠ Var[ n ] Conséquences des conditions (A1)-(A2)-(A3) : p E[Y. ] = x.. β. c-à-d ∀i ∈ J1, nK, E[Yi ] = β0 + ∑ βj xij j=1 p Cela entraîne : E[Y. ] = β0 + ∑ βj x.j j=1 Var[Y. ] = Cov(Y. , . ) = σ 2 In ; c-à-d ∀i, i ′ ∈ J1, nK, ↝ si i = i ′ , Cov(Yi , Yi ′ ) = Var[Yi ] = σ 2 , ↝ si i ≠ i ′ , Cov(Yi , Yi ′ ) = Cov(Yi , i ′ ) = 0 Cela entraîne : ∀i ∈ J1, nK, Var[Y. ] = Cov(Yi , Y. ) = σ 2 /n Exercice : prouvez ces égalités en justifiant chaque étape de calcul 48/97 . ⋮ Estimation des paramètres théoriques β0, β1,. , βp D’après le A), on construit à partir de l’observation ((x1. , Y1 ),... , (xn. , Yn )) les estimateurs aléatoires de β0 , β1 ,. , βp suivants : ⎛βˆ 0 ⎞ ⎛βˆ 0 ((x1. , Y1 ),... , (xn. , Yn ))⎞ ⎜βˆ ⎟ ⎜βˆ ((x , Y ),... , (x , Y ))⎟ ⎜ ⎟ ⎜ n. n ⎟ ⎟ = (t x.. x.. )−1 t x.. Y. βˆ. = ⎜ 1 ⎟ = ⎜ 1 1. 1 ⎟ ⎜ ⋮ ⎟ ⎜ ⋮ ⎟ ⎜ ⎟ ⎜ ⎝βˆ p ⎠ ⎝βˆ p ((x1. , Y1 ),... , (xn. , Yn ))⎠ Pr une réalisation ((x1. , y1 ),... , (xn. , yn )) de l’observation, on retrouve les valeurs ⎛β̂0 ⎞ ⎛β̂0 ((x1. , y1 ),... , (xn. , yn ))⎞ ⎜β̂1 ⎟ ⎜β̂1 ((x1. , y1 ),... , (xn. , yn ))⎟ t ⎟ ⎜ ⎟ = ( x.. x.. )−1 t x.. y. β̂. = ⎜ ⎜ ⋮ ⎟=⎜ ⎟ ⋮ ⎜ ⎟ ⎜ ⎟ ⎝β̂p ⎠ ⎝β̂p ((x1. , y1 ),... , (xn. , yn ))⎠ 49/97 Pour la RLS, d’après le A), on construit à partir de l’observation ((x1 , Y1 ),... , (xn , Yn )) les estimateurs aléatoires de β0 , β1 suivants : Y. − β̂1 x. ⎛ ⎞ n ˆ ⎞ ⎛β ˆ ((x , Y ),... , (x , Y ))⎞ ⎜ 1 ∑ ⎟ ⎛β (xi −x. )(Yi −Y. ) ⎟ ⎜ n n n 1 1 ˆ 0 0 ĉx ,Y ⎟ = ⎜ i=1 β. = ˆ = ˆ σ̂Y... ⎜ = σ̂ ρ̂x. ,Y. ⎟ = ⎝β ⎠ ⎝β ((x1 , Y1 ),... , (xn , Yn ))⎠ ⎜ n ⎟ 2 x. 1 1 σ̂x 1 ∑ (x −x )2. i. ⎝ ⎠ n i=1 et pr une réalisation ((x1 , y1 ),... , (xn , yn )) de l’observation, on retrouve les valeurs y. − β̂1 x. ⎞ ⎛ n ⎟ ⎜ 1 ∑ (x −x. )(yi −y. ) β̂ β̂ ((x1. , y1 ),... , (xn. , yn )) ⎟ ⎜ ĉx. ,y. σ̂y. ⎟ β̂. = ( 0 ) = ( 0 ) = ⎜ n i=1 i ⎟ ⎜ = = ρ̂ β̂1 β̂1 ((x1. , y1 ),... , (xn. , yn )) x ,y..⎟ n ⎜ σ̂x. σ̂x2 1 ∑ (x −x )2. i. ⎠ ⎝ n i=1 50/97 Propriétés de l’estimateur βˆ. de β. Propriétés de l’estimateur βˆ. de β. (A SAVOIR PAR COEUR) Sous les conditions (A0)-(A1)-(A2)-(A3) : pour tous (β. , σ 2 ) ∈ Rp+1 × R⋆+ E(β ,σ2 ) [βˆ. ] = β. ↝ βˆ. estimateur sans biais de β. ∀j ∈ J0, pK, E(β ,σ2 ) [βˆ j ] = βj ↝ βˆ j estimateur sans biais de βj.. Rq : reste vrai pr la RLS avec p = 1 Var(β ,σ2 ) [βˆ. ] = σ 2 (t x.. x.. )−1. Rq : pr la RLS, cela donne Var d’où Var (β. ,σ 2 ) (β. ,σ 2 ) (x )2 ˆ ] = σ2 ( 1 + [β n 0 nσ̂x2 2 2 ˆ ] = σ2 (σ̂x + (x ) [β. nσ̂x2 −x ) , Var (β. ,σ 2 ) ˆ ]= [β 1 σ2 nσ̂x2 −x ) 1 et Cov (β. ,σ 2 ) ˆ ,β ˆ ) = − σ2 x (β 2 0 1 nσ̂x Propriété de Gauss-Markov (ADMIS, pas à savoir) : βˆ. est de variance minimale parmi les estimateurs sans biais de β. linéaires en Y1 ,... , Yn Ds la suite, on omettra l’indice pr les IC et les tests 51/97 (β. ,σ 2 ) ; on reprendra cette notation plus loin Valeurs prédites & Résidus aléatoires : définition p ∀i ∈ J1, nK, Ŷi ∶= xi.βˆ. = βˆ 0 + ∑ βˆ j xij v.a.r. dite valeur prédite pour Yi j=1 ↝ Ŷi = Ŷi ((x1. , Y1 ),... , (xn. , Yn )) ; pr une réal. ((x1. , y1 ),... , (xn. , yn )) de p l’obs. ((x1. , Y1 ),... , (xn. , Yn )), Ŷi ((x1. , y1 ),... , (xn. , yn )) = β̂0 + ∑ β̂j xij = ŷi j=1 ˆ i ∶= Yi − Ŷi v.a.r. dite résidu ds la prédiction de Yi par Ŷi ∀i ∈ J1, nK, ˆi = ˆ i ((x1. , Y1 ),... , (xn. , Yn )) ; pr une réal. ((x1. , y1 ),... , (xn. , yn )) de ↝ ˆ i ((x1. , y1 ),... , (xn. , yn )) = yi − ŷi = ϵ̂i l’obs. ((x1. , Y1 ),... , (xn. , Yn )), p p n n n 2 2 2 ˆ i v.a.r. min ∑ (Yi − (β0 + ∑ βj xij )) = ∑ (Yi − (βˆ 0 + ∑ βˆ j xij )) = ∑ p+1 β. ∈R i=1 j=1 i=1 j=1 i=1 dite somme des carrés des résidus (SCR) ↝ SCR = SCR((x1. , Y1 ),... , (xn. , Yn )) et pr une réal. n ((x1. , y1 ),... , (xn. , yn )), SCR((x1. , y1 ),... , (xn. , yn )) = ∑ ϵ̂2i = scr i=1 52/97 Valeurs prédites & Résidus aléatoires : définition pour la RLS ˆ x v.a.r. dite valeur prédite pour Y ˆ =β ˆ +β ∀i ∈ J1, nK, Ŷi ∶= xi. β i. 0 1 i ↝ Ŷi = Ŷi ((x1 , Y1 ),... , (xn , Yn )) ; pr une réal. ((x1 , y1 ),... , (xn , yn )) de l’obs. ((x1 , Y1 ),... , (xn , Yn )), Ŷi ((x1 , y1 ),... , (xn , yn )) = β̂0 + β̂1 xi = ŷi ∀i ∈ J1, nK, ˆ i ∶= Yi − Ŷi v.a.r. dite résidu ds la prédiction de Yi par Ŷi ↝ ˆ i = ˆ i ((x1 , Y1 ),... , (xn , Yn )) ; pr une réal. ((x1 , y1 ),... , (xn , yn )) de l’obs. ((x1 , Y1 ),... , (xn , Yn )), ˆ i ((x1 , y1 ),... , (xn , yn )) = yi − ŷi = ϵ̂i min n n n 2 ˆ +β ˆ x ))2 = ∑ ˆ 2 v.a.r. dite somme des carrés des résidus (SCR) ∑ (Yi − (β0 + β1 xi )) = ∑ (Yi − (β i 0 1 i i=1 i=1 β. ∈Rp+1 i=1 n ↝ SCR = SCR((x1 , Y1 ),... , (xn , Yn )) et pr une réal. ((x1 , y1 ),... , (xn , yn )), SCR((x1 , y1 ),... , (xn , yn )) = ∑ ϵ̂2i = scr i=1 53/97 Valeurs prédites & Résidus aléatoires : propriétés Avec les notations matricielles Ŷ. ∶= ⎛Ŷ1 ⎞ ⎜ ⋮ ⎟ ⎝Ŷn ⎠ ˆ ˆ. ∶= ⎛⎜⋮1 ⎞⎟ , on a les propriétés : , Ŷ. = x..βˆ. = x.. (t x.. x.. )−1 t x.. Y. = ΠIm(x.. ) (Y. ) ⎝ ˆ n⎠ E[Ŷ. ] = x.. β. = E[Y. ] ↝ Ŷ. estimateur sans biais E[Y. ] p ∀i ∈ J1, nK, E[Ŷi ] = β0 + ∑ βj xij = E[Yi ] j=1 ↝ Ŷi estimateur sans biais de E[Yi ] Var[Ŷ. ] = σ 2 x.. (t x.. x.. )−1 t x.. ˆ. = (In − x.. (t x.. x.. )−1 t x.. )Y. = Π ⊥ (Y. ) (Im(x.. )) ⊥ (. ) = (In − x.. (t x.. x.. )−1 t x.. ). = Π (Im(x.. )) ˆ. ] = 0n,1 et Var[ ˆ. ] = σ 2 (In − x.. (t x.. x.. )−1 t x.. ) E[ ˆ. ) = 0n,n Cov(Ŷ. , Exercice : prouvez ces égalités 54/97 Valeurs prédites & Résidus aléatoires : propriétés pour la RLS E[Ŷi ] = β0 + β1 xi = E[Yi ] ↝ Ŷi est un estimateur sans biais de E[Yi ] ; ↝ E[ ˆ i ] = 0 car ˆ i = Yi − Ŷi par déf Var[Ŷi ] = σ 2 ( n1 + E[ n1 n (xi −x )2 nσ̂x2 ) et (x −x )2 Var[ ˆ i ] = σ 2 (1 − n1 − i 2 ) nσ̂x 2 σ2 ∑ ˆ i ] = n−2 n i=1 ′ ∀i, i ∈ J1, nK, Cov(Ŷi , ˆ i ′ ) = 0 On peut également montrer les propriétés secondaires suivantes : n E[ n1 ∑ Ŷi ] = β0 + β1 x = E[Y. ] n et i=1 2 Var[ n1 ∑ Ŷi ] = σn i=1 2 Cov(Ŷi , Ŷi ′ ) = Cov(Yi , Ŷi ′ ) = σ 2 (σ̂x2 + (xi − x )(xi ′ − x )) nσ̂ x 2 ↝ Cov( ˆ i , ˆ i ′ ) = σ 2 δii ′ − σ 2 (σ̂x2 + (xi − x )(xi ′ − x )) nσ̂x ′ (où δii ′ = 1 si i = i et 0 sinon) Rq : on a utilisé pour cela les propriétés intermédiaires suivantes : ˆ ) = σ2 (σ̂ 2 − x (x − x )) Cov(Yi , β i x 0 nσ̂ 2 ˆ )= Cov(Yi , β 1 55/97 x σ2 nσ̂x2 (xi − x ) et et ˆ ) = σ2 ; Cov(Y. , β n 0 ˆ )=0 Cov(Y. , β 1 Estimation du paramètre σ 2 Par déf, σ 2 = Var[i ] = Var[Yi ] = E[(Yi − E[Yi ])2 ]. Comme Ŷi estimateur de E[Yi ], on estime σ 2 par une quantité proportionnelle à Ŝ 2 = Ŝ 2 ((x1. , Y1 ),. , (xn. , Yn )) ∶= = 1 n n ∑ (Yi − Ŷi )2 : i=1 n 1 n 2 ∑(Yi − Ŷi ) n − (p + 1) n i=1 n SCR 1 ˆ 2i = ∑ n − (p + 1) i=1 n − (p + 1) Pr une réal. ((x1. , y1 ),. , (xn. , yn )) de l’obs., ŝ 2 ∶= Ŝ 2 ((x1. , y1 ),. , (xn. , yn )) = 1 n−(p+1) n ∑ ϵ̂2i = i=1 scr n−(p+1) 2 n 1 Cas de la RLS : Ŝ 2 = Ŝ 2 ((x1 , Y1 ),. , (xn , Yn )) ∶= n−2 et pr une réal. ((x1 , y1 ),. , (xn , yn )) de l’obs., ∑i=1 ˆ i = SCR n−2 n 2 1 scr ŝ 2 ∶= Ŝ 2 ((x1 , y1 ),. , (xn , yn )) = n−2 ∑ ϵ̂i = n−2 i=1 E[Ŝ 2 ] = σ 2 ↝ Ŝ 2 estimateur sans biais de σ 2 Rq : la perte de (p + 1) degrés de liberté dans l’expression de Ŝ 2 est le "coût" de l’estimation de β. ) 56/97 Prévision de l’étiquette d’une nouvelle donnée Supposons que l’on dispose d’une nouvelle donnée, non étiquetée, (xn+1,1 ,... , xn+1,p ) , mais qu’a priori l’on peut supposer son comportement "proche" de celui des données déjà connues. p Autrement dit, l’étiquette yn+1 de xn+1 doit vérifier yn+1 ≃ β0 + ∑ βj xn+1,j j=1 On utilise alors comme prévision de yn+1 la valeur p prev ∶= β̂0 + ∑βj xn+1, j = xn+1,. β̂. avec xn+1,. ∶= (1, xn+1,1 ,... , xn+1,p ). ŷn+1 j=1 Si l’on suppose que yn+1 est, comme les étiquettes précédentes, une p réalisation de Yn+1 = β0 + ∑ βj xn+1,j + n+1 = xn+1,. β. + n+1 , on définit j=1 p prev la prévision de Yn+1 par Ŷn+1 ∶= βˆ 0 + ∑βˆ j xn+1,j = xn+1,.βˆ. j=1 prev prev ˆ prev et l’erreur dans la prévision de Yn+1 par Ŷn+1 par n+1 ∶= Yn+1 − Ŷn+1 57/97 prev ˆ prev et Propriétés de Ŷn+1 n+1 p prev ] = β0 + ∑ βj xn+1, j = E[Yn+1 ] Si E[n+1 ] = 0, alors E[Ŷn+1 j=1 ↝ prev Ŷn+1 ˆ n+1 ] = 0) estimateur sans biais de E[Yn+1 ] (d’où E[ prev prev Cas de la RLS : E[Ŷn+1 ] = β0 + β1 xn+1 = E[Yn+1 ] Si de plus cov(n+1 , i ′ ) = { 0 si i ′ ∈ J1, nK , alors σ 2 si i ′ = n + 1 prev Var[Ŷn+1 ] = σ 2 xn+1,. (t x.. x.. )−1 t xn+1,. prev (x −x )2 Cas de la RLS : Var[Ŷn+1 ] = σ 2 ( n1 + n+1 2 ) nσ̂x Par conséquent, la variance augmente quand xn+1 s’éloigne du centre de gravité x des étiquettes des points. Effectuer une prévision de yn+1 lorsque xn+1 est loin de x est risqué, car dans cette situation, la variance de l’erreur de prévision peut être très grande... ˆ n+1 ] = σ 2 + Var[Ŷ prev ] Var[ n+1 prev (x −x )2 prev Cas de la RLS : Var[ ˆ n+1 ] = σ 2 (1 + n1 + n+1 2 ) dépend nσ̂x 2 de la variabilité intrinsèque σ de Yn+1 ↝ cette source de variabilité ne peut pas être réduite de la variabilité due à l’imprécision des estimations de β0 , β1 ↝ cette source de variabilité peut être réduite en augmentant le nombre n de données. 58/97 En supposant que les erreurs entre le modèle théorique et nos données sont centrées, de même variance et deux à deux non corrélées, on a déterminé l’espérance et la variance des prédictions et des prévisions. On a donc obtenu des informations supplémentaires par rapport à la partie A) où l’on obtenait juste des valeurs numériques ponctuelles pour les prédictions et les prévisions. Comment faire pour quantifier encore plus précisément le risque pris en déviant légèrement de ces valeurs numériques ? Nous avons besoin de leurs lois, et pour cela également de la loi de βˆ. et de Ŝ 2 ! 59/97 C) Modèle linéaire gaussien Nous nous plaçons enfin sous les conditions (A0) et (A4) : . ∼ Nn (0n,1 , σ2 In ), ou encore 1 ,. , n i.i.d. ∼ N (0, σ 2 ) qui équivaut à Y. ∼ Nn (x.. β. , σ 2 In ) ou encore à p Y1 ,. , Yn indépendants, avec ∀i ∈ J1, nK, Yi ∼ N (β0 + ∑βj xij , σ 2 ) j=1 Cas de la RLS : cette condition équivaut à Y1 ,. , Yn indépendants, avec ∀i ∈ J1, nK, Yi ∼ N (β0 + β1 xi , σ 2 ) Remarque : (A4) implique (A1)-(A2)-(A3) 60/97 C) 1. IC & Test de significativité sur βˆ j , j ∈ J0, pK Propriétés supplémentaires de βˆ. , Ŝ 2 Sous les conditions (A0) & (A4), βˆ. ∼ Np+1 (β. , σ 2 (t x.. x.. )−1 ) ∀j ∈ J0, pK, βˆ j ∼ N (βj , σ 2 ((t x.. x.. )−1 )(j+1),(j+1) ) ; (n−(p+1))Ŝ 2 σ2 βˆ. á ∼ χ2 (n − (p + 1)) (n−(p+1))Ŝ 2 σ2 ˆ β j −βj 61/97 ¿ Á Á Àσ 2 (t x x )−1.... ∀j ∈ J0, pK, √ ( )(j+1),(j+1) (n−(p+1))Ŝ 2 /(n−(p+1)) σ2 ˆ = √ Ŝ β j −βj ((t x.. x.. )−1 )(j+1),(j+1) ∼ t(n − (p + 1)) Cas de la RLS : Sous la condition (A4), ˆ ∼ N (β , σ 2 ( 1 + (x )2 )) β 0 n 0 nσ̂x2 (n−2)Ŝ 2 σ2 ˆ ∼ N (β , σ2 ) et β 1 nσ̂ 2 1 x ∼ χ2 (n − 2) ˆ ,β ˆ ) á Ŝ 2 (β 0 1 βˆ −β0 ¿0 Á 1 (x )2 À + Ŝ Á n nσ̂ 2 x βˆ 1 −β1 √ Ŝ nσ̂x2 √ ∼ t(n − 2) ( ↝ Ŝ ∼ t(n − 2) ( ↝ √ Ŝ nσ̂x2 1 n + (x )2 nσ̂x2 est un estimateur de est un estimateur de Application numérique pour l’exemple des appartements : √ ŝ nσ̂x2 √ ˆ ] ) Var[β 0 √ ˆ ] ) Var[β 1 ≃ 0.392 s2 = scr / (n-2) estimateursdhatbeta1 = np.sqrt( s2 / ( n * np.var(surface) ) ) ou bien, à l’aide de la fonction OLS : estimateursdhatbeta = reg.bse matvar = scipy.linalg.inv(np.dot(np.transpose(surfaceaugm), surfaceaugm)) estimateursdhatbeta1 = np.sqrt(s2 * np.diag(matvar)) 62/97 Intervalle de confiance pour βj , j ∈ J0, pK Intervalle de confiance de niveau exact 1 − α de βj : √ (β ) t(n−(p+1)) ((t x.. x.. )−1 )(j+1),(j+1) ] IC1−αj ((x1. , Y1 ),. , (xn. , Yn )) ∶= [ β̂j ± q1− α Ŝ 2 (β ) ↝ ∀(β. , σ 2 ) ∈ Rp+1 × R⋆+ , P(β ,σ2 ) ( βj ∈ IC1−αj ((x1. , Y1 ),. , (xn. , Yn )) ) = 1 − α. ↝ pour une réal. ((x1. , y1 ),. , (xn. , yn )) de l’obs. ((x1. , Y1 ),. , (xn. , Yn )), (β ) (β ) ic1−αj ∶= IC1−αj ((x1. , y1 ),. , (xn. , yn )) √ t(n−(p+1)) ((t x.. x.. )−1 )(j+1),(j+1) ] ŝ = [ β̂j ± q1− α 2 (β1 ) ˆ ± q t(n−2) √ Ŝ ] ((x1 , Y1 ),. , (xn , Yn )) ∶= [β Cas de la RLS : IC1−α 1 1− α nσ̂x2 2 Ex pr p = 1 (appartements) : (β ) 1 ≃ [3.024, 4.672] avec alpha = 5% ic0.95 qalpha = scipy.stats.t.ppf(1-alpha/2, df = n-2) precision = qalpha * np.sqrt(s2) /np.sqrt(n*np.var(surface)) ic1moinsalpha = hatbeta1 + np.array([-1,1]) * precision 63/97 Ex pr p = 2 (basketteurs) : IC de niveau 1 − α = 95% avec OLS : reg.conf− int(alpha) variable taille age intercept borne inf 5.946877 0.347830 -327.648146 borne sup 6.706443 0.954364 -265.549080 Test de Student de significativité de βj , j ∈ J0, pK On procède au test bilatère de H0 ∶ βj = 0 contre H1 ∶ βj ≠ 0 Observation : ((x1. , Y1 ),. , (xn. , Yn )) ∈ (Rp × R)n Y. ∼ N (x.. β. , σ 2 In ) β0 , β1 ,. , βp ∈ R, σ ∈ 2 ⎛Y1 ⎞ ⎛1 ( rappel : Y. = ⎜ ⋮ ⎟, x.. = ⎜ ⋮ ⎝Yn ⎠ ⎝1 R⋆+ (noté (Y1 ,. , Yn )) x11 ⋮ xn1...... ⎛β0 ⎞ x1p ⎞ ⎜β ⎟ ⋮ ⎟, β. = ⎜ 1 ⎟ ) ⎜ ⋮ ⎟ xnp ⎠ ⎝βp ⎠ inconnus ˆ Statistique utilisée Tn = Tn (Y1 ,. , Yn ) ∶= √ Ŝ βj ((t x.. x.. )−1 )(j+1),(j+1) sous H0 ∶ βj = 0, Tn ∼ t(n − (p + 1)) sous H1 ∶ βj ≠ 0, ∣Tn ∣ aura tendance à prendre des valeurs plus grandes que sous H0 64/97 Test exact de taille α ∈]0, 1[ : φα (Y1 ,. , Yn ) ∶= 1∣T On a φα (Y1 ,. , Yn ) = 1 t(n−(p+1)) n ∣>q1− α 2 j ((x1. ,Y1 ),. ,(xn. ,Yn )) 0 ∈ IC1−α (β ) p-valeur associée à une réalisation (y1 ,. , yn ) : π(y1 ,. , yn ) = Pβj =0 (∣Tn (Y1 ,. , Yn )∣ ≥ ∣Tn (y1 ,. , yn )∣) = 2(1 − Ft(n−(p+1)) (∣Tn (y1 ,. , yn )∣)) Remarques : rejeter H0 signifie que βj est significativement non nul ; βj s’interprète comme le taux d’accroissement moyen de y en fonction d’une variation de xj à x1 ,. , xj−1 , xj+1 ,. , xp fixés Exemple pour p = 2 (basketteurs) : avec OLS, on trouve variable j = 1 taille j = 2 age j = 0 intercept 65/97 reg.tvalues statistique 32.729203 4.218099 -18.767675 reg.pvalues p-valeur 1.312927e-126 2.923456e-05 6.229313e-60 donc, pr j = 0, 1, 2, on rejette H0 au niveau α = 5% pr le test de significativé de βj Rq : si on ne regarde que les 50 premiers basketteurs, les p-valeurs sont resp. 0.461129, 0.132534, 0.354148 et donc, pr j = 0, 1, 2, on conserve H0 au niveau α = 5% pr le test de significativé de βj Cas de la RLS : Test de Student de significativité de β1 Remarque : cela revient à tester l’absence de liaison linéaire entre les xi et les yi , puisque celle-ci se traduit par la nullité du coefficient directeur β1. On procède donc au test bilatère de H0 ∶ β1 = 0 contre H1 ∶ β1 ≠ 0 Observation : ((x1 , Y1 ),. , (xn , Yn )) ∈ R2n (notation abrégée (Y1 ,. , Yn )) ∀i ∈ J1, nK, Yi ∼ N (β0 + β1 xi , σ 2 ) et Y1 ,... , Yn á β0 , β1 ∈ R, σ 2 ∈ R⋆ + inconnus Statistique utilisée T (Y1 ,. , Yn ) ∶= βˆ 1 √ Ŝ nσ̂x2 sous H0 ∶ β1 = 0, T (Y1 ,. , Yn ) ∼ t(n − 2) sous H1 ∶ β1 ≠ 0, ∣T (Y1 ,. , Yn )∣ aura tendance à prendre des valeurs plus grandes que sous H0 Test statistique exact, dit de Student, de taille α ∈]0, 1[ : φα (Y1 ,. , Yn ) ∶= 1 t(n−2) ∣T (Y1 ,. ,Yn )∣>q 1− α 2 Application numérique dans notre exemple : scr ≃ 36476.879 et ŝ 2 ≃ 2026.493, d’où T (y1 ,. , yn ) ≃ 9.81 ; t(18) au niveau α = 5%, on a donc φ0.05 (y1 ,. , yn ) = 0 car T (y1 ,. , yn ) ≃ 9.81 > q0.975 ≃ 2.101 et l’on rejette donc H0 ∶ β1 = 0 au niveau 5% p-valeur associée à une réalisation (y1 ,. , yn ) : π(y1 ,. , yn ) = Pβ1 =0 (∣T (Y1 ,. , Yn )∣ ≥ ∣T (y1 ,. , yn )∣) = 2(1 − Ft(n−2) (∣T (y1 ,. , yn )∣)) Application numérique dans notre exemple : π(y1 ,. , yn ) ≃ 2(1 − Ft(18) (∣9.81∣)) ≃ 1.197e − 08 qui est très petit, la réalisation de T (Y1 ,. Yn ) que l’on a obtenue serait donc très aberrante sous H0 , et l’on est donc très confiant dans notre décision de rejeter H0. 66/97 n = len(surface) statStudent = hatbeta1 / ( np.sqrt(s2) / ( np.sqrt(n) * np.std(surface) ) ) alpha = 0.05 quantilealpha = scipy.stats.t.ppf(1 - alpha/2, df = n-2) phialpha = ( np.abs(statStudent) > quantilealpha ) if phialpha : print("On rejette H0 au niveau ", 100*alpha, " else : print("On conserve H0 au niveau ", 100*alpha, " pvalStudent = 2 * ( 1 - scipy.stats.t.cdf(np.abs(statStudent), df = n-2) ) ou bien, avec la fonction OLS statStudent = reg.tvalues pvalStudent = reg.pvalues 67/97 C)2. IP de Yn+1 & IC de E[Yn+1] Pr (xn+1,1 ,. , xn+1,p ) nvelle donnée non étiquetée, d’étiquette modélisée par p Yn+1 = β0 + ∑ βj xn+1, j + n+1 , où j=1 n+1 centrée et décorrélée de 1 ,. , n , de p prev prev prev ]. ] et Var[Ŷn+1 ∶= βˆ 0 + ∑ βˆ j xn+1,j on a déterminé E[Ŷn+1 prévision Ŷn+1 j=1 Si ∼ N (0, σ 2 ), avec xn+1,. ∶= (1, xn+1,1 ,. , xn+1,p ), on a : 1 ,. , n , n+1 i.i.d. p prev Ŷn+1 ∼ N ( β0 + ∑ βj xn+1, j = xn+1,. β. , σ 2 xn+1,. (t x.. x.. )−1 t xn+1,. ) j=1 prev Cas de la RLS : avec xn+1 ∈ R la nvelle donnée non étiquetée, Ŷn+1 68/97 σ Ŝ √ √ prev Ŷn+1 −xn+1,. β. xn+1,. (t x.. x.. )−1 t xn+1,. prev Ŷn+1 −xn+1,. β. xn+1,. (t x.. x.. )−1 t xn+1,. ∼ N (0, 1) ∼ N ( β0 + β1 xn+1 , σ 2 ( n1 + (Cas de la RLS : ∼ t(n − (p + 1)) prev −(β0 +β1 xn+1 ) Ŷ n+1 ¿ Á 1 (xn+1 −x )2 Á À σ + n nσ̂x2 (Cas de la RLS : (xn+1 −x )2 nσ̂x2 ∼ N (0, 1) ) prev −(β0 +β1 xn+1 ) n+1 ¿ Á 1 (xn+1 −x )2 À + Ŝ Á n nσ̂x2 Ŷ )) ∼ t(n − 2) ) prev ˆ prev L’erreur de prévision est n+1 ∶= Yn+1 − Ŷn+1. On a : 2 t −1 t ˆ prev xn+1,. ) ) n+1 ∼ N ( 0 , σ ( 1 + xn+1,. ( x.. x.. ) (x −x )2 prev (Cas de la RLS : ˆ n+1 ∼ N ( 0 , σ 2 ( 1 + n1 + n+1 2 ))) nσ̂ σ √ (Cas de la RLS : Ŝ √ ˆ prev n+1 ¿ Á (x −x )2 1 Á À σ 1+ + n+1 2 n nσ̂x ˆ prev n+1 ∼ N (0, 1) ∼ N (0, 1) ) 1 + xn+1,. (t x.. x.. )−1 t xn+1,. (Cas de la RLS : 69/97 x ˆ prev n+1 1 + xn+1,. (t x.. x.. )−1 t xn+1,. ˆ prev n+1 ¿ 2 Á À 1+ 1 + (xn+1 −x ) Ŝ Á n nσ̂x2 ∼ t(n − (p + 1)) ∼ t(n − 2) ) Intervalle de prévision de Yn+1 L’intervalle de prévision de Yn+1 au niveau 1 − α est (Y ) IP1−αn+1 ((x1. ,Y1 ),. , (xn. , Yn )) t(n−(p+1)) prev ± q1− α ∶= [Ŷn+1 2 (Y Ŝ √ 1 + xn+1,. (t x.. x.. )−1 t xn+1,. ] ) ↝ P( Yn+1 ∈ IP1−αn+1 ((x1. , Y1 ),. , (xn. , Yn )) ) = 1 − α ↝ pr une réal. ((x1. , y1 ),. , (xn. , yn )) de l’obs. ((x1. , Y1 ),. , (xn. , Yn )), (Y ) (Y ) ip1−αn+1 ∶= IP1−αn+1 ((x1. , y1 ),. , (xn. , yn )) √ t(n−(p+1)) prev ± q1− α = [ŷn+1 ŝ 1 + xn+1,. (t x.. x.. )−1 t xn+1,. ] 2 70/97 (Y √ ) t(n−2) Ŝ 1− α 2 (Yn+1 ) IP1−α ((x1 , y1 ), n+1 ((x , Y ),. , (x , Y )) ∶= [Ŷ prev ± q RLS : IP1−α n n 1 1 n+1 (Y ) n+1 ∶= de l’obs. ((x1 , Y1 ),. , (xn , Yn )), ip1−α (x −x )2 ] et pr une réal. ((x1 , y1 ),. , (xn , yn )) 1 + n1 + n+1 2 nσ̂x √ (x −x )2 t(n−2) prev. , (xn , yn )) = [ŷn+1 ± q α ŝ 1 + n1 + n+1 2 ] 1− prev Application numérique pour notre exemple des appartements : pour xn+1 = 77, ŷn+1 (Y 2 nσ̂x ≃ 329.926, et pour α = 5%, ) n+1 ≃ [231.595, 428.257] ip0.95 xnplus1 = 77 alpha = 0.05 hatprevynplus1 = hatbeta0 + hatbeta1 * xnplus1 qalpha = scipy.stats.t.ppf(1 - alpha/2, df=n-2) precision = qalpha * np.sqrt(s2 * ( 1 + 1/n + ( xnplus1 - np.mean(surface) )**2 / (n * np.var(surface)))) intervalleprev = hatprevynplus1 + np.array([-1,1]) * precision ou bien, avec la fonction OLS : xnplus1 = 77 xnplus1augm = np.array([1, xnplus1]) hatprevynplus1 = reg.predict(xnplus1augm) regnouv = reg.get− prediction(xnplus1augm) intervalleprev = regnouv.conf− int (obs=True, alpha=0.05) 71/97 On obtient une "région de prévision" pour y autour de l’hyperplan affine de p régression linéaire estimé (x1 ,. , xp ) ↦ βˆ 0 + ∑ βˆ j xj : j=1 p t(n−(p+1)) (x1 ,. , xp ) ↦ [ βˆ 0 + ∑βˆ j xj ± q1− α Ŝ j=1 2 √ 1 + xn+1,. (t x.. x.. )−1 t xn+1,. ] ˆ +β ˆ x Cas de la RLS : on obtient une "région de prévision" pour y autour de la droite de régression linéaire estimée x ↦ β 0 1 ¿ 2 ˆ +β ˆ x ± q t(n−2) Ŝ Á Á À 1 + 1 + (x − x ) ] x↦[β 0 1 1− α n nσ̂x2 2 xnplus = np.linspace(10, 120, num=1000) alpha = 0.05 avec la fonction OLS : xnplusaugm = statsmodels.api.add− constant(xnplus) hatprevynplus = reg.predict(xnplusaugm) regnouv = reg.get− prediction(xnplusaugm) regionprevisionynplus = regnouv.conf− int(obs=True, alpha=0.05) borneinfprev = regionprevisionynplus[:,0] bornesupprev = regionprevisionynplus[:,1] 72/97 : ou bien : hatprevynplus = hatbeta0 + hatbeta1 * xnplus qalpha = scipy.stats.t.ppf(1 - alpha/2, df=n-2) precision = qalpha * np.sqrt(s2 * (1 + 1/n + (xnplus-np.mean(surface))**2 / (n * np.var(surface)))) borneinfprev = hatprevynplus - precision bornesupprev = hatprevynplus + precision ou encore, avec les notations matricielles : tsurfaceaugmm = np.transpose(surfaceaugmm) approxvarhatbetam = s2 * scipy.linalg.inv(np.dot(tsurfaceaugmm, surfaceaugmm)) nnouv = len(xnplus) xnplusaugmm = np.ones((nnouv, 2)) xnplusaugmm[:,1] = xnplus hatprevynplusm = np.dot(xnplusaugmm, hatbetam) txnplusaugmm = np.transpose(xnplusaugmm) approxvarprevm = np.dot(np.dot(xnplusaugmm, approxvarhatbetam), txnplusaugmm) + s2 * np.identity(nnouv) precisionm = qalpha * np.sqrt(np.diag(approxvarprevm)) borneinfprevm = hatprevynplusm - precisionm bornesupprevm = hatprevynplusm + precisionm 73/97 Intervalle de confiance de E[Yn+1] p Intervalle de confiance de E[Yn+1 ] = β0 + ∑ βj xn+1, j de niveau 1 − α j=1 (E[Yn+1 ]) IC1−α ((x1. ,Y1 ),. , (xn. , Yn )) t(n−(p+1)) prev ± q1− α ∶= [Ŷn+1 2 (E[Yn+1 ]) ↝ P( E[Yn+1 ] ∈ IC1−α Ŝ √ xn+1,. (t x.. x.. )−1 t xn+1,. ] ((x1. , Y1 ),. , (xn. , Yn )) ) = 1 − α ↝ pr une réal. ((x1. , y1 ),. , (xn. , yn )) de l’obs. ((x1. , Y1 ),. , (xn. , Yn )), (E[Yn+1 ]) ic1−α (E[Yn+1 ]) ∶= IC1−α prev = [ ŷn+1 74/97 ((x1. , y1 ),. , (xn. , yn )) √ t(n−(p+1)) ± q1− α ŝ xn+1,. (t x.. x.. )−1 t xn+1,. ] 2 RLS : Intervalle de confiance de E[Yn+1 ] = β0 + β1 xn+1 de niveau 1 − α √ 2 (E[Y ]) t(n−2) prev 1 + (xn+1 −x ) ] IC1−α n+1 ((x1 , Y1 ),. , (xn , Yn )) ∶= [Ŷn+1 ± q α Ŝ 2 n 1− nσ̂x 2 ↝ pour une réalisation ((x1 , y1 ),. , (xn , yn )) de l’observation ((x1 , Y1 ),. , (xn , Yn )), √ 2 (E[Y ]) (E[Y ]) t(n−2) prev 1 + (xn+1 −x ) ic1−α n+1 ∶= IC1−α n+1 ((x1 , y1 ),. , (xn , yn )) = [ ŷn+1 ± q α ŝ ] 2 n 1− nσ̂x 2 prev Application numérique pour notre exemple des appartements : pour xn+1 = 77, ŷn+1 (E[Y ]) ic0.95 n+1 ≃ 329.926, et pour α = 5%, ≃ [303.014, 356.838] xnplus1 = 77 alpha = 0.05 hatprevynplus1 = hatbeta0 + hatbeta1 * xnplus1 qalpha = scipy.stats.t.ppf(1 - alpha/2, df=n-2) precision = qalpha * np.sqrt(s2 * ( 1/n + ( xnplus1 - np.mean(surface) )**2 / (n * np.var(surface)))) intervalleconf = hatprevynplus1 + np.array([-1,1]) * precision ou bien, avec la fonction OLS : xnplus1 = 77 xnplus1augm = np.array([1, xnplus1]) hatprevynplus1 = reg.predict(xnplus1augm) regnouv = reg.get− prediction(xnplus1augm) intervalleconf = regnouv.conf− int (obs=False, alpha=0.05) 75/97 On obtient une "région de confiance" pour la valeur moyenne de y autour de p l’hyperplan affine de régression linéaire estimé (x ,. , x ) ↦ βˆ + ∑ βˆ x : 1 p t(n−(p+1)) (x1 ,. , xp ) ↦ [ βˆ 0 + ∑βˆ j xj ± q1− α Ŝ j=1 2 √ p 0 j=1 j j xn+1,. (t x.. x.. )−1 t xn+1,. ] Cas de la RLS : On obtient une "région de confiance" pour la valeur moyenne de y autour de la droite de régression linéaire ˆ +β ˆ x par estimée x ↦ β 0 1 ¿ 2 ˆ +β ˆ x ± q t(n−2) Ŝ Á Á À 1 + (x − x n ) ] x↦[β 0 1 1− α n nσ̂x2 2 xnplus = np.linspace(10, 120, num=1000) alpha = 0.05 avec la fonction OLS : xnplusaugm = statsmodels.api.add− constant(xnplus) hatprevynplus = reg.predict(xnplusaugm) regnouv = reg.get− prediction(xnplusaugm) regionconfianceespynplus = regnouv.conf− int(obs=False, alpha=0.05) borneinfconf = regionconfianceespynplus[:,0] bornesupconf = regionconfianceespynplus[:,1] 76/97 ou bien : hatprevynplus = hatbeta0 + hatbeta1 * xnplus qalpha = scipy.stats.t.ppf(1 - alpha/2, df=n-2) precision = qalpha * np.sqrt(s2 * ( 1/n + (xnplus-np.mean(surface))**2 / (n * np.var(surface)))) borneinfconf = hatprevynplus - precision bornesupconf = hatprevynplus + precision ou encore, avec les notations matricielles : tsurfaceaugmm = np.transpose(surfaceaugmm) approxvarhatbetam = s2 * scipy.linalg.inv(np.dot(tsurfaceaugmm, surfaceaugmm)) nnouv = len(xnplus) xnplusaugmm = np.ones((nnouv, 2)) xnplusaugmm[:,1] = xnplus hatprevynplusm = np.dot(xnplusaugmm, hatbetam) txnplusaugmm = np.transpose(xnplusaugmm) approxvarconfm = np.dot(np.dot(xnplusaugmm, approxvarhatbetam), txnplusaugmm) precisionm = qalpha * np.sqrt(np.diag(approxvarconfm)) borneinfconfm = hatprevynplusm - precisionm bornesupconfm = hatprevynplusm + precisionm 77/97 C)3. Pertinence du modèle Variabilité Somme des carrés n Expliquée ddl 2 SCE = ∑ ( Ŷi − Y. ) (degrés de liberté) p i=1 1 pr la RLS Table d’ANOVA : n Résiduelle 2 SCR = ∑ ( Yi − Ŷi ) i=1 n − (p + 1) n − 2 pr la RLS n Totale SCT = ∑ ( Yi − Y. ) 2 n−1 i=1 variabilité totale = variabilité expliquée + variabilité résiduelle ´¹¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¶ ´¹¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¶ ´¹¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¸ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¶ SCT SCE SCR Le coefficient de détermination est la v.a.r. SCE SCR R 2 = R 2 (Y1 ,. , Yn ) ∶= =1− SCT SCT Pr une réal. ((x1. , y1 ),... , (xn. , yn )) de l’obs. ((x1. , Y1 ),... , (xn. , Yn )), on a r 2 = R 2 (y1 ,... , yn ) = sce = 1 − scr sct sct 78/97 Test de Fisher issu de l’ANOVA de l’absence de liaison linéaire entre les xi et les yi L’absence de liaison linéaire entre les xi et les yi se traduit par la nullité de β1 ,. , βp. On procède donc au test bilatère : H0 ∶ β1 =... = βp = 0 contre H1 ∶ ∃j ∈ J1, pK, βj ≠ 0 Observation : ((x1. , Y1 ),. , (xn. , Yn )) ∈ (Rp × R)n Y. ∼ N (x.. β. , σ 2 In ) β0 , β1 ,. , βp ∈ R, σ ∈ 2 ⎛Y1 ⎞ ⎛1 ( rappel : Y. = ⎜ ⋮ ⎟, x.. = ⎜ ⋮ ⎝Yn ⎠ ⎝1 R⋆+ (noté (Y1 ,. , Yn )) x11 ⋮ xn1...... ⎛β0 ⎞ x1p ⎞ ⎜β ⎟ ⋮ ⎟, β. = ⎜ 1 ⎟ ) ⎜ ⋮ ⎟ xnp ⎠ ⎝βp ⎠ inconnus Statistique de Fisher : Tn = Tn (Y1 ,. , Yn ) ∶= SCE /p SCR/(n−(p+1)) sous H0 , Tn ∼ F(p, n − (p + 1)) (F étant la loi de Fisher). sous H1 , Tn aura tendance à prendre des valeurs plus grandes que sous H0 79/97 Test statistique exact, dit de Fisher, de taille α ∈]0, 1[ : φα (Y1 ,. , Yn ) = 1T (Y ,. ,Y )>qF (p,n−(p+1)) n 1 n 1−α p-valeur associée à une réalisation (y1 ,. , yn ) : π(y1 ,. , yn ) = Pβ1 =0,. ,βp =0 (Tn (Y1 ,. , Yn ) ≥ Tn (y1 ,. , yn )) = 1 − FF (p,n−(p+1)) (Tn (y1 ,. , yn )) Ex pr p = 2 (basketteurs) : on se restreint aux n = 50 premiers, on trouve T50 (y1 ,. , y50 ) ≃ 1.615 F (p,n−(p+1)) F (2,47) pr α = 5%, q1−α = q0.95 ≃ 3.195 donc F (p,n−(p+1)) T50 (y1 ,. , y50 ) ≤ q1−α et l’on conserve H0 au niveau α = 5% on trouve comme p-valeur : π(y1 ,. , y50 ) ≃ 0.21 donc pour α < 0.21 on conserve H0 au niveau α pour α > 0.21 on rejette H0 au niveau α et on accepte H1 stat− Fisher = reg.fvalue seuil− alpha = stats.f.ppf(1 - alpha, dfn=p, dfd = n-(p+1)) pval− Fisher = reg.f− pvalue 80/97 Cas de la RLS : Test statistique de Fisher issu de l’ANOVA de significativité du modèle, autrement dit de l’absence de liaison linéaire entre les xi et les yi , c-à-d test bilatère de H0 ∶ β1 = 0 contre H1 ∶ β1 ≠ 0 Observation : ((x1 , Y1 ),. , (xn , Yn )) ∈ (R × R)n (noté (Y1 ,. , Yn )) Y. ∼ N (x.. β. , σ 2 In ) ⎛ Y1 ⎞ ⎛1 ( rappel : Y. = ⎜ ⋮ ⎟, x.. = ⎜ ⋮ ⎝Yn ⎠ ⎝1 β0 , β1 ∈ R, σ 2 ∈ R⋆ + inconnus x1 ⎞ β ⋮ ⎟, β. = ( 0 ) ) β1 xn ⎠ SCE Statistique utilisée : T (Y1 ,. , Yn ) ∶= SCR n−2 sous H0 , T (Y1 ,. , Yn ) ∼ F (1, n − 2) (loi de Fisher) βˆ 1 2 SCE = ( √ Rq : SCR ) et, par définition de la loi de Ŝ/ nσ̂ 2 n−2 x Fisher, le carré d’une v.a.r. de loi t(n − 2) est de loi F (1, n − 2) sous H1 , T (Y1 ,. , Yn ) aura tendance à prendre des valeurs plus grandes que sous H0 Test exact, dit de Fisher, de taille α ∈]0, 1[ : φα (Y1 ,. , Yn ) = 1 F (1,n−2) T (Y1 ,. ,Yn )>q 1−α Application numérique dans notre exemple : T (y1 ,. , y20 ) = sce scr ≃ 96.259, et pour α = 5%, F (1,n−2) q1−α F (1,18) = q0.95 n−2 ≃ 4.414 donc on retrouve le fait que l’on rejette H0 ∶ β1 = 0 au niveau 5% p-valeur associée à une réalisation (y1 ,. , yn ) : π(y1 ,. , yn ) = Pβ1 =0 (T (Y1 ,. , Yn ) ≥ T (y1 ,. , yn )) = 1 − FF (1,n−2) (Tn (y1 ,. , yn )) Application numérique dans notre exemple : π(y1 ,. , y20 ) ≃ 1 − FF (1,18) (96.259) ≃ 1.197e − 08 statFisher = sce/s2 quantilealpha = scipy.stats.f.ppf(1 - alpha, dfn=1, dfd = n-2) phialpha = (statFisher > quantilealpha) pvalFisher = 1 - scipy.stats.f.cdf(statFisher, dfn=1, dfd = n-2) Remarque : ce test de Fisher est ÉGAL au test de Student de nullité de β1 (attention, c’est spécifique à la RLS) 81/97 D) Limites du modèle & Sélection de variables Nous allons enfin, à la lumière des résultats obtenus, procéder à un retour critique sur la qualité de notre modèle et sur les conditions que nous avions imposées initialement. Objectifs : Etude de la qualité du modèle Retour critique sur les conditions (A0)-(A1)-(A2)-(A3)-(A4) Voir si nous ne pouvons pas supprimer certaines des p variables explicatives pour "alléger" le modèle Essayer de déterminer si notre modèle est suffisamment "robuste" par rapport aux données, en regardant quels sont les points leviers et s’il y a des outliers. 82/97 D)1. Adéquation du modèle On peut évaluer graphiquement l’adéquation du modèle en traçant le graphique (yi , ŷi )i∈J1,nK puisque les ŷi devaient "approximer" les yi : les points doivent être alignés le long de la première bissectrice, cela signifie que l’adéquation du modèle aux données est correcte ↝ si ce n’est pas le cas il y a une mauvaise adéquation du modèle aux données et il faut changer de modèle Exemple pour p = 1 (appartements) : ↝ les points sont bien proches de la 1ère bissectrice 83/97 Exemple pour p = 2 (basketteurs) : ↝ les points sont plutôt proches de la 1ère bissectrice, mais pas les points extrêmes D)2. Retour critique sur les conditions (A0-A4) Condition (A0) : n > p + 1 ↝ vérification immédiate t x.. x.. inversible ↝ on vérifie que son déterminant est non nul : linalg.det(np.dot(np.transpose(x− tab), x− tab)) != 0 Conditions (A1)-(A2)-(A3) : 1 ,., n devaient être aléatoires, centrés, de même variance et non-corrélés entre eux : nous n’avons pas ˆi accès aux i , mais nous avons accès aux ↝ les résidus doivent avoir un "comportement aléatoire, sans structure évidente" (pente, entonnoir,...) on trace le graphe des (i, ϵ̂i )i∈J1,nK ↝ si ce graphe présente une structure, alors une des conditions (A1)-(A2)-(A3) peut ne pas être vérifiée : il peut par exemple y avoir une autocorrélation (t ext) on trace le graphes des (ŷi , ϵ̂i )i∈J1,nK et des (ŷi , ϵ̂i )i∈J1,nK ↝ une structure particulière dans ces graphes peut indiquer que l’homoscédasticité n’est pas vérifiée 84/97 S’il y a une structure évidente, il faut changer de modèle pour essayer de prendre en compte cette structure (par exemple rajouter un terme quadratique x2 ou polynomial dans les variables explicatives,...) Exemple pour p = 1 Exemple pour p = 2 ↝ la variance semble augmenter avec i (les joueurs sont classés par taille croissante) ↝ il semblerait que les ϵ̂i aient tendance à décroître... ↝ alors qu’ici il n’y a pas de structure particulière à signaler ↝ structure en entonnoir, la variance augmente avec les ŷi (elle dépend donc des xi ) : il y a hétéroscédasticité et il faut changer de modèle pour essayer de la prendre en compte 85/97 Notations supplémentaires Pour i ∈ J1, nK, on note hii = ( ΠIm(x.. ) )ii = ( x.. (t x.. x.. )−1 t x.. )ii la i ème valeur sur la diagonale de la matrice de projection orthogonale sur Im(x.. ) (= hat matrix ) ↝ hii représente le poids de l’obs. i sur sa propre estimation ˆi ˆ (1) ˆ i ] = σ 2 (1 − hii )) (rappel : Var[ i ème résidu normalisé : ∶= σ√1−h i ii (attention, ils ne sont pas observables à cause de σ inconnu) ˆ (t) i ème résidu standardisé : i ∶= ˆ √ i Ŝ 1−hii (t) ( Cas de la RLS : ˆ i ∶= ˆ i ¿ 2 Á À1− 1 − (xi −x ) Ŝ Á n nσ̂x2 ) ˆ1 ,... , ˆ n ne sont pas indépendants et ne peuvent donc Attention : pas être représentatifs d’une absence/présence de structuration par autocorrélation. √ (t ext) (t) n−(p+2) ˆ ˆ i ème résidu studentisé (externe) : i ∶= i 2 ˆ (t) n−(p+1)−( i ) ↝ reg.outlier.test()[:,0] (t) (t) Tous ces résidus sont de même variance 1. 86/97 Condition (A4) : 1 ,... , n ∼ N (0, σ 2 ) pour un certain σ 2 c-à-d que nous avons supposé que les erreurs sont aléatoires, indépendantes, centrées, gaussiennes, de même variance. Nous n’avons pas accès aux i , mais nous avons accès aux réalisations ext) ˆ i et aux résidus studentisés ˆ (t des résidus . i Or, si (A4) est vérifiée, et si retirer la i ème ligne de la matrice x.. ne ext) ˆ (t change pas son rang, ∼ t(n − (p + 2)) i Nous allons donc vérifier sur ces résidus que la condition (A4) sur les ˆ i n’est pas aberrante i.i.d. Nous utilisons pour cela : des méthodes graphiques des tests 87/97 D)3. Sélection de variables On dispose d’un modèle de RLM à p variables explicatives. On souhaite simplifier notre modèle et ne garder que certaines variables explicatives. Comment les choisir au mieux ? Nous avons besoin : d’un critère de qualité d’un modèle, permettant de comparer deux modèles n’ayant pas nécessairement le même nombre de variables explicatives d’une procédure de choix de modèles, qui permet de choisir, parmi tous les modèles, le meilleur au sens de ce critère. Problème de complexité : le nombre de modèles à considérer est p q ∑ Cp = 2p − 1 , qui croît exponentiellement avec p. Par exemple, si p = 30, q=1 on devrait considérer 230 = 109 modèles... En pratique : on utilise des heuristiques dont les plus simples sont les procédures pas à pas ascendante/descendante 88/97 Critères de qualité d’un modèle Le premier critère de qualité qui mesure l’ajustement du modèle aux données que nous avons vu est le coefficient de détermination R2 = 1 − σ̃ SCR =1−. SCT σ̃Y. Problème : le R 2 augmente lorsque le nombre de variables explicatives incluses dans le modèle augmente et il ne permet donc pas de comparer deux modèles n’ayant pas le même nombre de variables explicatives, seulement deux modèles ayant le même nombre de variables explicatives. Pour pallier cette difficulté, on introduit le coefficient de détermination ajusté n−1 SCR/(n − p − 1) n−p−1 σ̃. 2 Rajusté =1− =1− SCT /(n − 1) σ̃Y. 2 ↝ Rajusté n’augmente pas forcément lorsque le nombre de variables incluses dans le modèle augmente, et permet de comparer des modèles ayant un nombre de variables différent 89/97 Si notre modèle de RLM dispose d’une structure permettant de lui associer une vraisemblance, ce qui est le cas dans le cadre du C) Modèle Linéaire Gaussien, on dispose également d’un critère de vraisemblance pénalisé, le Critère d’Information de Akaike (Akaike Information Criterion) : AIC ∶= −2ℓ⋆n + 2k où ℓ⋆n est la log-vraisemblance maximisée k est le nombre de paramètres libres du modèle : ds le modèle de RLM à q variables explicatives, il y a q + 2 paramètres, β0 , β1 ,... , βq , σ 2 , et une équation, y = β0 + β1 x1 +... + βp xp , donc k = (q + 2) − 1 = q + 1 paramètres libres la fct de vraisemblance associée à une réal. ((x1. , y1 ),... , (xn. , yn )) de l’obs. ((x1. , Y1 ),... , (xn. , Yn )) est la densité jointe de ((x1. , Y1 ),... , (xn. , Yn )) au point ((x1. , y1 ),... , (xn. , yn )) Ln (. ; ((x1. , y1 ),... , (xn. , yn )) ) ∶ (β. , σ 2 ) ∈ Rp+1 × R⋆+ ↦ f((x 2 1. ,Y1 ),...,(xn. ,Yn ));(β. ,σ ) 90/97 ((x1. , y1 ),... , (xn. , yn )) = e − 1 2 ∣∣y. −x..β. ∣∣2 n R 2σ (2πσ 2 )n/2 l’estimateur du maximum de vraisemblance de (β. , σ 2 ) est (β̂.MV , σ̂ 2 MV ) = ( (t x.. x.. )−1 t x.. Y. , SCR ) n Preuve (même méthode que Ex 3 Feuille 1 de TD) : on commence par optimiser en β. : pour σ02 ∈ R⋆ + fixé, β. ∈ Rp+1 ↦ ℓn ( β. , σ02 ; ((x1. , y1 ),... , (xn. , yn )) ) = − n2 ln(2π) − n2 ln(σ02 ) − 1 2 ∣∣y. − x..β. ∣∣2 n est deux fois dérivable et concave, de plus ∇β. ℓn ( β. , σ02 ; ((x1. , y1 ),... , (xn. , yn )) ) = 2σ 0 1 t x (y... 2 σ 0 R − x.. β. ), donc ℓn (. , σ02 ; ((x1. , y1 ),... , (xn. , yn )) ) atteint son maximum en son unique point critique β̂. = (t x.. x.. )−1 t x.. y. 2 On optimise maintenant la fonction σ 2 ∈ R⋆ + ↦ ℓn ( β̂. , σ ; ((x1. , y1 ),... , (xn. , yn )) ) = − n2 ln(2π) − n2 ln(σ 2 ) − 1 2 ∣∣y. − x..β̂. ∣∣2 n = − n2 ln(2π) − n2 ln(σ 2 ) − 1 2 scr 2σ R 2σ 1 ( − nσ 2 + scr ) qui est ≥ 0 ssi σ 2 ≤ scr On a ∂σ2 ℓn ( β̂. , σ 2 ; ((x1. , y1 ),... , (xn. , yn )) ) = − n 2 + scr2 2 = n 2σ 2(σ ) 2(σ 2 )2 donc ℓn ( β̂. ,. ; ((x1. , y1 ),... , (xn. , yn )) ) atteint son maximum en σ̂ 2 = scr n le maximum de la log-vraisemblance est : n ℓ⋆n = ℓn (β̂.MV , σ̂ 2 MV ; ((x1. , y1 ),... , (xn. , yn )) ) = − ( ln(2πσ̂ 2 MV ) + 1) 2 On obtient donc AIC = n ln(SCR) + 2k + cste 91/97 2 Ces critères, du coefficient de détermination ajusté Rajusté ou d’information de Akaike AIC, doivent maintenant être optimisés dans une procédure de choix de modèle / sélection de variables. Nous présentons les deux plus simples, la procédure pas à pas ascendante et la procédure pas à pas descendante. 92/97 Procédure pas à pas ascendante Procédure pas à pas ascendante (= forward stepwise) : On choisit un critère On part du modèle nul sans variable On effectue les p RLS et on sélectionne le modèle qui optimise le 2 critère choisi c-à-d qui maximise le Rajusté (resp. minimise le AIC) si tel est le critère choisi ↝ on sélectionne ainsi une première variable Vi1 On effectue les p − 1 RLM avec Vi1 et une autre variable explicative choisie parmi les p − 1 restantes et on sélectionne le modèle qui optimise le critère choisi ↝ on sélectionne ainsi une deuxième variable Vi2 On recommence jusqu’à ce qu’on ne puisse plus optimiser le critère 2 c-à-d que le Rajusté ne puisse plus augmenter (resp. le AIC ne puisse plus diminuer) si tel est le critère choisi 93/97 Procédure pas à pas descendante Procédure pas à pas descendante (= backward stepwise) : On choisit un critère On part du modèle complet à p variables On effectue les p RLM à p − 1 variables et on sélectionne le modèle qui 2 optimise le critère choisi c-à-d qui maximise le Rajusté (resp. minimise le AIC) si tel est le critère choisi ↝ on supprime ainsi une première variable Vj1 On effectue les p − 1 RLM à p − 2 variables choisies parmi celles différentes de Vj1 et on sélectionne le modèle qui optimise le critère choisi ↝ on supprime ainsi une deuxième variable Vj2 On recommence jusqu’à ce qu’on ne puisse plus optimiser le critère 2 c-à-d que le Rajusté ne puisse plus augmenter (resp. le AIC ne puisse plus diminuer) si tel est le critère choisi 94/97 D)4. Robustesse Pour i ∈ J1, nK, on dit que le point (xi. , yi ) est un point influent s’il a un poids disproportionné, qu’il influe beaucoup sur les paramètres estimés (quand on l’enlève, la droite de régression estimée change beaucoup) (t) hii ϵ̂ élevée (> 1) ↝ il est caractérisé par une distance de Cook 21 1−h ii i Un tel point influent (xi. , yi ) est soit un point levier ↝ en général, point dans la continuité des autres mais "loin" des autres s’il vérifie un des critères suivants : critères de Welsch : hii > 2p n ou (hii > critère de Huber : hii > 0.5 3p n soit un outlier = point aberrant/atypique pr p > 6 et n − p > 12) ↝ en général, point qui n’est pas "dans la continuité des autres", et est éloigné des régions de confiance/prédiction (t ext) si ∣ϵ̂i t(n−(p+2)) ∣ > q1− 1 n ↝ outlier.test statsmodels.stats.outliersinfluence.OLSInfluence(reg) soit les deux à la fois ! 95/97 Annexe : Commande OLS de statsmodels surfaceaugm = statsmodels.api.add− constant(surface) ↝ rajout de la colonne de 1 modelelin = statsmodels.api.OLS(prix,surfaceaugm) reg = modelelin.fit() ↝ initialise le modèle print(reg.summary()) ↝ OLS Regression Results ============================================================================== Dep. Variable: y R-squared: 0.842 Model: OLS Adj. R-squared: 0.834 Method: Least Squares F-statistic: 96.26 Date: Fri, 10 Mar 2023 Prob (F-statistic): 1.20e-08 Time: 18:28:07 Log-Likelihood: -103.47 No. Observations: 20 AIC: 210.9 Df Residuals: 18 BIC: 212.9 Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] const 33.6438 24.445 1.376 0.186 -17.713 85.001 x1 3.8478 0.392 9.811 0.000 3.024 4.672 Omnibus: 0.506 Durbin-Watson: 0.990 Prob(Omnibus): 0.777 Jarque-Bera (JB): 0.583 Skew: -0.125 Prob(JB): 0.747 Kurtosis: 2.202 Cond. No. 151 96/97 reg.params ↝ contient le vecteur des paramètres estimés β̂. reg.conf− int(alpha)[j] ↝ contient l’intervalle de confiance de niveau 1 − α de βj reg.tvalues[j] ↝ contient la valeur de la statistique du test de Student de significativité de βj reg.pvalues[j] ↝ contient la p-valeur du test de Student de significativité de βj reg.fittedvalues ↝ contient le vecteur des prédictions ŷ. (la valeur de ŷi est à l’indice i − 1) reg.resid ↝ contient le vecteur des résidus ϵ̂. (la valeur de ϵ̂i est à l’indice i − 1) reg.ssr ↝ contient la valeur de la somme des carrés des résidus scr reg.ess ↝ contient la valeur de la somme des carrés expliquée sce reg.centered− tss ↝ contient la valeur de la somme des carrés totale sct reg.rsquared ↝ contient la valeur du coefficient de détermination r 2 reg.get− prediction(x− new).conf− int(obs=True, alpha=0.05) ↝ contient l’intervalle de prévision de Yn+1 pour x− new= (1, xn+1,1 ,. xn+1,p ) une nvelle donnée non étiquetée reg.get− prediction(x− new).conf− int(obs=False,alpha=0.05) ↝ contient l’intervalle de confiance de E[Yn+1 ] pour x− new= (1, xn+1,1 ,. xn+1,p ) une nvelle donnée non étiquetée reg.fvalue ↝ contient la valeur de la statistique du test de Fisher de significativité du modèle reg.f− pvalue ↝ contient la p-valeur du test de Fisher de significativité du modèle 97/97

Cours Chapitre 3 Régression Linéaire PDF A3S6 2023-2024

Document Details

Tags

Related

Summary

Full Transcript

Upgrade to continue