Tests de comparaison pour l’augmentation du volume de précipitation

(last modiﬁcation date: October 10, 2017)

Version pdf de ce document
Version sans bandeaux

1 Peut-on augmenter la quantité de pluie?
2 Un test de comparaison non paramétrique
3 Modèle gaussien
4 Et si on testait si X_i et Y _j ont même loi

1 Peut-on augmenter la quantité de pluie?
2 Un test de comparaison non paramétrique
3 Modèle gaussien
3.1 Comparaison des moyennes
3.2 (Facultatif) Comparaison des variances des 2 échantillons gaussiens
4 Et si on testait si X_i et Y _j ont même loi
4.1 Un peu de théorie
4.2 La pratique
4.3 En Scilab

1 Peut-on augmenter la quantité de pluie?

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 solidiﬁer sur les cristaux de glace qui grossissent et précipitent vers le sol en se réchauﬀant. En eﬀet, 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 artiﬁciellement 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

Expérience i	1	2	3	4	5	6	7	8
Date T_i	0	1	3	4	6	9	18	25
Ensemencement E_i	0	1	1	0	1	0	0	0
Pluies Z_i	12.85	5.52	6.29	6.11	2.45	3.61	0.47	4.56

Expérience i	9	10	11	12	13	14	15	16
Date T_i	27	8	29	32	33	35	38	39
Ensemencement E_i	0	1	1	1	0	1	1	0
Pluies Z_i	6.35	5.06	2.76	4.05	5.74	4.84	11.86	4.45

Expérience i	17	18	19	20	21	22	23	24
Date T_i	53	55	56	59	65	68	82	83
Ensemencement E_i	0	1	0	1	1	0	1	0
Pluies Z_i	3.66	4.22	1.16	5.45	2.02	0.82	1.09	0.28

T=[0,1,3,4,6,9,18,25,27,8,29,32,33,35,38,39,53,55,56,59,65,68,82,83];
E=[0,1,1,0,1,0,0,0,0,1,1,1,0,1,1,0,0,1,0,1,1,0,1,0];
Z=[12.85,5.52,6.29,6.11,2.45,3.61,0.47,4.56,6.35,5.06,2.76,4.05,5.74];
Z=[Z,4.84,11.86,4.45,3.66,4.22,1.16,5.45,2.02,0.82,1.09,0.28];

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 eﬃcace, le procédé permettrait d’augmenter la production.

Question 1 Qu’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 X₁,…,X_n₁ les hauteurs de pluies avec ensemencement, et Y ₁,…,Y _n₂ les hauteurs de pluies sans ensemencement.

X=Z(E==1);
Y=Z(E==0);

On suppose que les variables aléatoires Y ₁,…,Y _n₂ sont indépendantes et de même fonction de répartition F. On suppose que les variables aléatoires X₁,…,X_n₁ sont indépendantes et de même fonction de répartition F_Δ = F(⋅− Δ). Le but est donc de tester l’hypothèse H′₀ = {Δ ≤ 0}, ou plus exactement H₀ = {Δ = 0}, contre son alternative H₁ = {Δ > 0}. Remarquons que l’on ne fait pas d’hypothèse sur la loi de Y _i. 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éﬁnit

∑ ∑n1 U = 1{Yj<Xi } et W = Ri, 1≤i≤n1,1≤j≤n2 i=1

où R_i est le rang de X_i dans le réordonnement croissant de la population totale (X,Y ). On rappelle la relation entre W et U :

n1(n1-+-1)- U = W − 2 .

De plus les lois sous H₀ de U et W ne dépendent pas de la fonction de répartition F.

Question 2 Quelles sont les valeurs maximales et minimales de W et U.

Sous H₀, la loi de (R₁,…,R_n₁) est la loi de (σ₁,…,σ_n₁), où la permutation aléatoire σ = (σ₁,…,σ_n₁,σ_n₁+1,…,σ_n₁+n₂) suit la loi uniforme sur l’ensemble des permutations de {1,…,n₁ + n₂}.

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 :

function[U]=simulation_loi_U(n1,n2,N)
W=simulation_loi_W(n1,n2,N);
U=W-0.5*n1*(n1+1);
endfunction;

On rappelle que sous H₀, on a

1- -1- 𝔼 [U ] = 2 n1n2 et V ar(U ) = 12 n1n2 (n1 + n2 + 1 ).

De plus :

1 -----U--−-2-n1n2------- loi min(nl1im,n2)→∞ ∘ -1------------------ = 𝒩 (0,1). 12 n1n2 (n1 + n2 + 1)

On trace la densité gaussienne sur l’histogramme de l’approximation de la loi de U aﬁn 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 H₁, en revanche, les variables aléatoires X_i prennent en moyenne des valeurs plus grandes que les variables aléatoires Y _i. Bien que la loi de U sous H₁ ne soit pas accessible, car elle dépend de F, on remarque que W prend en moyenne des valeurs plus grandes sous H₁ que sous H₀. En particulier, on peut montrer que sous H₁

1 lim ∘----U--−-2-n1n2------- p=.s. +∞. min(n1,n2)→ ∞ -1n n (n + n + 1) 12 1 2 1 2

Question 3 Construire un test convergent à partir de la statistique U pour tester H₀ contre H₁. Déterminer la région critique associée à ce test de niveau α = 5%. On utilisera l’approche asymptotique et l’approche directe. Dans ce dernier cas, comme la fonction de répartition de U n’est pas explicite, on utilisera la fonction cdfUX qui fournit par simulation les quantiles de 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 4 Calculer 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 IC_Q l’intervalle de conﬁance 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 5 Calculer la p-valeur pour les données du problème, à l’aide de l’approximation gaussienne. Conclusion.

Remarquons que ne pas rejeter H₀ ne signiﬁe pas que H₀ soit vraie. Ainsi si on considère H′₀ = {Δ ≤ 0}, au lieu de H₀, on obtient la même région critique pour la statistique de test U. En particulier on ne peut pas distinguer H₀ de H₀′ et en fait même H₀ de H₀′′ = {ℙ(X₁ ≥ x) ≤ ℙ(Y ₁ ≥ x) ∀x ∈ ℝ} ou de H₀′′′ = {ℙ(X₁ ≤ Y ₁) ≤ 1∕2}.

3 Modèle gaussien

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

X1,...,Xn i.i.d. 𝒩 (μ1, σ2), 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?

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

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

1 ∑n1 1 ∑n2 ¯X = --- Xi, Y¯ = --- Yi, n1 i=1 n2 i=1

ainsi que les variances empiriques

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

// 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 2
(n1-−-1)S-X
σ2 + 2
(n2-−-1)S-Y
σ2 est la loi du χ² à n₁ + n₂ − 2 degrés de liberté.

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

∘ -n1n2--- ( ¯X − Y¯) T = --------∘-------2--------2- n1 + n2 (n1−-1)SnX++n-(n−2−21)SY 1 2

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 6 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 7 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 ([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 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 9 Reprendre les questions 6 et 7, si l’hypothèse nulle est {μ₁ ≤ μ₂}.

3.2 (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 10 (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 11 (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.

4 Et si on testait si X_i et Y _j 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

1--n∑1 1-∑n2 F Xn1(x ) = n 1{Xi≤x} et F Yn2(x ) = n 1{Yi≤x}. 1 i=1 2 i=1

On déﬁnit la statistique du test

∘ -------- -n1n2--- ζn1,n2 = n1 + n2 sux∈pℝ |F Xn1(x) − F Yn2(x)|.

Si l’hypothèse nulle H₀ = {Loi(X_i) = Loi(Y _i)} est vraie alors ζ_n₁,n₂ 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

+∑∞ 22 K (s) = (− 1)ke−2ks s > 0. k=−∞

Si l’hypothèse nulle est fausse alors ζ_n₁,n₂ diverge presque sûrement.

4.2 La pratique

Utilisation de ce résultat :

on dispose de l’observation des échantillons X₁ = x₁,…,X_n₁ = x_n₁ et Y ₁ = y₁,…,Y _n₂ = y_n₂
on calcule les fonctions de répartition empiriques FX_n₁ et FY _n₂ à partir des échantillons
on calcule la valeur de ζ_n₁,n₂ = z_n₁,n₂
on calcule la p-valeur α_n₁,n₂ = ℙ(S ≥ z_n₁,n₂).
si α_n₁,n₂ prend des valeurs faibles inférieures à α = 5% typiquement alors on rejette l’hypothèse H₀ sinon on l’accepte.
α représente le niveau de conﬁance 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 ﬁchier 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 : ζ_n₁,n₂∕
q : ℙ = ℙ(S ≥ ζ_n₁,n₂).

Question 12 Le procédé d’ensemencement proposé change-t-il signiﬁcativement la loi des hauteurs de pluie?