Fermer X

Contents

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 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

Expérience i 1 2 3 4 5 6 7 8
Date Ti 0 1 3 4 6 9 18 25
Ensemencement Ei0 1 1 0 1 0 0 0
Pluies Zi 12.855.526.296.112.453.610.474.56
Expérience i 9 10 11 12 13 14 15 16
Date Ti 27 8 29 32 33 35 38 39
Ensemencement Ei0 1 1 1 0 1 1 0
Pluies Zi 6.355.062.764.055.744.8411.864.45
Expérience i 17 18 19 20 21 22 23 24
Date Ti 53 55 56 59 65 68 82 83
Ensemencement Ei0 1 0 1 1 0 1 0
Pluies Zi 3.664.221.165.452.020.821.090.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 efficace, 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 X1,,Xn1 les hauteurs de pluies avec ensemencement, et Y 1,,Y n2 les hauteurs de pluies sans ensemencement.

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

On suppose que les variables aléatoires Y 1,,Y n2 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 H0 = {Δ 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 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éfinit

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

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 :

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

De plus les lois sous H0 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 H0, la loi de (R1,,Rn1) est la loi de (σ1,n1), où la permutation aléatoire σ = (σ1,n1n1+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);


PIC

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 H0, 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 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;  

PIC

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 Y i. 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

                      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 H0 contre H1. 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 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 5 Calculer la p-valeur pour les données du problème, à l’aide de l’approximation gaussienne. Conclusion.

Remarquons que ne pas rejeter H0 ne signifie pas que H0 soit vraie. Ainsi si on considère H0 = {Δ 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 H0et en fait même H0 de H0′′ = {(X1 x) (Y 1 x) x } ou de H0′′′ = {(X1 Y 1) 12}.

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

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

  1. Les données peuvent-elles être considérées comme des réalisations de lois gaussienne pour les paramètres (moyenne, variance) appropriés?
  2. 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

      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 (n1 1)SX2∕σ2 (resp. (n 2 1)SY 2∕σ2) sont indépendantes et de loi 𝒩(μ12∕n 1) (resp. 𝒩(μ22∕n 2)) et χ2(n 1 1) (resp. χ2(n 2 1)). Comme SX2 et S Y 2 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 χ2 à n 1 + n2 2 degrés de liberté.

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

    ∘  -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 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)→∞∘ --------
   -n1n2---
   n1 + n2 = +, on en déduit que sous H1, la statistique T diverge vers + quand min(n1,n2) →∞.

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éfinie 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

  1. Faire varier Y , les autres valeurs étant fixées. Pour quelles valeurs de Y rejetez vous H0?
  2. Avec la valeur de Y ainsi déterminée, qu’observez vous si SY 2 augmente, si S Y 2 diminue?
  3. En reprenant les valeurs numériques du problème, qu’observez vous si n1 et n2 augmentent, si n1 et n2 diminuent jusqu’aux valeurs limites?

Question 9 Reprendre les questions 6 et 7, si l’hypothèse nulle est {μ1 μ2}.

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 diffé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, 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∕S Y 2 sous l’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é 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 Xi 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éfinit la statistique du test

        ∘ --------
           -n1n2---
ζn1,n2 =    n1 + n2 sux∈pℝ |F Xn1(x) − F Yn2(x)|.
Si l’hypothèse nulle H0 = {Loi(Xi) = Loi(Y i)} 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
         +∑∞            22
K (s) =      (− 1)ke−2ks   s > 0.
        k=−∞
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 Y 1 = y1,,Y n2 = yn2
  • on calcule les fonctions de répartition empiriques FXn1 et FY n2 à 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 :

  1. D : ζn1,n2∘ ------
  -n1n2-
  n1+n2
  2. q : (              )
      ∘ -n1n2-
 S ≥    n1+n2D = (S ζn1,n2).

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

L'École des Ponts ParisTech est membre fondateur de

L'École des Ponts ParisTech est certifiée 

ParisTech ParisEst ParisTech

 

Copyright 2014
École des Ponts ParisTech
All rights reserved