Modèle linéaire gaussien

Didier Chauveau et Jean-François Delmas
(last modiﬁcation date: October 10, 2017)

Version pdf de ce document
Version sans bandeaux

1 Peut-on augmenter la quantité de pluie?
2 (Facultatif) Analyse de la variance à un facteur
3 Régression linéaire

1 Peut-on augmenter la quantité de pluie?
1.1 La problématique
1.2 Les données
1.3 Le modèle
1.4 Comparaison des moyennes
1.5 (Facultatif) Comparaison des variances des 2 échantillons gaussiens
2 (Facultatif) Analyse de la variance à un facteur
2.1 Les données
2.2 Le modèle
3 Régression linéaire
3.1 Les données
3.2 Régression linéaire simple
3.3 Régression linéaire multiple

1 Peut-on augmenter la quantité de pluie?

1.1 La problématique

Il existe deux types de nuages qui donnent lieu à des précipitations: les nuages chauds et les nuages froids. Ces derniers possèdent une température maximale de l’ordre de -10°C à -25°C. Ils sont composés de cristaux de glace et de gouttelettes d’eau. Ces gouttelettes d’eau subsistent alors que la température ambiante est inférieure à la température de fusion. On parle d’eau surfondue. Leur état est instable. De fait, quand une particule de glace rencontre une gouttelette d’eau, elles s’aggrègent pour ne former qu’une seule particule de glace. Les particules de glace, plus lourdes que les gouttelettes, tombent sous l’action de la gravité. Enﬁn si les températures des couches d’air inférieures sont suﬃsamment élevées, les particules de glace fondent au cours de leur chute formant ainsi de la pluie.

En l’absence d’un nombre suﬃsant de cristaux de glace pour initier le phénomène décrit ci-dessus, on peut injecter dans le nuage froid des particules qui ont une structure cristalline proche de la glace, par exemple de l’iodure d’argent (environ 100 à 1000 grammes par nuage). Autour de ces particules, on observe la formation de cristaux de glace, ce qui permet, on l’espère, de déclencher ou d’augmenter les précipitations. Il s’agit de l’ensemencement des nuages. Signalons que cette méthode est également utilisée pour limiter le risque de grêle.

Il est évident que la possibilité de modiﬁer ainsi les précipitations présente un grand intérêt pour l’agriculture. De nombreuses études ont été et sont encore consacrées à l’étude de l’eﬃcacité de l’ensemencement des nuages dans divers pays. L’étude de cette eﬃcacité est cruciale et délicate. Le débat est encore d’actualité.

1.2 Les données

On dispose des données concernant l’ensemencement par iodure d’argent des nuages en Floride (1975)¹ .

Les volumes de pluie déversées en 10⁷m³ (cf. les deux tableaux ci-dessous) concernent 23 jours similaires dont n₁ = 11 jours avec ensemencement correspondant aux réalisations des variables aléatoires X₁,…,X_n₁ et n₂ = 12 jours sans ensemencement correspondant aux réalisations des variables aléatoires Y ₁,…,Y _n₂.


i	1	2	3	4	5	6	7	8	9	10	11

X_i	7.45	4.70	7.30	4.05	4.46	9.70	15.10	8.51	8.13	2.20	2.16

Table 1: Volume de pluie en 10⁷m³ déversée avec ensemencement


j	1	2	3	4	5	6	7	8	9	10	11	12

Y _j	15.53	10.39	4.50	3.44	5.70	8.24	6.73	6.21	7.58	4.17	1.09	3.50

Table 2: Volume de pluie en 10⁷m³ déversée sans ensemencement

X=[7.45, 4.70, 7.30, 4.05, 4.46, 9.70, 15.10, 8.51, 8.13, 2.20, 2.16];
Y=[15.53, 10.39, 4.50, 3.44, 5.70, 8.24, 6.73, 6.21, 7.58, 4.17, 1.09, 3.50];

1.3 Le modèle

On suppose que les volumes de pluies suivent une loi gaussienne et que l’eﬀet de l’ensemencement ne modiﬁe pas la variance des lois gaussiennes.

Ceci revient à dire dans ce contexte paramétrique que les observations sont constituées de deux échantillons gaussiens indépendants

X ,...,X i.i.d. 𝒩 (μ ,σ2), 1 n1 1 2 Y1,...,Yn2 i.i.d. 𝒩 (μ2, σ ).

On souhaite tester l’hypothèse nulle H₀: “le procédé n’a pas d’eﬀet” contre “le procédé produit une augmentation signiﬁcative de la quantité de pluie”, autrement dit

H = {μ = μ } contre H = {μ > μ }. 0 1 2 1 1 2

Il est indispensable de s’assurer d’abord que les hypothèses de modèle que l’on a faites sont raisonnables:

Les données peuvent-elles être considérées comme des réalisations de lois gaussienne pour les paramètres (moyenne, variance) appropriés?
Peut-on considérer que les variances sont les mêmes?

Des éléments de réponses au point (1) seront abordés ultérieurement dans le cours. Le point (2) est abordé dans le paragraphe facultatif 1.5.

Nous admettons donc dans un premier temps ces deux hypothèses.

1.4 Comparaison des moyennes

Aﬁn de se faire une idée des données et de l’écart entre les deux populations, on peut commencer par calculer les moyennes empiriques

n n ¯ -1-∑ 1 ¯ -1-∑ 2 X = n Xi, Y = n Yi, 1 i=1 2 i=1

ainsi que les variances empiriques

[ ] 1 ∑n1 n1 1 n∑1 ( )2 S2X = ------- (Xi − X¯)2 = ------- --- X2i − X¯ , n1 − 1 i=1 n1 − 1 n1 i=1 ∑n2 [ ∑n2 ( ) ] S2 = ---1--- (Yi − ¯Y )2 = --n2--- -1- Y 2 − ¯Y 2 . Y n2 − 1 i=1 n2 − 1 n2 i=1 i

// taille des echantillons
n1 = length(X);
n2 = length(Y);
// calcul des moyennes
Xbar=sum(X)/n1
Ybar=sum(Y)/n2
// calcul de la somme des carres
SSX=sum(X.^2);
SSY=sum(Y.^2);
// calcul des variances empiriques sans biais
SX2 = SSX/(n1-1)- n1/(n1-1)*Xbar^2
SY2 = SSY/(n2-1)- n2/(n2-1)*Ybar^2

On rappelle que les variables aléatoires X (resp. Y ) et (n₁ − 1)S_X²∕σ² (resp. (n₂ − 1)S_Y²∕σ²) sont indépendantes et de loi 𝒩(μ₁,σ²∕n₁) (resp. 𝒩(μ₂,σ²∕n₂)) et χ²(n₁ − 1) (resp. χ²(n₂ − 1)). Comme S_X² et S_Y² sont indépendants, on en déduit que la loi de (n − 1)S2
---1--2---X
σ + (n − 1)S2
--2---2---Y
σ est la loi du χ² à n₁ + n₂ − 2 degrés de liberté.

Sous H₀, la loi de ∘ --------
-n1n2---
n1 + n2 ( ¯X-−-¯Y)-
σ est la loi gaussienne centrée réduite. On en déduit que sous H₀,

∘ -------- -n1n2--------( ¯X-−-Y¯)----- T = n1 + n2 ∘ (n1−-1)S2+-(n2−1)S2- -----nX1+n2−-2---Y

suit une loi de Student de paramètre n₁ + n₂ − 2.

Comme sous H₁, X −Y converge p.s. vers μ₁ − μ₂ > 0 quand min(n₁,n₂) →∞, et que lim _{min(n₁,n₂)→∞} ∘ --------
-n1n2---
n1 + n2 = +∞, on en déduit que sous H₁, la statistique T diverge vers +∞ quand min(n₁,n₂) →∞.

Question 1 Construire, à partir de la statistique T, un test d’égalité des moyennes en déterminant sa région critique, la valeur critique au niveau α et la p-valeur.

Question 2 En utilisant la fonction ttest2 déﬁnie ci-dessous, calculer numériquement la valeur de la statistique de test, la valeur critique au niveau 5%, et la p-valeur pour les données de pluie. Conclusion.

La fonction ttest2 (ttest2(n1,Xbar,SX2,n2,Ybar,SY2,alpha)) retourne la valeur de la statistique T observée, la valeur critique associée au niveau α et la p-valeur du test suggéré dans la question précédente. Cliquer sur le lien ci-dessus pour obtenir le code de la fonction, le sauvegarder, sous le nom ttest2.sce, dans le répertoire où vous utilisez scilab. Pour charger la fonction, utiliser la commande: getf 'ttest2.sce'. Pour appliquer la fonction, utiliser la commande:

ttest2(n1,Xbar,SX2,n2,Ybar,SY2,alpha)

On donne quelques indications pour la compréhension du code Scilab de ttest2 qui utilise la fonction de répartion (“cumulative distribution function” en anglais) de la loi de Student:

c = cdft("T",k,p,1-p) donne à c la valeur du quantile d’ordre p de la loi de Student de paramètre k: p = ℙ(T ≤ c), où T suit la loi de Student de paramètre k.
p = cdft("PQ",t,k) donne à p la valeur p = ℙ(T ≤ t), où T suit la loi de Student de paramètre k.

Pour avoir plus d’information sur la fonction cdft, on peut consulter le manuel en utilisant la commande help cdft.

Le but de ce qui suit est d’observer les résultats du tests quand les valeurs observées ou la taille de l’échantillon varient.

Question 3

Faire varier Y , les autres valeurs étant ﬁxées. Pour quelles valeurs de Y rejetez vous H₀?
Avec la valeur de Y ainsi déterminée, qu’observez vous si S_Y² augmente, si S_Y² diminue?
En reprenant les valeurs numériques du problème, qu’observez vous si n₁ et n₂ augmentent, si n₁ et n₂ diminuent jusqu’aux valeurs limites?

Question 4 Reprendre les questions 1 et 2, si l’hypothèse nulle est {μ₁ ≤ μ₂}.

1.5 (Facultatif) Comparaison des variances des 2 échantillons gaussiens

Pour s’assurer du point (2), on peut tout d’abord supposer que ces variances sont diﬀérentes, i.e.

2 Xi ∼ 𝒩 (μ1,σ1), i = 1,...,n1, Yj ∼ 𝒩 (μ2,σ2), j = 1,...,n2, 2

puis construire un test d’égalité des variances, aﬁn de voir si la diﬀérence sur les variances (ou les écart-types) observée est signiﬁcative.

Question 5 (Facultatif) Quelle est la loi de la statistique F = S_X²∕S_Y² sous l’hypothèse nulle d’égalité des variances? Quel est son comportement sous H₁, quand n₁ → ∞ et n₂ →∞?

Question 6 (Facultatif) Construire, à partir de la statistique F, un test d’égalité des variances. Déterminer la région critique du test, évaluer la valeur critique au niveau α = 5% et la p-valeur pour les données de pluie? Conclusion.

On donne les commandes Scilab suivantes concernant la fonction de répartion de la loi de Fisher:

c = cdff("F",k1,k2,p,1-p) donne à c la valeur du quantile d’ordre p de la loi de de Fisher de paramètres (k1,k2): p = ℙ(F ≤ c), où F suit la loi de de Fisher de paramètres (k1,k2).
[p,q] = cdff("PQ",f,k1,k2) donne à p et q les valeurs p = 1 −q = ℙ(F ≤ f), où F suit la loi de Fisher de paramètres (k1,k2).

Pour avoir plus d’information sur la fonction cdff, on peut consulter le manuel en utilisant la commande help cdff.

2 (Facultatif) Analyse de la variance à un facteur

2.1 Les données

On dispose de données (réelles) relatives à 6 jeux de notes pour un même contrôle, corrigé par un unique correcteur, mais concernant six petites classes (et donc six enseignants) diﬀérentes. On désire savoir si les résultats des petites classes sont signiﬁcativement diﬀérents.

|------|-------|------|-------|------|-------| |Note--|Classe-|Note--|Classe-|-Note-|-Classe-| |15.50 | 1 |15.50 | 1 | 12.00 | 5 | |16.00 | 3 |11.75 | 2 | 18.50 | 2 | |14.00 | 5 |18.50 | 4 | 8.00 | 6 | |13.75 | 2 | 8.00 | 5 | 10.50 | 3 | | | | | | | | |17.00 | 2 |15.00 | 4 | 7.00 | 4 | |11.00 | 6 |12.50 | 2 | 11.00 | 6 | |16.50 | 1 |15.50 | 4 | 12.50 | 5 | |12.75 | 5 |19.00 | 1 | 15.00 | 5 | |15.00 | 6 |19.00 | 1 | 12.50 | 4 | |13.50 | 5 | 9.50 | 2 | 16.50 | 2 | |14.00 | 2 | 8.00 | 5 | 16.50 | 1 | | | | | | | | |19.50 | 1 |16.50 | 1 | 18.00 | 4 | |15.50 | 5 | 9.50 | 6 | 9.50 | 5 | |14.50 | 6 |19.00 | 1 | 13.50 | 6 | |13.00 | 2 |17.00 | 6 | 13.00 | 3 | |18.00 | 4 |17.00 | 5 | 16.00 | 2 | |15.00 | 5 | 9.50 | 5 | 16.00 | 4 | |14.50 | 6 |15.50 | 2 | 12.00 | 3 | |14.50 | 5 |15.00 | 4 | 18.50 | 4 | | | | | | | | |18.50 | 3 |17.00 | 2 | 12.00 | 3 | |17.50 | 1 |19.50 | 3 | 13.00 | 6 | |18.00 | 3 |16.50 | 6 | 14.50 | 3 | |11.00 | 3 | 7.00 | 4 | 12.00 | 3 | |10.50 | 5 | 8.50 | 2 | 19.50 | 1 | |10.00 | 1 |12.50 | 4 | 12.00 | 4 | |16.00 | 3 |18.75 | 6 | 10.50 | 1 | | | | | | | | |12.25 | 4 |16.00 | 4 | 15.00 | 6 | |9.50 | 1 |19.50 | 3 | 15.50 | 2 | |12.00 | 5 |17.50 | 1 | 14.25 | 3 | |13.50 | 3 |14.50 | 4 | 14.00 | 2 | |19.50 | 1 | 9.00 | 6 | 15.50 | 3 | |18.50 | 6 | | | | | |------|-------------------------------------- | |

2.2 Le modèle

On fait l’hypothèse que les notes sont distribuées suivant une loi gaussienne de même variance, et de moyenne dépendant de l’enseignant. Le modèle s’écrit donc

Xi,j = mi + 𝜀i,j, i = 1,...,6, j = 1,...,ni,

où m_i est la note moyenne (inconnue) liée à l’enseignant i, et les 𝜀_i,j sont des v.a. indépendantes de loi 𝒩(0,σ²), où σ² > 0 est inconnu également.

// notes
X=[
15.50; 15.50; 12.00; 16.00; 11.75; 18.50; 14.00; 18.50;  8.00; 13.75;
8.00; 10.50; 17.00; 15.00;  7.00; 11.00; 12.50; 11.00; 16.50; 15.50;
12.50; 12.75; 19.00; 15.00; 15.00; 19.00; 12.50; 13.50;  9.50; 16.50;
14.00;  8.00; 16.50; 19.50; 16.50; 18.00; 15.50;  9.50;  9.50; 14.50;
19.00; 13.50; 13.00; 17.00; 13.00; 18.00; 17.00; 16.00; 15.00;  9.50;
16.00; 14.50; 15.50; 12.00; 14.50; 15.00; 18.50; 18.50; 17.00; 12.00;
17.50; 19.50; 13.00; 18.00; 16.50; 14.50; 11.00;  7.00; 12.00; 10.50;
8.50; 19.50; 10.00; 12.50; 12.00; 16.00; 18.75; 10.50; 12.25; 16.00;
15.00;  9.50; 19.50; 15.50; 12.00; 17.50; 14.25; 13.50; 14.50; 14.00;
19.50;  9.00; 15.50; 18.50];
// groupe
Grp=[
1; 1; 5; 3; 2; 2; 5; 4; 6; 2; 5; 3; 2; 4; 4; 6; 2; 6; 1; 4; 5; 5; 1; 5;
6; 1; 4; 5; 2; 2; 2; 5; 1; 1; 1; 4; 5; 6; 5; 6; 1; 6; 2; 6; 3; 4; 5; 2;
5; 5; 4; 6; 2; 3; 5; 4; 4; 3; 2; 3; 1; 3; 6; 3; 6; 3; 3; 4; 3; 5; 2; 1;
1; 4; 4; 3; 6; 1; 4; 4; 6; 1; 3; 2; 5; 1; 3; 3; 4; 2; 1; 6; 3; 6];

On peut résumer les données.

// nombre de groupes
k = max(Grp);
// calcul de l'effectif n(i), de la moyenne moy(i)
// et de la variance s2(i) du groupe i
for i=1:k,
n(i) = length(X(Grp==i));
moy(i) = sum(X(Grp==i))/n(i);
ss(i)=sum(X(Grp==i).^2);
s2(i)=ss(i)/(n(i)-1)- n(i)/(n(i)-1)*moy(i)^2;
end;
[n,moy,s2] // effectif, moyenne et variance par groupe

Question 7 Quelles sont les estimations des m_i, i = 1,…, 6? Quelle est l’estimation de σ²?

Question 8 Peut-on considérer qu’il existe une diﬀérence signiﬁcative entre les notes des classes, dues aux qualités pédagogiques des enseignants? Eﬀectuer le test de l’hypothèse nulle H₀ : “les enseignants ne présentent pas de diﬀérence” contre “c’est faux”. Écrire la table d’analyse de la variance et conclure.

On pourra utiliser les commandes suivantes:

// calcul de la moyenne generale
mg=sum(X)/length(X);
// calcul des sommes de carres
SSM = sum(n.*(moy-mg).^2);
SSE = sum(s2.*(n-1));

La fonction anova (anova(n,moy,s2,alpha)) retourne la table d’analyse de la variance associée. Les paramètres sont n, vecteur des tailles de chacun des groupes; moy, vecteur des moyennes de chacun des groupes; s2, vecteur des variances de chacun des groupes; alpha, niveau du test d’analyse de la variance. Cliquer sur le lien ci-dessus pour obtenir le code de la fonction, le sauvegarder sous le nom anova.sce dans le répertoire où vous utilisez scilab. Pour charger la fonction, utiliser la commande: getf 'anova.sce'. Pour appliquer la fonction, utiliser la commande:

anova(n,moy,s2, alpha)

Le but de ce qui suit est d’observer les résultats du tests quand les valeurs observées ou la taille de l’échantillon varient.

Question 9

Faire varier le vecteur moy, les autres valeurs étant ﬁxées. Pour quelles valeurs de moy rejetez vous H₀?
Qu’observez vous si le vecteur, s2, des variances augmente, s’il diminue?
En reprenant les valeurs numériques du problème, qu’observez vous si le vecteur, n, des tailles des groupes augmente ou diminue?

Question 10 Si l’on rejette H₀, il est naturel de déterminer les classes qui sont signiﬁcativement diﬀérentes.

Proposer des intervalles de conﬁance pour les notes moyennes pour chaque classe.
Proposer un test de l’hypothèse nulle $H0 = {mi = mj } contre H1 = {mi ⁄= mj },$ fondé sur la loi de Student, et l’utiliser aﬁn de comparer quelques couples de moyennes.

3 Régression linéaire

3.1 Les données

L’exemple utilisé a déjà été proposé dans les exercices sur le chapitre 2. On dispose de données concernant l’âge (X₁), le kilométrage en milliers de kms (X₂), et le prix en milliers d’euros (Y ) pour un échantillon de voitures d’occasion d’un même type:

|----|-----|-----| |X1--|-X2--|-Y---| | 5 | 92 | 7.8 | | 4 | 64 | 9.5 | | | | | | 6 |124 | 6.4 | | 5 | 97 | 7.5 | | 5 | 79 | 8.1 | | 5 | 76 | 9.0 | | 6 | 93 | 6.1 | | | | | | 6 | 63 | 8.7 | | 2 | 13 |15.4 | | 7 |111 | 6.4 | | 7 |143 | 4.4 | |----|----------- | |

Si l’on veut étudier la liaison entre les variables potentiellement explicatives et la variable à expliquer (le prix Y ), on peut commencer par visualiser les nuages de points (X₁,Y ) et (X₂,Y ):

donnees = [
5 92 7.8; 4 64 9.5; 6 124 6.4; 5 97 7.5; 5 79 8.1; 5 76 9.0;
6 93 6.1; 6 63 8.7; 2 13 15.4; 7 111 6.4; 7 143 4.4];
// definition des variables
X1=donnees(:,1); X2=donnees(:,2); Y=donnees(:,3);
// allures des nuages
xbasc();
subplot(121), plot2d(X1,Y,-2,"111","Age",[1,0,8,16]);
subplot(122), plot2d(X2,Y,-3,"111","Kms",[0,0,150,16]);

Question 11 Que suggèrent les deux nuages?

3.2 Régression linéaire simple

On propose tout d’abord le modèle M1:

Y = β + α1X1 + 𝜀,

où 𝜀 est un échantillon de 𝒩(0,σ²).

Question 12 L’âge a-t-il une inﬂuence signiﬁcative sur le prix?.

La fonction regression (regression(M,Y,alpha)) calcule, à partir de la donnée du modèle sous la forme Y = M𝜃 + 𝜀, où Y est la variable à expliquer et M la matrice déﬁnissant le modèle, les éléments usuels:

L’estimateur du paramètre 𝜃: theta = (M^tM)⁻¹M^tY.
L’estimateur de σ²: sigma2 = MSE.
La table d’analyse de la variance pour l’hypothèse nulle “aucun régresseur n’est signiﬁcatif”.
La valeur critique associée au niveau alpha.

Cliquer sur le lien ci-dessus pour obtenir le code de la fonction, le sauvegarder sous le nom regression.sce dans le répertoire où vous utilisez scilab. Pour charger la fonction, utiliser la commande: getf 'regression.sce'. Pour appliquer la fonction, utiliser la commande:

regression(M,Y,alpha)

Pour le modèle, où on explique le prix en fonction de l’âge, on utilise les commaneds suivantes:

// MODELE (M1) Y = theta1(0) + theta1(1) X1 + eps
M1 = [ones(X1) X1]; // matrice des regresseurs
// Pour obtenir les estimations de theta, sigma2,
// et la p-valeur
[theta1,sigma21,pval1]=regression(M1,Y,alpha);

On représente enﬁn la droite de régression sur le nuage:

xbasc();
// trace de la droite
z=1:0.1:8; t=theta1(1)+theta1(2)*z; xbasc();
plot2d(X1,Y,-2,"111","Age",[1,0,8,16]);
plot2d(z,t,[1],"000");

3.3 Régression linéaire multiple

On propose d’essayer à présent le modèle “complet”

Y = γ + α1X1 + α2X2 + 𝜀.

Question 13 Donner, à l’aide de la fonction regression des estimations sans biais des paramètres (γ,α₁,α₂) et σ².

Question 14 Test de l’utilité d’un régresseur: Tester, au niveau 5%, l’hypothèse nulle “le kilométrage n’a pas d’eﬀet sur le prix” contre “c’est faux”.

La fonction regression n’eﬀectue pas le test demandé (quel test eﬀectue-t-elle?). Ecrire les commandes Scilab permettant de calculer la statistique et la p-valeur du test demandé, et conclure.

Question 15 A l’aide de la question précédente, déterminer quel sous-modèle, parmi les deux sous-modèles à un seul régresseur, explique le mieux le prix?