Chacun sait que l’eau gèle en dessous de 0°C. Ce que l’on sait moins c’est que la phase liquide
peut exister en-dessous de cette limite. Quand elle reste liquide en dessous de 0°C on dit que l’eau
est surfondue. Ainsi les nuages sont presque toujours constitués de gouttelettes d’eau surfondue
coexistant avec quelques cristaux de glace. C’est un état d’équilibre instable et tout apport
supplémentaire de cristaux de glace ou d’iodure d’argent provoque l’évaporation rapide
des gouttelettes en surfusion. La vapeur ainsi libérée va aussitôt se solidifier sur les
cristaux de glace qui grossissent et précipitent vers le sol en se réchauffant. En effet, aux
températures négatives, l’air atmosphérique est saturé de vapeur d’eau par rapport à
la glace avant de l’être par rapport à l’eau surfondue. C’est cette propriété qui est
exploitée pour forcer les précipitations. Ces pluies artificiellement provoquées peuvent
être intéressantes pour l’agriculture comme en témoignent les nombreuses expériences
d’ensemencement réalisées dans le monde au cours des années 1950 à 1970. C’est ainsi qu’entre le
1er juin et le 23 août 1975 (83 jours), une expérience visant à accroître le volume des
précipitations s’est tenue en Floride (USA). Les données suivantes sont extraites de
l’article de W.L. Woodley, J. Simpson, R. Biondini et J. Berkeley (1977) : Rainfall
Results 1970-75 : Florida Area Cumulus Experiment, publié dans Science, vol. 195, pp.
735-742.
Données d’ensemencement de nuages par iodure d’argent
Les sentiments mitigés des agriculteurs peuvent se résumer ainsi :
le procédé est cher ;
jusqu’à ce jour nous nous en sommes passés et nous pouvons continuer à nous en
passer ;
cependant si l’ensemencement de nuages par iodure d’argent est efficace, le procédé
permettrait d’augmenter la production.
Question 1Qu’attendent les paysans du statisticien?
Les paragraphes qui suivent ont pour but de fournir la réponse du statisticien à l’aide
d’une approche non-paramétrique (voir le livre de J. Bernier, E. Parent, JJ. Boreux,
2000, Stat pour l’environnement, Tec et Doc, où cet exemple est traité dans un cadre
paramétrique).
2 Un test de comparaison non paramétrique
On note X1,…,Xn1 les hauteurs de pluies avec ensemencement, et Y1,…,Yn2 les hauteurs de pluies
sans ensemencement.
X=Z(E==1); Y=Z(E==0);
On suppose que les variables aléatoires Y1,…,Yn2 sont indépendantes et de même fonction de
répartition F. On suppose que les variables aléatoires X1,…,Xn1 sont indépendantes et de
même fonction de répartition FΔ = F(⋅− Δ). Le but est donc de tester l’hypothèse
H′0 = {Δ ≤ 0}, ou plus exactement H0 = {Δ = 0}, contre son alternative H1 = {Δ > 0}.
Remarquons que l’on ne fait pas d’hypothèse sur la loi de Yi. Il existe beaucoup de
tests pour ce genre de problème. On se propose d’utiliser le test de Mann et Whitney
(1947) qui repose sur les rangs de l’échantillon X dans l’échantillon des deux populations
(X,Y ).
On définit
où Ri est le rang de Xi dans le réordonnement croissant de la population totale (X,Y ). On
rappelle la relation entre W et U :
De plus les lois sous H0 de U et W ne dépendent pas de la fonction de répartition
F.
Question 2Quelles sont les valeurs maximales et minimales de W et U.
Sous H0, la loi de (R1,…,Rn1) est la loi de (σ1,…,σn1), où la permutation aléatoire
σ = (σ1,…,σn1,σn1+1,…,σn1+n2) suit la loi uniforme sur l’ensemble des permutations de
{1,…,n1 + n2}.
Il existe des tables pour la loi de W. Toutefois, il est aussi rapide d’en donner une
approximation par simulation. On peut par exemple donner une estimation de ℙ(W > w) à l’aide
de simulations de N variables aléatoires indépendantes de loi W. On commence par simuler la loi
de W :
//n1 longueur de la population X //n2 longueur de la population Y //N nombre de simulations function[W]=simulation_loi_W(n1,n2,N) //longueur de la population totale n=n1+n2; vect=[1:n]'; //matrice dont chaque colonne est une permutation de {1,...,n} R=grand(N,'prm',vect); //on prend les n1 premierséléments de chaque permutation R=R([1:n1],:); W=sum(R,'r'); endfunction;
On fait N = 10000 simulations de W, W=simulation_loi_W(n1,n2,10000);, et on trace un
histogramme (graphique 1) des valeurs obtenues pour visualiser une approximation de la loi de W
:
xbasc();
histplot([0:2:n1*n2+n1*(n1+1)/2],W);
Figure 1: Simulation de la loi de W
De même on simule la loi de U en se basant sur la relation entre U et W :
On trace la densité gaussienne sur l’histogramme de l’approximation de la loi de U afin de
visualiser la convergence en loi ci-dessus. Pour obtenir le graphique 2, on utilise la fonction
Scilab
xbasc(); //nettoie la fenêtre graphique
comparaison_U_Normale(n1,n2,N);
avec par exemple N = 10000 simulations de U.
function comparaison_U_Normale(n1,n2,N) U=simulation_loi_U(n1,n2,N); t=(U-0.5*n1*n2)/sqrt(n1*n2*(n1+n2+1)/12); x=[-4:2/sqrt(n1*n2*(n1+n2+1)/12):4]; histplot(x,t); z=[-4:0.1:4]; y=exp(-z.^2/2)/sqrt(2*%pi); plot2d(z,y); endfunction;
Figure 2: Comparaison entre la loi de U recentrée et renormalisée et une gaussienne 𝒩(0, 1).
Sous H1, en revanche, les variables aléatoires Xi prennent en moyenne des valeurs plus grandes
que les variables aléatoires Yi. Bien que la loi de U sous H1 ne soit pas accessible, car elle dépend
de F, on remarque que W prend en moyenne des valeurs plus grandes sous H1 que sous H0. En
particulier, on peut montrer que sous H1
Question 3Construire un test convergent à partir de la statistique U pour tester H0contreH1. Déterminer la région critique associée à ce test de niveau α = 5%. On utilisera l’approcheasymptotique et l’approche directe. Dans ce dernier cas, comme la fonction de répartition deU n’est pas explicite, on utilisera la fonction cdfUXqui fournit par simulation les quantilesde la loi de U.
Cette fonction calcule U à partir de deux échantillons X et Y (elle ne tient pas compte des
ex-aequo).
//calcul de U function[U]=Mann_Whitney(X,Y) n1=length(X); //classe est la population totale (X,Y) triée par ordre croissant //permut donne les rangs deséléments du vecteur classe dans //la population (X,Y) [classe,permut]=gsort([X Y],'g','i'); //rangs stocke les rangs deséléments de X //dans la population triée rangs=find(permut<=n1); //on fait la somme de ces rangs, cela donne W W=sum(rangs); //on calcule U U=W-n1*(n1+1)/2; endfunction;
Voici la fonction qui calcule par simulation les quantiles de la loi de U :
//renvoie x le quantile d'ordre p de la loi U : p=P(U<x) //IC intervalle de confianceà 5% sur x //N nombre de simulations pour estimer x //n1 et n2 paramètres de U function[x,IC]=cdfUX(n1,n2,p,N) //simule N va de loi U U=simulation_loi_U(n1,n2,N); //tri par ordre croissant [s,k]=gsort(U,'g','i'); //estimateur x=s(floor(p*N)); inf_IC=s(floor(p*N-sqrt(N)*1.96*sqrt(p*(1-p)))); sup_IC=s(floor(p*N+sqrt(N)*1.96*sqrt(p*(1-p)))); IC=[inf_IC,sup_IC]; endfunction;
Question 4Calculer la p-valeur associée au test précédent pour les données du problème,à l’aide de la fonction Scilab ci-dessous cdfU. Conclusion.
La fonction [P,Q,IC_Q]=cdfU(u,n1,n2,N) renvoie une estimation par simulation de
P = ℙ(𝕌 ≤ ≊), Q = 1 −P, ainsi que ICQ l’intervalle de confiance asymptotique sur Q à 95%. On
prendra N = 10000.
//n1 et n2 cardinaux deséchantillons //N nbombre de simulations //P=P(U<=u) //Q=1-P //IC_Q intervalle de confiance sur Q function[P,Q,IC_Q]=cdfU(u,n1,n2,N) U=simulation_loi_U(n1,n2,N); u=u*ones(U); P=sum(U<=u)/N; Q=sum(U>u)/N; IC_Q=[Q-1.96*sqrt(Q*(1-Q))/sqrt(N),Q+1.96*sqrt(Q*(1-Q))/sqrt(N)]; endfunction;
Question 5Calculer la p-valeur pour les données du problème, à l’aide de l’approximationgaussienne. Conclusion.
Remarquons que ne pas rejeter H0 ne signifie pas que H0 soit vraie. Ainsi si on considère
H′0 = {Δ ≤ 0}, au lieu de H0, on obtient la même région critique pour la statistique
de test U. En particulier on ne peut pas distinguer H0 de H0′ et en fait même H0 de
H0′′ = {ℙ(X1≥ x) ≤ ℙ(Y1≥ x) ∀x ∈ ℝ} ou de H0′′′ = {ℙ(X1≤ Y1) ≤ 1∕2}.
3 Modèle gaussien
On suppose que les volumes de pluies suivent une loi gaussienne et que l’effet de l’ensemencement
ne modifie 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
On souhaite tester l’hypothèse nulle H0: “le procédé n’a pas d’effet” contre “le procédé produit
une augmentation significative de la quantité de pluie”, autrement dit
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?
On peut répondre au point 1. à l’aide du test de Kolmogorov quand la moyenne et la variance sont
connues a priori. Sinon, plusieurs tests de normalité sont détaillés dans la littérature. Le point 2.
est abordé dans le paragraphe facultatif 3.2.
Nous admettons donc dans un premier temps ces deux hypothèses.
3.1 Comparaison des moyennes
Afin 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
ainsi
que les variances empiriques
// 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 (n1− 1)SX2∕σ2 (resp. (n2− 1)SY2∕σ2)
sont indépendantes et de loi 𝒩(μ1,σ2∕n1) (resp. 𝒩(μ2,σ2∕n2)) et χ2(n1− 1) (resp. χ2(n2− 1)).
Comme SX2 et SY2 sont indépendants, on en déduit que la loi de + est la
loi du χ2 à n1 + n2− 2 degrés de liberté.
Sous H0, la loi de est la loi gaussienne centrée réduite. On en déduit que
sous H0,
suit une loi de Student de paramètre n1 + n2− 2.
Comme sous H1, X−Y converge p.s. vers μ1− μ2> 0 quand min(n1,n2) →∞, et que
lim min(n1,n2)→∞ = +∞, on en déduit que sous H1, la statistique T diverge vers +∞
quand min(n1,n2) →∞.
Question 6Construire, à partir de la statistique T, un test d’égalité des moyennes endéterminant sa région critique, la valeur critique au niveau α et la p-valeur.
Question 7En utilisant la fonction ttest2définie ci-dessous, calculer numériquementla valeur de la statistique de test, la valeur critique au niveau 5%, et la p-valeur pour lesdonnées de pluie. Conclusion.
La fonction ttest2 ([Tobs,vc,pval]=ttest2(n1,Xbar,SX2,n2,Ybar,SY2,alpha)) retourne la
valeur de la statistique T observée, Tobs, la valeur critique associée au niveau α, vc, et la p-valeur,
pval, 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 8
Faire varierY , les autres valeurs étant fixées. Pour quelles valeurs deY rejetez vousH0?
Avec la valeur deY ainsi déterminée, qu’observez vous si SY2augmente, si SY2diminue?
En reprenant les valeurs numériques du problème, qu’observez vous si n1et n2augmentent, si n1et n2diminuent jusqu’aux valeurs limites?
Question 9Reprendre les questions 6et 7, si l’hypothèse nulle est {μ1≤ μ2}.
Pour s’assurer du point (2), on peut tout d’abord supposer que ces variances sont différentes, i.e.
puis construire un test d’égalité des variances, afin de voir si la différence sur les variances (ou les
écart-types) observée est significative.
Question 10 (Facultatif)Quelle est la loi de la statistique F = SX2∕SY2sousl’hypothèse nulle d’égalité des variances? Quel est son comportement sous H1, quand n1→∞et n2→∞?
Question 11 (Facultatif)Construire, à partir de la statistique F, un test d’égalité desvariances. 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.
4 Et si on testait si Xi et Yj ont même loi
Ce test est une variante du test de Kolmogorov Smirnov pour deux échantillons indépendants.
4.1 Un peu de théorie
On dispose aussi dans ce cas d’un test de Kolmogorov Smirnov. Il s’agit de comparer les fonctions
de répartition empirique
On définit la statistique du test
Si
l’hypothèse nulle H0 = {Loi(Xi) = Loi(Yi)} est vraie alors ζn1,n2 tend en loi vers une variable
aléatoire S dont la loi est indépendante de celles de X et Y . De plus on connaît sa fonction de
répartition
Si
l’hypothèse nulle est fausse alors ζn1,n2 diverge presque sûrement.
4.2 La pratique
Utilisation de ce résultat :
on dispose de l’observation des échantillons X1 = x1,…,Xn1 = xn1 et Y1 = y1,…,Yn2 =
yn2
on calcule les fonctions de répartition empiriques FXn1 et FYn2 à partir des
échantillons
on calcule la valeur de ζn1,n2 = zn1,n2
on calcule la p-valeur αn1,n2 = ℙ(S ≥ zn1,n2).
si αn1,n2 prend des valeurs faibles inférieures à α = 5% typiquement alors on rejette
l’hypothèse H0 sinon on l’accepte.
α représente le niveau de confiance du test.
4.3 En Scilab
La fonction ks2 ([D,q]=ks2(X,Y)) permet de réaliser le test précédemment décrit. Pour utiliser
cette fonction, cliquer sur le lien ci-dessus et enregistrer el fichier sous le nom ks2.sce dans le
répertoire dans lequel vous travaillez. Pour charger la fonction, utiliser la commande :
exec 'ks2.sce'.
Explications des paramètres :
x : vecteur ligne des observations du premier échantillon
y : vecteur ligne des observations du second échantillon
Cette fonction renvoie :
D : ζn1,n2∕
q : ℙ = ℙ(S ≥ ζn1,n2).
Question 12Le procédé d’ensemencement proposé change-t-il significativement la loi deshauteurs de pluie?