Cours complet pour débutants pour apprendre la programmation en R

C. Planification d'une simulation en R▲

Il existe de multiples façons de réaliser la mise en œuvre informatique d'une simulation, mais certaines sont plus efficaces que d'autres. Cette annexe passe en revue diverses façons de faire des simulations avec R à l'aide d'un exemple simple de nature statistique.

C-1. Contexte▲

Soit kitxmlcodeinlinelatexdvpX_1, \dots, X_nfinkitxmlcodeinlinelatexdvp un échantillon aléatoire tiré d'une population distribuée selon une loi uniforme sur l'intervalle kitxmlcodeinlinelatexdvp\left(\theta-\frac{1}{2}, \theta + \frac{1}{2} \right )finkitxmlcodeinlinelatexdvp. On considère trois estimateurs sans biais du paramètre inconnu :

La moyenne arithmétique

kitxmlcodeinlinelatexdvp\hat{\theta_1} = \frac{1}{n}\sum_{i=1}^n X_ifinkitxmlcodeinlinelatexdvp
La médiane empirique

kitxmlcodeinlinelatexdvp\hat{\theta_2} = \left\{\begin{matrix} X_{\left(\frac{n+1}{2} \right )}, & n$ impair$\\ \frac{1}{2}\left( X_{\left(\frac{n}{2} \right )} + X_{\left(\frac{n}{2}+1 \right )}\right ), & n$ pair$ \end{matrix}\right.finkitxmlcodeinlinelatexdvp

où kitxmlcodeinlinelatexdvpX_{(k)}finkitxmlcodeinlinelatexdvp est la k^e statistique d'ordre de l'échantillon aléatoire ;
La mi-étendue kitxmlcodeinlinelatexdvp\hat{\theta_3} = \frac{X_{(1)}+X_{(n)}}{2}finkitxmlcodeinlinelatexdvp

À l'aide de la simulation on veut, d'une part, vérifier si les trois estimateurs sont bel et bien sans biais et, d'autre part, déterminer lequel a la plus faible variance.

Pour ce faire, on doit d'abord simuler un grand nombre d'échantillons aléatoires de taille kitxmlcodeinlinelatexdvpnfinkitxmlcodeinlinelatexdvp d'une distribution kitxmlcodeinlinelatexdvpU\left(\theta-\frac{1}{2}, \theta + \frac{1}{2} \right )finkitxmlcodeinlinelatexdvp pour une valeur de kitxmlcodeinlinelatexdvp\thetafinkitxmlcodeinlinelatexdvp choisie. Pour chaque échantillon, on calculera ensuite les trois estimateurs ci-dessus, puis la moyenne et la variance, par type d'estimateur, de tous les estimateurs obtenus. Si la moyenne des kitxmlcodeinlinelatexdvpNfinkitxmlcodeinlinelatexdvp estimateurs kitxmlcodeinlinelatexdvp\hat{\theta_i},\ i=1,2,3finkitxmlcodeinlinelatexdvp est près de kitxmlcodeinlinelatexdvp\thetafinkitxmlcodeinlinelatexdvp alors on pourra conclure que kitxmlcodeinlinelatexdvp\hat{\theta_i}finkitxmlcodeinlinelatexdvp est sans biais. De même, on déterminera lequel des trois estimateurs a la plus faible variance selon le classement des variances empiriques.

C-2. Première approche : avec une boucle▲

La façon la plus intuitive de mettre en œuvre cette étude de simulation en R consiste à utiliser une boucle for. Avec cette approche, il est nécessaire d'initialiser une matrice de 3 lignes et kitxmlcodeinlinelatexdvpNfinkitxmlcodeinlinelatexdvp colonnes (ou l'inverse) dans laquelle seront stockées les valeurs des trois estimateurs pour chaque simulation. Une fois la matrice remplie dans la boucle, il ne reste plus qu'à calculer la moyenne et la variance par ligne pour obtenir les résultats souhaités.

La figure C.1 présente un exemple de code adéquat pour réaliser la simulation à l'aide d'une boucle.

Si l'on souhaite pouvoir exécuter le code de la figure C.1 facilement à l'aide d'une seule expression, il suffit de placer l'ensemble du code dans une fonction. La fonction simul1 de la figure C.2 reprend le code de la figure C.1, sans les commentaires. On a alors :

Sélectionnez

1.
2.
3.
4.
5.
6.
7.
8.

> simul1(10000, 100, 0)
$biais
     Moyenne       Mediane   Mi-etendue
6.945790e-05 -5.293742e-04 9.199309e-05

$variances
     Moyenne      Mediane   Mi-etendue
8.128808e-04 2.396603e-03 4.918073e-05

C-3. Seconde approche : avec sapply▲

On le sait, les boucles sont inefficaces en R. Il est en général plus efficace de déléguer les boucles aux fonctions lapply et sapply (section 7.37.3 Pour en savoir plus). On rappelle que la syntaxe de ces fonctions est

Fig. C.1 - Code pour la simulation utilisant une boucle for

Sélectionnez

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.

## Bonne habitude à prendre : stocker les constantes dans
## des variables faciles à modifier au lieu de les écrire
## explicitement dans le code.
size <- 100             # taille de chaque échantillon
nsimul <- 10000         # nombre de simulations
theta <- 0              # la valeur du paramètre

## Les lignes ci-dessous éviteront de faire deux additions
## 'nsimul' fois.
a <- theta - 0.5        # borne inférieure de l'uniforme
b <- theta + 0.5        # borne supérieure de l'uniforme

## Initialisation de la matrice dans laquelle seront
## stockées les valeurs des estimateurs. On donne également
## des noms aux lignes de la matrice afin de facilement
## identifier les estimateurs.
x <- matrix(0, nrow = 3, ncol = nsimul)
rownames(x) <- c("Moyenne", "Mediane", "Mi-etendue")

## Simulation comme telle.
for (i in 1:nsimul)
{
    u <- runif(size, a, b)
    x[, i] <- c(mean(u),  # moyenne
                median(u) # médiane
                mean(range(u))) # mi-étendue
}

## On peut maintenant calculer la moyenne et la variance
## par ligne.
rowMeans(x) - theta     # vérification du biais
apply(x, 1, var)        # comparaison des variances

Fig. C.2 - Définition de la fonction simul1

Sélectionnez

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.

simul1 <- function(nsimul, size, theta)
{
    a <- theta - 0.5
    b <- theta + 0.5
    x <- matrix(0, nrow = 3, ncol = nsimul)
    rownames(x) <- c("Moyenne", "Mediane", "Mi-etendue")
    for (i in 1:nsimul)
    {
        u <- runif(size, a, b)
        x[, i] <- c(mean(u), median(u), mean(range(u)))
    }
    list(biais = rowMeans(x) - theta,
         variances = apply(x, 1, var))
}

Sélectionnez

1.
2.

lapply(x, FUN, ...)
sapply(x, FUN, ...)

Ces fonctions appliquent la fonction FUN à tous les éléments de la liste ou du vecteur x et retournent les résultats sous forme de liste (lapply) ou, lorsque c'est possible, de vecteur ou de matrice (sapply). Il est important de noter que les valeurs successives de x seront passées comme premier argument à la fonction FUN. Le cas échéant, les autres arguments de FUN sont spécifiés dans le champ ....

Pour pouvoir utiliser ces fonctions dans le cadre d'une simulation comme celle dont il est question ici, il s'agit de définir une fonction qui fera tous les calculs pour une simulation, puis de la passer à sapply pour obtenir les résultats de simulations. La figure C.3 présente une première version d'une telle fonction. On remarquera que l'argument i ne joue aucun rôle dans la fonction. Voici un exemple d'utilisation pour un petit nombre de simulations (4) :

Fig. C.3 - Définitions des fonctions fun1 et simul2

Sélectionnez

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.

fun1 <- function(i, size, a, b)
{
    u <- runif(size, a, b)
    c(Moyenne = mean(u),
      Mediane = median(u),
      "Mi-etendue" = mean(range(u)))
}
simul2 <- function(nsimul, size, theta)
{
    a <- theta - 0.5
    b <- theta + 0.5

    x <- sapply(1:nsimul, fun1, size, a, b)

    list(biais = rowMeans(x) - theta,
         variances = apply(x, 1, var))
}

Sélectionnez

1.
2.
3.
4.
5.
6.

> sapply(1:4, fun1, size = 10, a = -0.5, b = 0.5)

                   [,1]        [,2]         [,3]        [,4]
   Moyenne -0.076212566 -0.07839662 0.0006719503  0.02844156
   Mediane -0.117479511 -0.01346815 0.0403468686  0.02980248
Mi-etendue -0.005960669 -0.05199959 0.0395837893 -0.03937694

On remarque donc que les résultats de chaque simulation se trouvent dans les colonnes de la matrice obtenue avec sapply.

Pour compléter l'analyse, on englobe le tout dans une fonction simul2, dont le code se trouve à la figure C.3 :

Sélectionnez

1.
2.
3.
4.
5.
6.
7.
8.
9.

> simul2(10000, 100, 0)

$biais
     Moyenne      Mediane   Mi-etendue
0.0003148615 0.0002347348 0.0001190661

$variances
     Moyenne      Mediane   Mi-etendue
8.272098e-04 2.401754e-03 4.875476e-05

Il est généralement plus facile de déboguer le code avec cette approche puisque l'on peut rapidement circonscrire un éventuel problème à fun1 ou simul2.

C-4. Variante de la seconde approche▲

Une chose manque d'élégance dans la seconde approche : l'obligation d'inclure un argument factice dans la fonction fun1. La fonction replicate (section 6.5) permet toutefois de passer outre cette contrainte. En effet, cette fonction exécute un nombre donné de fois une expression quelconque.

Fig. C.4 - Définitions des fonctions fun2 et simul3

Sélectionnez

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.

fun2 <- function(size, a, b)
{
    u <- runif(size, a, b)
    c(Moyenne = mean(u),
      Mediane = median(u),
      "Mi-etendue" = mean(range(u)))
}

simul3 <- function(nsimul, size, theta)
{
    a <- theta - 0.5
    b <- theta + 0.5

    x <- replicate(nsimul, fun2(size, a, b))

    list(biais = rowMeans(x) - theta,
         variances = apply(x, 1, var))
}

Les fonctions fun2 et simul3 de la figure C.4 sont des versions légèrement modifiées de fun1 et simul2 pour utilisation avec replicate. On a alors

Sélectionnez

1.
2.

     Moyenne      Mediane   Mi-etendue
8.578549e-04 2.483565e-03 4.928925e-05

C-5. Gestion des fichiers▲

Pour un petit projet comme celui utilisé en exemple ici, il est simple et pratique de placer tout le code informatique dans un seul fichier de script. Pour un plus gros projet, cependant, il vaut souvent mieux avoir recours à plusieurs fichiers différents. Le présent auteur utilise pour sa part un fichier par fonction.

À des fins d'illustration, supposons que l'on utilise l'approche de la section C.4Variante de la seconde approche avec la fonction replicate et que le code des fonctions fun2 et simul3 est sauvegardé dans des fichiers fun2.R et simul3.R, dans l'ordre. Si l'on crée un autre fichier, disons go.R, ne contenant que des expressions source pour lire les autres fichiers, il est alors possible de démarrer des simulations en exécutant ce seul fichier. Dans notre exemple, le fichier go.R contiendrait les lignes suivantes :

Sélectionnez

source("fun2.R")
source("simul3.R")
simul3(10000, 100, 0)

Une simple commande

Sélectionnez

> source("go.R")

exécutera alors une simulation complète.

C-6. Exécution en lot▲

Les utilisateurs plus avancés pourront vouloir exécuter leur simulation R en lot (batch) pour en accélérer le traitement. Dans ce mode, aucune interface graphique n'est démarrée et tous les résultats sont redirigés vers un fichier pour consultation ultérieure. Pour les simulations demandant un long temps de calcul, c'est très pratique.

On exécute R en lot depuis la ligne de commande (Invite de commande sous Windows, Terminal sous OS X ou Linux). Une fois placé dans le répertoire contenant les fichiers de script, il suffit d'entrer à la ligne de commande

Sélectionnez

R CMD BATCH go.R

La sortie de cette commande (et donc tous les résultats des expressions R du fichier go.R) sera placée par défaut dans le fichier go.Rout. Sous Windows, le dossier d'installation de R peut ne pas se trouver dans la variable d'environnement %PATH%, auquel cas il faut spécifier le chemin d'accès complet de l'exécutable à la ligne de commande :

Sélectionnez

"c:\Program Files\R\R-x.y.z\bin\R" CMD BATCH go.R

Remplacer R-x.y.z par le numéro de version courant de R.

C-7. Conclusion▲

Le nombre de simulations, kitxmlcodeinlinelatexdvpNfinkitxmlcodeinlinelatexdvp et la taille de l'échantillon, kitxmlcodeinlinelatexdvpnfinkitxmlcodeinlinelatexdvp, ont tous deux un impact sur la qualité des résultats, mais de manière différente. Quand kitxmlcodeinlinelatexdvpnfinkitxmlcodeinlinelatexdvp augmente, la précision des estimateurs augmente. Ainsi, dans l'exemple ci-dessus, le biais et la variance des estimateurs de kitxmlcodeinlinelatexdvp\thetafinkitxmlcodeinlinelatexdvp seront plus faibles. D'autre part, l'augmentation du nombre de simulations diminue l'impact des échantillons aléatoires individuels et, de ce fait, améliore la fiabilité des conclusions de l'étude.

D'ailleurs, les conclusions de l'étude de simulation sur le biais et la variance des trois estimateurs de la moyenne d'une loi uniforme sont les suivantes : les trois estimateurs sont sans biais et la mi-étendue a la plus faible variance. En effet, on peut démontrer mathématiquement que, pour kitxmlcodeinlinelatexdvpnfinkitxmlcodeinlinelatexdvp impair,

kitxmlcodelatexdvpVar[\hat{\theta_1}] = \frac{1}{12n}finkitxmlcodelatexdvp kitxmlcodelatexdvpVar[\hat{\theta_2}] = \frac{1}{4n+2}finkitxmlcodelatexdvp kitxmlcodelatexdvpVar[\hat{\theta_3}] = \frac{1}{2(n+1)(n+2)}finkitxmlcodelatexdvp

et donc

kitxmlcodelatexdvpVar[\hat{\theta_3}] \le Var[\hat{\theta_1}] \le Var[\hat{\theta_2}]finkitxmlcodelatexdvp

pour tout kitxmlcodeinlinelatexdvpn \ge 2finkitxmlcodeinlinelatexdvp.