Programmation R

Cours d'Introduction à la Programmation en R :

Cette page web ne fait que repprendre les principaux éléments du cours de M1.

Manipulations élémentaires

Lancement de R :

  • Sous Linux : dans un terminal tapez R (puis entrée). Utilisez ensuite l'édideur de votre choix pour écrire les scripts. La transmisson a R, se fera soit par copier-coller soit en utilisant la commande source('votre-fichier.R')dans R. Pour les utilisateur d'Emacs, il existe un package automatisant tout cela : ESS.
  • Sous Windows : Utilisez Tinn-R .Dans le menu 'R' faites 'Start R gui' et utilisez les commandes de Tinn-R pour transmettre vos scripts à R. Il est possible aussi d'utiliser directement R et d'appeller un éditeur intégré par la commande edit('nom-de-fichier.R').

Une grosse calculatrice

Les opérations classiques * / + - sont disponibles. Il convient de noter la présence de la division entière %/% et du modulo %%. Les exemples suivant sont a taper directement dans R. Il n'y a pas a mettre le symbole > qui est juste une invitation de R à entrer une nouvelle commande. Si une ligne comporte un #, tous les caractères situé après sont ignorés par R (il s'agit d'un commentaire).

Un aide mémoire est disponible sur http://cran.r-project.org/doc/contrib/Short-refcard.pdf, fournissant nombre de commandes utiles.

Un excellent site plein de recettes toutes prêtes pour R (in english) : http://www.statmethods.net/index.html.

> 3 + 2
[1] 5
> 3*5
[1] 15
> 50%%30
[1] 20
> 51%/%2
[1] 25

Le [1] qui apparait devant chaque résultat indique que la réponse est un vecteur dont le premier élément est donné après le [1]. De manière générale vous pouvez faire comme si ce chiffre n'était pas la. Il est possible de stocker un résulat dans une variable :

> maVariable <- 10
> maVariable
[1] 10
> maVariable * 2
[1] 20
> maVariable * maVariable
[1] 100
> maVariable <- maVariable + 1 # Attention point important à saisir
> maVariable
[1] 11

Vecteurs

R est en conçu pour travailler facilement avec des vecteurs (des sortes de tableaux de nombres). Observez les résulats suivants :

> 1:10
 [1]  1  2  3  4  5  6  7  8  9 10
> a <- 1:10 # On met dans a les nombres de 1 à 10
> a + a
 [1]  2  4  6  8 10 12 14 16 18 20
> a * 10 
 [1]  10  20  30  40  50  60  70  80  90 100
> b <- 5:8 # les nombres de 5 a 8 dans b
> b # affichage de b pour vérifier
[1] 5 6 7 8
> b[1] # le premier élément de b
[1] 5
> b[4] # le 4e élément de b
[1] 8
> b[5] # le 5e (il n'existe pas, b n'a que 4 éléments)
[1] NA
> b[7]<-9 # On demande a mettre 9 en 7e position dans b
> b # affichage de b
[1]  5  6  7  8 NA NA  9
> seq(5,30,by=5) # Les nombres de 5 à 30 de 5 en en 5
[1]  5 10 15 20 25 30
> a>5 # Opérateur de comparaison sur un tableau entier
 [1] FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
> a[a>5] # Sélection des éléments supérieurs à 5
[1]  6  7  8  9 10

Fonctions

Il est possible de crééer ses propres fonctions :

f <- function(x) { # f(x) = 3*x^2
    return (3*x^2)
    }

g <- function(x,a) { # g(x,a) = a*x^2
    return (a*x^2)
    }
    
> f(3)
[1] 27
> g(3,2)    
[1] 18    

Il est possible de définir une valeur par défaut à un paramètre :


g <- function(x,a=2) { # g(x,a) = a*x^2
    return (a*x^2)
    }
    
> g(3)
[1] 18
> g(3,2)    
[1] 18    

Structures de contrôle

if ... else

a <- 3
if (a > 5) {
    print('a est plus grand que 5')
} else {
    print ('a est plus petit que 5')    
}

Produit :

[1] "a est plus petit que 5"

Boucle for

for (i in 1:3) {
     print('hello world') 
     }

Va produire :

[1] "hello world"
[1] "hello world"
[1] "hello world"
Exercice :
  • Ecrire une fonction qui prend un nombre x et renvoie 1 si x est divisible par 3. et 0 sinon
    • divisibleTrois <- function(x) {
		if (x%%3 == 0 ) {
			return(1)
		} else {
			return(0)
		}
	}

    En fait par convention 1 vaut TRUE et 0 vaut FALSE. On pourrait donc écrire :

    divisibleTrois <- function(x) {
	return(x%%3==0)
}
  • En utilisant la boucle for, écrire un programme qui calcule la somme des n premiers entiers. Pour tester on prendra n egal à 100. A noter qu'en R on effectue généralement cette opération par sum(1:n) (ou mieux en faisant appel a ses souvenir de maths par n*(n+1)/2). 
    • n <- 100
s <- 0
for (i in 1:n) {
	s <- s + i
	}
print(s)
  • En utilisant la boucle for et un if, faire de même pour les n premiers entiers pairs. Plus simplement on peut faire sum((1:(n/2))*2) ou en étant un peu tordu par sum((1:n)[(1:n%%2)==0]).
    • n <- 100
s <- 0
for (i in 1:n) {
	if ((i %% 2) == 0) {
		s <- s + i
		}
	}
print(s)
  • Ecrire une fonction qui recherche le plus petit élément dans un tableau. (la fonction min() de R)
    • minimum <- function(x) {
	mini <- x[1]
	for (candidat in x[2:length(x)]) {
			if (candidat < mini) {
				mini <- candidat	
				}
		}
	return(mini)	
	}
	
x <- runif(100) # Initialisation de x avec 100 valeurs tirées au hasard dans [0,1] uniformément
minimum(x) == min(x) # Va renvoyer TRUE	indiquant que le résultat correspond avec celui donné par min
  • Ecrire une fonction qui renvoie le second plus petit élément (sans trier tout le tableau)
    • secondMinimum <- function(x) {
	mini <- Inf #  Inf correspond à l'infini en R
	secondMini <- Inf
	for (candidat in x) {
		if (candidat < secondMini) {
			if (candidat < mini) {
				secondMini <- mini
				mini <- candidat
			} else {
				secondMini <- candidat
			} 
		}			
	}
	return(secondMini)	
	}
	
x <- runif(100) # Initialisation de x 
secondMinimum(x) == sort(x)[2] # Doit renvoyer TRUE (comparaison du résultat avec le second élément de x trié)
  • Même chose pour le ième minimum. Pour bien faire il vaut mieux écrire une fonction chargée d'insérer un nombre dans un tableau trié. Ici pour faire plus simple, on fera appel à sort (uniquement sur des tableaux de taille i).
    • iMinimum <- function(x,i) {
	iMini <- (1:i) * Inf 
	for (candidat in x) {
		if (candidat < iMini[i]) {
			iMini[i] <- candidat
			iMini <- sort(iMini)
		}			
	}
	return(iMini[i])	
	}
	
x <- runif(100) # Initialisation de x
i <- 8
iMinimum(x,i) == sort(x)[i] # Doit renvoyer TRUE
  • Mêmes questions mais en recherchant la position dans le tableau initial du minimum (ou du ième minimum). On peut obtenir directement ce genre d'information avec les fonctions which.min (which.max) et order.
  • Ecrire une fonction qui prend un paramètre n et génère l'affichage suivant (pour n = 5). Vous aurez besoin de la fonction cat() et du fait qu'un saut de ligne se traduit par le caractère spécial \n.
    • *
    **
    ***
    ****
    *****
     
    etoiles <- function(n) {
	for (i in 1:n) { 
	 	for (j in 1:i) { 
	 		cat('*') 
	 	}
		cat('\n') 
	}	
}

etoiles(5)
  • Afficher le cercle unité dans une fenêtre graphique (nécessite de lire l'aide sur les commandes plot.new() plot.window() et points() , on peut aussi utiliser simplement plot)
    • x <- numeric (2*pi/.01)
y <- numeric (2*pi/.01)
i <- 1
for (theta in seq (from = 0, to = 2*pi, by = .01)) {
  x [i] <- cos (theta)
  y [i] <- sin (theta)
  i <- i + 1
}
plot.new()
plot.window (xlim = c (-1, 1), ylim = c (-1, 1), asp = 1)
points (x, y)

Boucle while

i <- 0
while (i < 3) {
     print('hello world')
     i <- i + 1 
     }

Va produire :

[1] "hello world"
[1] "hello world"
[1] "hello world"
Exercice :
  • Ecrire une fonction qui recherche la position du premier élément négatif dans un vecteur. Ecrire deux version : une avec un for, l'autre avec un while.
    • #Solution avec for
positionPremierNegatif <- function(x) { 
	for (i in 1:length(x)){
			if (x[i]<0) return(i)
		}
	return('pas de négatif')	
}
		
#Solution avec while		
positionPremierNegatif <- function(x) {
	pos <- 1
	while (x[pos]>0) {
		pos <- pos + 1
		if (pos > length(x)) return('pas de négatif')
	}
	return(pos)
}

#Vérification
x <- runif(10,-1,1)
positionPremierNegatif(x) == which(x<0)[1]
  • En utilisant l'algorithme d'Euclide, écrire une fonction qui calcule le pgcd de deux entiers naturels (en utilisant while) a et b. 
    • pgcd <- function (a,b) {
	while ( a%%b != 0) {
		reste <- a%%b
		a <- b
		b <- reste
		}
	return(b)	
	}
	
>pgcd(50,30)		
[1] 10
    ou plus rigolo récursivement (mais sans utiliser while): 
    pgcdr <- function (a,b) {
	if (b == 0) return(a)
	return(pgcdr(b,a%%b))	
	}
	
>pgcdr(50,30)		
[1] 10
  • Ecrire une fonction qui prend en argument deux entiers naturels x,b et qui tranforme l'écriture de x pour le mettre en base b. On représentera le nombre en écrivant chacun des chiffres dans un tableau.
  • Programmer un crible d'Erathostène. L'utiliser pour afficher tous les nombres premiers inférieurs à 100. La fonction which peut vous être utile mais n'est pas indispensable.
    • cribleErathostene <- function(n) {
	x <- 1:n
	positionsBarres <- rep(FALSE,n)
	for (i in 2:length(x)) {
		if ( (x[i]*2 < n) && (!positionsBarres[i]) ){ 
			for (j in seq(2*x[i],n,by=x[i]) ) {
				positionsBarres[j] <- TRUE	
			}
		} 
	}
	return(x[!positionsBarres])
}
cribleErathostene(100)

La suite…

Les matrices

Il est très fréquent en programmation d'avoir besoin de manipuler des tableaux à plus d'une dimension. En R, les tableaux à deux dimensions s'utilisent très simplement par l'intermédiaire de la classe matrix

> m <- matrix(1,3,5) # matrice de 3 lignes 5 colonnes remplie de 1
> m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    1    1    1    1
[3,]    1    1    1    1    1
> m + m # Les opérations classiques fonctionnent (* N'est PAS la multiplication matricielle)
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    2    2    2    2
[2,]    2    2    2    2    2
[3,]    2    2    2    2    2
> m[3,4] <- 8 # Remplacement de l'élément de la 3ème ligne 4ème colonne
> m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    1    1    1    1
[3,]    1    1    1    8    1
> sum(m)
[1] 22
> max(m)
[1] 8
> m[,2] <- rep(4,3) # Remplacement de la 2e colonne par une série de 4 (vecteur de 4 de longueur 3)
> m
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    4    1    1    1
[2,]    1    4    1    1    1
[3,]    1    4    1    8    1

Les tableaux multidimensionnels

Plus généralement, il est possible de définir des tableaux à plus de 2 dimensions

> m <- array(1:8,c(2,2,2)) # Tableau à 3 dimensions (longueur de 2 pour chaque dimension)
> m
, , 1

     [,1] [,2]
[1,]    1    3
[2,]    2    4

, , 2

     [,1] [,2]
[1,]    5    7
[2,]    6    8

Exercices

  • Ecrire une fonction qui prend en argument une matrice et qui renvoie vecteur dont les composantes sont la somme de chacune des colonnes.
  • Ecrire une fonction qui renvoie les n premières lignes du triangle de Pascal dans une matrice.
  • Ecrire une fonction qui détermine si une matrice représente un carré magique.
  • Ecrire une fonction qui fabrique un carré magique (pour les dimensions impaires).
  • Coder un Tic-Tac-Toe (scan et/ou menu peuvent vous être utiles).
  • Coder un jeu de la vie de Conway.
  • Coder et représenter graphiquement un système dynamique de type prédateur-proie. Chercher son attracteur de Lorentz. On utilisera la dynamique suivante (Lotka-Volterra):
    • Formula : \frac{\mathrm{d}x(t)}{\mathrm{d}t} = x(t)(\alpha - \beta y(t))
    • Formula : \frac{\mathrm{d}y(t)}{\mathrm{d}t} = -y(t)(\gamma - \delta x(t))
    où x(t) est l'effectif des proies; y(t) est l'effectif des prédateurs; t est le temps; et les autres variables sont des paramètres propres aux espèces considérées.
  • Calculer 1000 trajectoires du jeu de la roulette avec une banque de 1000 euros et une mise initiale de 1 euro. On double la mise chaque fois que l'on perd. On vous offre 25 euros pour jouer mais on ne vous autorise a retirer l'argent qu'après avoir joué au moins 30 fois 25 euros. Donner une estimation de l'espérance de gain. Essayez de la maximiser.

Et pour aller plus loin :

Des problèmes de programmation de plus en plus difficiles : http://projecteuler.net/. A noter que pour certains de ces problèmes R n'est pas toujours le meilleur choix de langage de programmation.