Simulations avec Splus
Voici les différents sujets abordés dans
cette rubrique:
1. Introduction:
Une des forces de S-Plus est qu'il est facile pour un utilisateur de créer ses propres fonctions et, ainsi, d'augmenter les possibilités de S-Plus. La syntaxe de base pour créer sa propre fonction est la suivante:
ma_fonction <- function( arguments )
{
instructions
}
où arguments sont les différents paramètres ou variables utilisées par la fonction pour effectuer les calculs qui lui sont demandés. Il faut séparer chaque argument par une virgule. Le bloc instructions contient les différentes commandes S-plus pour effectuer les calculs proprement dits. On tape exactement comme si on utilisait S-Plus de façon interactive, i.e. comme à l'habitude. Bien qu'il est possible de taper directement les instructions dans la fenêtre de S-Plus, il est fortement recommandé d'utiliser un éditeur de texte et, ensuite, de faire du "cut and paste". En effet, ceci permet de corriger les erreurs beaucoup plus facilement, car l'éditeur de texte par défaut de S-Plus n'est pas fameux et de garder une copie de nos fonctions dans un autre fichier.
Par exemple, parmi les fonctions de S-Plus, on ne retrouve pas de fonction pour calculer le coefficient d'asymétrie ("skewness"). Alors, nous avons simplement à la programmer en tapant les quelques lignes suivantes:
skewness <- function(x)
{
n <- length(x)
1/n * sum((x-mean(x))^3) / ( (n-1)/n*var(x) )^(3/2)
}
Ensuite, pour exécuter la fonction nous n'avons qu'à taper skewness(x) où x est un vecteur de données (par exemple, x <- rnorm(250)). De plus, pour voir le contenu de la fonction "skewness", taper tout simplement son nom sans arguments, ni parenthèses après le prompt de S-Plus (i.e., ">").
Par défaut, la fonction retourne la dernière ligne de commande comme sortie, voir l'exemple ci-haut. Cependant, il est possible que l'on désire retourner plusieurs valeurs ou des sorties différentes selon un certain critère. Pour ce faire l'utilisation de la commande if et d'une liste sont des outils très utiles. Par exemple on peut utiliser la syntaxe suivante:
ma_fonction <- function( arguments, out=T )
{
instructions
if (out) output <- list(valeurs1)
else output <- list(valeurs2)
output
}
Voici un autre exemple combinant le "skewness" (le coefficient d'asymétrie) et le "kurtosis" (le coefficient d'aplatissement) .
skew.kurt <- function(x)
{
n <- length(x)
skew <- 1/n * sum((x-mean(x))^3) / ( (n-1)/n*var(x) )^(3/2)
kurt <- 1/n * sum((x-mean(x))^4) / ( (n-1)/n*var(x) )^2 - 3
list(skewness=skew, kurtosis=kurt)
}
Envoyer une simulation en Batch
Si le temps d'exécution de notre programme de S-Plus est long il peut être utile, voir indispensable, d'envoyer sa simulation en backgroup (i.e., en Batch). Notons qu'il n'est pas nécessaire d'avoir créé ses propres fonctions S-Plus pour envoyer quelque chose en Batch, on peut simplement taper quelques lignes de commandes. Pour utiliser S-Plus en Batch, on doit taper la commande suivante après le prompt de Unix:
tukey% Splus BATCH fichier_entree fichier_sortie
Le fichier d'entrée contient la liste des commandes S-Plus que l'on désire exécuter. Encore une fois, on tape exactement comme si on utilisait S-Plus de façon interactive. Alors que le fichier de sortie est utilisé par S-Plus pour écrire les différentes sorties qu'il afficherait normalement à l'écran. Remarquer que BATCH doit être écrit en majuscule.
En ce qui concerne sur quel ordinateur envoyer vos simulations, voici quelques conseils et réglements d'ordre général:
1) Vous pouvez envoyer une simulation sur une station (i.e., belle, andre, neyman, ...) si vous travaillez présentement sur cette station. Cependant, soyez conscient que ces stations n'ont pas la puissance de tukey, en particulier, elles n'ont que 64 Meg de Ram, alors que tukey en a 768 Meg. Par conséquent, selon la taille et la complexité de votre simulation, la station peut être plus lente ou même très lente. Soulignons que si votre processus demande plus de 64 Meg de Ram (utiliser la commande top de Unix), la station sera très très lente par la suite. Comme l'envoi d'un processus en Batch peut ralentir passablement une station (donc l'usagé qui l'utilise), vous comprendrez qu'il est formellement interdit d'envoyer un processus en Batch sur une station, de se débrancher et de s'en aller ailleurs (sans oublier que si vous utilisez startx, ceci causera des problèmes et vous ne pourrez pas vous débranchez). En résumé, utiliser les stations est recommandé pour les petits programmes ou les essais du début.
2) Tukey est capable de prendre des simulations beaucoup plus volumineuses, mais ces ressources à lui non plus ne sont pas illimitées, bien que plus grandes que cramer. Il est maintenant interdit d'envoyer ses simulations sur cramer. Par conséquent, il est conseillé de regarder (en utilisant la commande top de Unix) quelles ressources sont présentement utilisées et, s'il y a déjà 2 ou 3 usagés qui font rouler des simulations (S-Plus ou autre), d'envoyer sa simulation ailleurs ou d'attendre un peu. De toute façon, tout ce qui va arriver c'est que vous allez vous ralentir entre vous, car tukey doit partager son temps de CPU entre les différents processus, ainsi il y a de grosses pertes de temps associées à l'échange du CPU entre les processus.
3) Tout ceci nous amène à parler du serveur maths1. En effet, maths1 est le serveur dédié pour les simulations et le calcul intensif au Département. Notons, en autres, qu'il a 8 CPU. Si vous n'avez pas déjà un compte sur maths1, je vous encourage fortement à aller voir votre co-administratrice (local 4249) pour en obtenir un.
Comme vous le savez tous probablement S-Plus sauvegarde vos différents objets dans un répertoire nommé ".Data" (pour le voir faire la commande ls -a de Unix dans votre répertoire). De plus, S-Plus utilise une partie de la mémoire vive de l'ordinateur pour garder les valeurs de certains objets temporaires. Ceci est appelé, en langage de S-Plus, le "frame 0". Par conséquent, si vous envoyez un simulation en Batch et que vous travaillez avec S-Plus de façon interactive tout ça dans le même répertoire ".Data" et avec le même "frame 0" (les deux vont de paire), ceci créera des conflits et des mélanges. Prenons le cas où vous accédez à un objet de votre répertoire ".Data" en même temps que votre simulation qui roule en arrière plan.
La solution est tout simplement, de se créer d'autres répertoire ".Data". Pour ce faire il faut d'abord se créer un nouveau répertoire que l'on peut nommer simulation, par exemple (taper mkdir simulation après le prompt de Unix). Ensuite, se déplacer de ce répertoire (cd simulation) et finalement se créer un nouveau ".Data" (mkdir .Data). Alors, il ne nous reste plus qu'à lancer la commande Splus BATCH fichier_entree fichier_sortie lorsque l'on se trouve dans le répertoire simulation et le tour est joué.
Ainsi, on peut continuer à travailler dans le ".Data" de notre répertoire principal. Notons, que ces commentaires ne s'appliquent pas seulement, si on travail avec S-Plus de façon interactive alors que nous avons une simulation qui roule en arrière plan, mais également au cas où nous avons plusieurs simulations en Batch.
2. Boucles et instructions de contrôle:
L'utilisation des boucles for est très populaire dans tout les langages de programmation et S-Plus ne fait pas exception. La syntaxe à suivre est la suivante:
for(compteur in valeurs)
{
instructions
}
Comme vous le savez probablement, la boucle for fonctionne de la façon suivante: la variable compteur (notée généralement par "i" ou "j") prend successivement la valeur des éléments du vecteur valeurs et la partie instructions est alors exécutée de façon indépendante pour chaque élément du vecteur valeurs. Malheureusement, les boucles for ne sont pas très efficaces en S-Plus. En effet, si le nombre d'itérations est élevé elles deviennent rapidement lentes et gourmandes en mémoire. Nous aborderons ce problème plus en détail au chapître suivant.
Voici une fonction utilisant une boucle for réalisée dans le cadre de mon cours de Bootstrap. Elle permet d'estimer la variance et le biais si la médiane est utilisée comme estimateur.
bootstrap <- function(x, B)
{
median.star <- rep(0, B)
for(b in 1:B)
{
median.star[b]<- median(sample(x, size=length(x), replace=T))
}
var.boot <- var(median.star)
biais.boot <- mean(median.star) - median(x)
list(variance=var.boot, biais=biais.boot)
}
La sélection avec l'instruction if
Il y a deux syntaxes possible pour cette commande, voici la première:
if (condition) { instructions}
La condition est évaluée, si elle est vraie alors la partie instructions est exécutée, sinon on passe tout simplement à l'instruction suivante du programme.
La seconde syntaxe possible est:
if (condition) { instructions }
else { instruction_secondaires }
Encore une fois, la condition est évaluée en premier, si elle est vraie alors la partie instructions est exécutée, sinon la partie instructions secondaires est exécutée. Comme en Fortran, en Pascal ou en C++, on peut également utiliser des if imbriqués.
Un élément essentiel que l'on retrouve lorsque l'on définit la condition dans notre if est l'opérateur logique. En voici une liste:
|
|
|
|
|
|
|
égale à |
|
et vectoriel |
|
|
différent de |
|
et logique |
|
|
supérieur à |
|
ou vectoriel |
|
|
supérieur ou égal à |
|
ou logique |
|
|
inférieur à |
|
non logique |
|
|
inférieur ou égal à |
Contrairement à plusieurs langages de programmation, la boucle repeat en S-Plus ne contient pas d'instruction until. Il faut donc l'utiliser en combinaison avec l'instruction break, sinon on se retrouve dans une boucle infinie. Cependant, après quelques utilisations on s'aperçoit que cette structure a l'avantage que l'on peut mettre la condition d'arrêt n'importe où dans la liste d'instructions et non pas seulement à la fin. La marche à suivre est la suivante:
repeat
{
bloc_d'instructions#1
if (condition) break
bloc_d'instructions#2
}
On effectue d'abord le bloc d'instructions #1 après, on évalue la condition, si elle est vrai on quitte la boucle, sinon on effectue le bloc d'instructions #2 pour ensuite recommencer l'exécution du bloc #1 et ainsi de suite...
Voici un exemple avec le calcul des pk pour le plan de Poisson en échantillonnage. On se rappelle qu'il faut que pk £ 1" k OE U (ce programme provient de mon projet en Stt6005).
mu284.dat <- read.table("/home/boudreau/mu284.dat")
mu284.dat <- mu284.dat[ , -1]
dimnames(mu284.dat) <- list( 1:284, c("P85", "P75","RMT85",
"CS82", "SS82", "S82", "ME84", "REV84", "REG", "CL"))
n <- 75
U.tp <- NULL
pik <- rep(0, 284)
repeat
{
pik <- (n - length(U.tp)) * mu284.dat[ ,"P75"] /
sum(mu284.dat[pik<1, "P75"])
pik[U.tp] <- 1
if ( sum(pik>1)==0 ) break
pik[pik>1] <- 1
U.tp <- (1:284)[pik>=1]
}
Autres instructions de contrôle
Parmi les autres commandes intéressantes lorsque l'on écrit une fonction en S-Plus, on retrouve:
® while, qui est semblable au repeat
® switch, qui est l'équivalent du case en Pascal et peut remplacer des if imbriqués
® ifelse, qui est une version vectoriel du if standard
® stop, qui permet d'arrêter une fonction S-Plus et de renvoyer un message d'erreur à l'utilisateur
® warning, qui permet de renvoyer un avertissement à l'utilisateur (faire attention à la différence avec la commande warnings qui permet d'afficher les avertissements en questions)
La première chose que l'on constate en utilisant les boucles for c'est que le temps d'exécution n'est pas linéaire. En effet, le temps d'exécution d'une boucle de taille 2000 n'est pas le double de celui d'une boucle de taille 1000. En réalité, le temps d'exécution augmente de façon quadratique, mais le facteur en question est très petit donc les effets ne se font sentir que le lorsque le nombre d'itérations est grand.
Le problème avec la boucle for est au niveau de l'allocation de la mémoire vive (i.e., Ram). En effet, cette allocation s'effectue à chaque itération, en fonction du travail demandé. Ainsi, plus le nombre d'itérations augmente plus le processus associé grossit en espace mémoire et devient lent. De plus, comme S-Plus ne se débarrasse pas très bien de la mémoire qui n'est plus utilisée, on se retrouve à gaspiller inutilement de l'espace mémoire. Pire, cette mémoire inutile n'est pas libérée tant que la boucle for n'est pas complétée.
Il est difficile de dire à partir de combien d'itérations le problème devient sérieux, car ceci dépend de la complexité des calculs et de la qualité de la programmation, mais les livres de S-Plus parle de 10 000 itérations et plus. Notons que les boucles for imbriquées sont extrêmement lentes et qu'il faut les éviter à tout prix.
Les fonctions apply, lapply et sapply
Un bon moyen d'éviter les boucles for est d'avoir recours aux fonctions apply, lapply et sapply. La fonction apply permet "d'appliquer" une certaine fonction sur les lignes ou colonnes d'une matrice. Les fonctions lapply et sapply sont similaires, mais sont destinées aux listes. La syntaxe de la fonction apply est la suivante:
apply(matrice, dimension, fonction, ...)
où matrice est la matrice sur laquelle on veut effectuer les calculs, dimension est remplacé par 1 pour indiquer que les calculs doivent être effectués sur les lignes alors que 2 est utilisé pour les colonnes, fonction est la fonction à utiliser en question et les trois " ..." peuvent être utilisés pour passer certains arguments supplémentaires à la fonction (voir plus loin pour plus d'explications sur les trois " ...").
Voici un exemple avec jeux de donnés "ethanol" que l'on retrouve dans S-Plus:
> apply(ethanol, 2, mean)
NOx C E
1.957375 12.03409 0.9264773
Un second exemple, serait celui de la fonction de Bootstrap de la page 5. En effet, en utilisant la fonction apply, il est possible de ne pas utiliser la boucle for, ce qui augmente beaucoup la rapidité d'exécution.
bootstrap <- function(x, B)
{
x.star <- matrix(sample(x, size=length(x)*B, replace=T), B, length(x))
median.star <- apply(x.star, 1, median)
var.boot <- var(median.star)
biais.boot <- mean(median.star) - median(x)
list(variance=var.boot, biais=biais.boot)
}
Comme le vecteur est l'objet de base en S-Plus, la majorité des fonctions performent leur calcul sur tout le vecteur. Il est donc inutile de faire un boucle for de façon explicite. A vous de faire bon usage de cette possibilité et d'utiliser les vecteurs au maximum. En particulier, notons qu'il est possible d'utiliser les "subscripts" (tant pour les vecteurs que pour les matrices) et, ainsi, réduire le temps d'exécution. Voici un exemple illustrant bien cette dernière possibilité:
for( i in 1:dim(x)[[1]] )
for( j in 1:dim(x)[[2]] )
if( x[i, j] < 0 ) x[i, j] <- 0
En effet, on peut remplacer les deux boucles for imbriquées, par cette simple ligne qui sera exécutée beaucoup plus rapidement par S-Plus.
x[x<0] <- 0
Si malgré tout il faut absolument utilisée une boucle for avec un nombre élevé d'itérations, il est alors beaucoup plus efficace d'utiliser la fonction For (F majuscule) de S-Plus que le traditionnel for. Voici la syntaxe à suivre:
For(i = 1:nombre_d'iterations,
{
bloc d'instructions
}, grain.size=voir_plus_bas)
Le grain.size défini le nombre d'itérations à effectuer, dans chaque processus indépendant de la fonction For. Par exemple, si nombre_d'itérations = 5000 et grain.size = 10, ceci veut dire que les 5000 itérations seraient divisées en 500 processus de 10 itérations, avec chaque processus utilisant sa propre mémoire allouée. De cette manière on obtient une économie sur la mémoire, augmentant la rapidité d'exécution du programme S-Plus. À noter qu'il est possible que l'object.size soit trop petit lorsque l'on utilise ce type de boucle. C'est qu'avec la boucle For, la mémoire est allouée dès le départ du processus S-Plus associé, et non pas au fur et à mesure. Donc, il est possible que la taille maximale d'un objet aie besoin d'être ajustée. Pour ce faire on tape la commande suivante:
options(object.size = taille_en_question)
Pour un exemple complet, voir le programme en annexe.
Cette section contient quelques conseils et recommandations d'ordre général pour bien programmer avec S-Plus. Je vous encourage donc à les mettre en pratique lorsque vous écrirez vos prochaines fonctions.
En effet, pour pouvoir se retrouver dans nos fonctions et savoir pourquoi nous avons écrit telle et telle lignes, il est utile de se laisser quelques commentaires. De plus, si une autre personne lit notre fonction, il saura plus facilement ce que l'on a fait.
Pour mettre un commentaire, il suffit de le précéder par un "#".
# ceci est un commentaire
Il est important et souvent pratique de savoir où notre simulation est rendue. Il est ainsi plus facile d'estimer le temps restant d'exécution et surtout d'attente après nos résultats. Un bon moyen est de rajouter la ligne suivante dans notre boucle for (ou For).
if ( i%%150==0 ) cat(i, ")", date(), "\n")
Il est également recommandé d'utiliser la commande top de Unix pour voir depuis combien de temps notre simulation roule, combien de mémoire vive elle prend et si elle est en mode sleep (ce qui n'est pas bon signe) ou cpu.
Pour être plus précis, ce qu'il faut éviter est d'avoir des objets S-Plus dont la taille augmente. Comme les fonctions cbind et rbind sont souvent associées à ce problème nous avons intitulé la section "Éviter les cbind et rbind". En effet, ceci oblige S-Plus a ré-allouer sa mémoire à chaque fois, ce qui gaspille temps et mémoire vive. Voici donc, un exemple à éviter:
grow <- function( )
{
x <- NULL
for ( i in 1:100)
x <- rbind(x, i:(i+9))
}
Réutiliser les calculs et les mêmes noms
Il arrive fréquemment que l'on doive conserver un résultat temporaire dans une variable. Cependant, l'exemple suivant illustre bien qu'il est avantageux d'utiliser le même nom de variable, autant que possible. En effet, la fonction g demande 6 millions de bytes de mémoire, alors que la fonction g1 (qui n'est q'une très légère modification de la première) n'en demande que 4.5 millions.
g <- function(n)
{
tmp <- runif(n)
tmp1 <- 2 * tmp + tmp^2
tmp2 <- tmp1 - trunc(tmp1)
mean(tmp2 > 0.5)
}
g1 <- function(n)
{
tmp <- runif(n)
tmp <- 2 * tmp + tmp^2
tmp <- tmp - trunc(tmp)
mean(tmp > 0.5)
}
Il est possible de donner des valeurs par défaut à un (des) argument(s) de notre fonction S-Plus. Pour ce faire, il suffit tout simplement d'utiliser le "=". Par exemple, regardons le cas de la fonction mean dont voici l'entête:
mean <- function(x, trim = 0, na.rm = F)
{
le reste de la fonction
}
Ainsi, on constate que la valeur par défaut de l'argument trim (le pourcentage de troncature) est 0, donc on travail avec la moyenne traditionnelle et non pas la moyenne tronquée. La valeur par défaut du 3e argument, na.rm, est F (pour False). Donc, si nous avons des valeurs manquantes parmi les données, la moyenne sera elle aussi une valeur manquante. Si on change la valeur pour T, alors on obtient la moyenne sur les données restantes.
Une autre possibilité intéressante de S-Plus est l'utilisation des trois ..., ceci remplace un nombre indéterminé d'arguments. Par exemple, si on veut une fonction Bootstrap qui permet d'utiliser différent estimateur.
bootstrap <- function(x, B, estimateur, ...)
{
x.star <- matrix(sample(x, size=length(x)*B, replace=T), B, length(x))
s.star <- apply(x.star, 1, estimateur, ...)
var.boot <- var(s.star)
biais.boot <- mean(s.star) - estimateur(x, ...)
list(variance=var.boot, biais=biais.boot)
}
options(object.size= 50920016)
nb.echantillon <- 75
nb.monte.carlo <- 5000
alpha <- 0.05
moyenne.pop <- mean(mu284.dat[ ,"REV84.corr"])
#initialisations#
s.var.dilat <- NULL
s.mean.dilat <- NULL
TRMC.dilat <- 0
s.var.reg <- NULL
s.mean.reg <- NULL
TRMC.reg <- 0
s.var.reg2 <- NULL #sans les poids gk
TRMC.reg2 <- 0
s.var.ra <- NULL
s.mean.ra <- NULL
TRMC.ra <- 0
For(i=1:nb.monte.carlo,
{
s <- mu284.dat[(sample(1:284, nb.echantillon, replace=F)),
c("REV84.corr", "P75", "RMT85")]
s.y <- s[ ,"REV84.corr"]
s.x <- matrix(unlist(s[ , c("P75", "RMT85")]), nb.echantillon, 2)
#estimateur par dilatation
s.mean.dilat[i] <- mean(s.y)
s.var.dilat[i] <- (1/nb.echantillon - 1/284)*var(s.y)
if (moyenne.pop >= s.mean.dilat[i] - qnorm(1-alpha/2)*sqrt(s.var.dilat[i])
& moyenne.pop <= s.mean.dilat[i] + qnorm(1-alpha/2)*sqrt(s.var.dilat[i]))
{TRMC.dilat <- TRMC.dilat + 1/nb.monte.carlo}
#estimateur par le ratio
R.hat <- mean(s.y) / mean(s.x[ , 1])
s.mean.ra[i] <- mean(mu284.dat[ , "P75"]) * R.hat
s.var.ra[i] <- (1/nb.echantillon - 1/284)*
mean(mu284.dat[ , "P75"])^2 / mean(s.x[ , 1])^2 *
1/(nb.echantillon - 1) * sum( (s.y-R.hat*s.x[ , 1])^2 )
if (moyenne.pop >= s.mean.ra[i] - qnorm(1-alpha/2)*sqrt(s.var.ra[i])
& moyenne.pop <= s.mean.ra[i] + qnorm(1-alpha/2)*sqrt(s.var.ra[i]))
{TRMC.ra <- TRMC.ra + 1/nb.monte.carlo}
#estimateur par regression
T.matrix <- 284/nb.echantillon * t(s.x)%*%s.x
Beta.hat <- as.vector(284/nb.echantillon * solve(T.matrix)%*%t(s.x)%*%s.y)
tempo <- apply(mu284.dat[ ,c("P75", "RMT85")], 2, mean) - apply(s.x, 2, mean)
s.mean.reg[i] <- mean(s.y) + Beta.hat %*% as.vector(tempo)
tempo.gk <- 1 + 284*tempo%*%solve(T.matrix)%*%t(s.x)
s.var.reg[i] <- (1/nb.echantillon - 1/284) *
var( as.vector(tempo.gk*as.vector(s.y - s.x%*%Beta.hat)) )
if (moyenne.pop >= s.mean.reg[i] - qnorm(1-alpha/2)*sqrt(s.var.reg[i])
& moyenne.pop <= s.mean.reg[i] + qnorm(1-alpha/2)*sqrt(s.var.reg[i]))
{TRMC.reg <- TRMC.reg + 1/nb.monte.carlo}
#estimateur par regression sans les poids gk
s.var.reg2[i] <- (1/nb.echantillon - 1/284) *
var( as.vector(s.y - s.x%*%Beta.hat) )
if (moyenne.pop >= s.mean.reg[i] - qnorm(1-alpha/2)*sqrt(s.var.reg2[i])
& moyenne.pop <= s.mean.reg[i] + qnorm(1-alpha/2)*sqrt(s.var.reg2[i]))
{TRMC.reg2 <- TRMC.reg2 + 1/nb.monte.carlo}
if (i%%250==0) {cat(i, ") ", date(), "\n")}
}, grain.size=15)
EMC.moyenne.dilat <- mean(s.mean.dilat)
VMC.dilat <- var(s.mean.dilat)
EMCV.dilat <- mean(s.var.dilat)
EMC.moyenne.reg <- mean(s.mean.reg)
VMC.reg <- var(s.mean.reg)
EMCV.reg <- mean(s.var.reg)
EMCV.reg2 <- mean(s.var.reg2)
EMC.moyenne.ra <- mean(s.mean.ra)
VMC.ra <- var(s.mean.ra)
EMCV.ra <- mean(s.var.ra)
#faire afficher les resultats
moyenne.pop
EMC.moyenne.dilat
EMC.moyenne.reg
EMC.moyenne.ra
TRMC.dilat
TRMC.reg
TRMC.reg2
TRMC.ra
VMC.dilat
EMCV.dilat
VMC.reg
EMCV.reg
EMCV.reg2
VMC.ra
EMCV.ra
#Voici une fonction qui choisit une modele de regression selon la methode
Cp de
#Mallow, Forward ou Backward et qui nous donne le nombre de fois que
le modele
#contient 1, 2, ... n variables.
#La première colonne de X est la variable dépendante(y).
Les suivantes sonts
#les variables indépendantes
ChoixBootstrap <- function(x, methode = "Cp", B){
#fonction ayant pour argument une matrice x, une methode de selection
(Cp,
#forward ou backward), un nombre de repetition B.
#Fixons la taille a NULL (i.e. taille de dimension 0)
#On cree les objets nb.row.x et nb.col.x contenant les dimensions de
x
taille_NULL
nb.row.x_nrow(x)
nb.col.x_ncol(x)
#De 1 a B
for(i in 1:B){
#On imprime la date a toutes les 50 iterations
if (i%%50==0) {cat(i, ")", date(), "\n")}
#On choisit un echantillon avec remise de taille n (nb d'observations)
etiquettes_sample(1:nb.row.x, size=nb.row.x, replace=T)
#On selectionne les lignes de x correspondant a cet echantillon
x.star_x[etiquettes,]
#On calcule selon la methode
if (methode=="Cp") {
choix_leaps(x.star[,2:ncol(x.star)],x.star[,1], method="Cp")
#On cree l'objet cp contenant le cp des modeles de taille 1 a n
cp_choix$Cp
#On choisit la taille associe au cp minimal
taille[i]_choix$size[cp==min(cp)]
}
else if (methode=="forward"){
choix_stepwise(x.star[,2:ncol(x)],x.star[,1], method="forward")
#On cree l'objet fstat contenant les statistiques f pour chacune des
tailles
#de modele
fstat_choix$f.stat
#On choisit la taille associee a la premiere stat F plus petite que
2 (Fin)
taille[i]_choix$size[fstat<2][1]
}
else if (methode=="backward"){
choix_stepwise(x.star[,2:ncol(x)],x.star[,1], method="backward")
#On cree l'objet fstat contenant les statistiques f pour chacune des
tailles
#de modele
fstat_choix$f.stat
#On choisit la taille associee a la premiere stat F plus grande que
2 (Fout)
taille[i]_choix$size[fstat>2][1]
}
#Si la methode ne correspond a aucune des 3 connues, on arrete
else stop("Methode inconnue, choisissez Cp, forward ou backward")
}
#On calcule le nombre de fois que le modele contient 1, 2, ..., n variables
table(taille)
}
N.B. Certaines conditions s'appliquent sur la matrice X lorsque
l'on utilise la fonction leaps