| Date : 14-11-2021.

Programmation avec R

Intervenant : Yousri Henchiri.

 Lois de probabilités et simulation d’échantillons

L’objectif dans ce TP est de se familiariser avec la manipulation des lois de probabilité classiques ainsi qu’avec leur simulation sur R et illustration de quelques jeux de hasard.

1 Générer d’échantillons aléatoires à partir d’un vecteur

⋆ fournit un échantillon aléatoire de 5 valeurs prises dans le vecteur donné, sans remise (la taille de l’échantillon doit être inférieure à la taille du vecteur).

 vect<-seq(3,12)
 sample(vect,5)

⋆ fournit un échantillon aléatoire de 5 valeurs prises dans le vecteur donné, mais avec remise. (replace = TRUE signifiant remettre = VRAI) :

 sample(vect,5,replace=TRUE)

⋆ fournit un échantillon aléatoire de length(vect) valeurs prises dans vect (c’est-à-dire une permutation).

 sample(vect)

⋆ fournit un échantillon aléatoire de 5 valeurs parmi les 10 premiers entiers (i.e. parmi seq(1,10)).

 sample(10,5)

⋆ fournit un échantillon aléatoire de 3 valeurs prises parmi le vecteur donné, chaque valeur étant choisie avec une probabilité proportionnelle aux valeurs données dans l’argument prob . La somme des poids n’a pas besoin d’être 1.

 sample(seq(1,4),3,prob=c(0.1,0.2,0.3,0.4),replace=TRUE)

| Date : 14-11-2021.

2 Simulation de jeux de hasard

2.1 Tirage d’un numéro à la roulette

Une roulette comporte 37 cases numérotées de 0 à 36 :

 0:36 -> roulette
roulette
[1] 0 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
[28] 27 28 29 30 31 32 33 34 35 36

On fait 20 tirages avec remise et on range le résultat dans l’objet RienNeVaPlus, signifiant qu’il est trop tard pour faire une mise :

 sample(roulette, 20, replace = TRUE) -> RienNeVaPlus
RienNeVaPlus
[1] 27 7 0 27 17 20 11 23 15 7 10 8 26 24 18 30 8 34 5 12

La fonction table() permet de compter le nombre de fois où chaque numéro a été obtenu :

 table(RienNeVaPlus)
RienNeVaPlus
0 5 7 8 10 11 12 15 17 18 20 23 24 26 27 30 34
1 1 2 2 1 1 1 1 1 1 1 1 1 1 2 1 1

2.2 Tirage pile ou face

Nous définissons un objet appelé pièce contenant les deux résultats possibles lors d’un lancer de pièce, puis nous l’utilisons comme une urne dans laquelle nous faisons un tirage :

 c("pile", "face") -> pièce
sample(pièce, 1)
[1] "pile"

On peut faire plusieurs tirages de suite mais il faut alors préciser que ces tirages se font avec remise, sans quoi notre urne sera rapidement vide puisqu’elle ne contient que deux éléments. Ceci se fait en ajoutant replace = TRUE :

 sample(pièce, 20, replace = TRUE)
[1] "face" "pile" "pile" "pile" "face" "face" "face" "pile" "face" "pile" "face"
[12] "face" "pile" "pile" "pile" "face" "face" "pile" "face" "pile"

2.3 Lancer d’un dé

Nous commençons par la simulation du lancer d’un dé à six faces. Pour générer tous les entiers de 1 à 6 nous utilisons le caractère de ponctuation deux points ( :) encadré par les valeurs de départ et d’arrivée :

| Date : 14-11-2021.

 1:6
[1] 1 2 3 4 5 6

Nous rangeons cet ensemble de valeurs dans un objet appelé dé :

 1:6 -> dé

Pour examiner le contenu d’un objet il suffit d’entrer son nom dans la console :

 dé
[1] 1 2 3 4 5 6

Nous utilisons maintenant cet objet dé comme une urne dans laquelle nous faisons un tirage aléatoire, c’est-à-dire au hasard, grâce à la fonction sample() :

 1:6
sample(dé, 1)
[1] 6

Le 1 après la virgule dans sample(dé, 1) signifie simplement que l’on ne veut faire qu’un seul tirage. Pour faire plusieurs lancers de dé sans avoir à retaper la commande en entier il suffit d’utiliser la touche ↑ pour invoquer une instruction antérieure.

 sample(dé, 1)
[1] 1
sample(dé, 1)
[1] 2
sample(dé, 1)
[1] 6

On peut également faire plusieurs tirages d’un coup :

 sample(dé, 3)
[1] 3 6 2

# **3 R est-il un bon simulateur ?**

Nous venons de simuler trois jeux de hasard avec un ordinateur, c’est-à-dire un système au comportement complètement déterministe. Il y a quelque chose d’un peu paradoxal ici. On peut se convaincre que les résultats donnés par R ne doivent rien au hasard en contrôlant la condition initiale du générateur de nombres aléatoires grâce à la fonction set.seed() :

 set.seed(1)
matrix(sample(dé), ncol = 2)
[,1] [,2]
[1,] 2 4

| Date : 14-11-2021.

[2,] 6 1
[3,] 3 5

 set.seed(1)
matrix(sample(dé), ncol = 2)
[,1] [,2]
[1,] 2 4
[2,] 6 1
[3,] 3 5

On constate que l’on a obtenu exactement le même tirage. C’est pourquoi on parle de tirage pseudo-aléatoires. Si on change la condition initiale on obtient des tirages différents :

 set.seed(2)
matrix(sample(dé), ncol = 2)
[,1] [,2]
[1,] 2 1
[2,] 4 6
[3,] 3 5

 set.seed(2)
matrix(sample(dé), ncol = 2)
[,1] [,2]
[1,] 2 1
[2,] 4 6
[3,] 3 5

! ! Propriétés souhaitables d’un tirage pseudo-aléatoire ! !

Ce que l’on souhaite c’est que les tirages pseudo-aléatoires se comportent comme si on avait un tirage réellement aléatoire. Reprenons le cas du lancer de dé. Il est raisonnable de souhaiter que le dé ne soit pas truqué et que toutes les faces aient la même probabilité, 1/6, d’être obtenues. Est-ce le cas ? Faisons une expérience, reproductible en invoquant set.seed(1) au préalable, en lançant le dé six-cent- mille (600000) fois :

 set.seed(1)
table(sample(dé, 600000, replace = TRUE)) -> obs
obs
1 2 3 4 5 6
100129 99924 99734 100389 100207 99617

Dans cette expérience nous avons donc obtenu :

⋆ La face 1 cent-mille-cent-vingt-neuf (100129) fois. ⋆ La face 2 quatre-vingt-dix-neuf-mille-neuf-cent-vingt-quatre (99924) fois. ⋆ La face 3 quatre-vingt-dix-neuf-mille-sept-cent-trente-quatre (99734) fois. ⋆ La face 4 cent-mille-trois-cent-quatre-vingt-neuf (100389) fois. ⋆ La face 5 cent-mille-deux-cent-sept (100207) fois. ⋆ La face 6 quatre-vingt-dix-neuf-mille-six-cent-dix-sept (99617) fois.

| Date : 14-11-2021.

Une représentation graphique valant généralement mieux qu’un long discours on peut également résumer les résultats de notre expérience avec un diagramme en bâtons :

 barplot(obs, main = "Résultat de l’expérience", xlab = "Face du dé", ylab ="Nombre
de tirages obtenus")

Nous pouvons donc constater que le dé n’est pas truqué puisque nous avons obtenu à peu près chaque face cent-mille (100000) fois, soit avec une probabilité voisine de 1/6, puisque nous avons fait six-cent-mille (600000) lancers en tout.

4 Loi, fonction de répartition et quantiles

Les fonctions usuelles pour les lois de probabilité sont composées de deux parties :

⋆ un préfixe : d pour la densité, p pour la fonction de répartition et q pour la fonction quantile. ⋆ un suffixe spécifiant la loi à laquelle on s’intéresse : binom pour la loi binomiale, norm pour la loi normale, etc.

Les paramètres des lois sont précisés à l’intérieur de la fonction :

⋆ size et prob pour la loi binomiale ; ⋆ mean et sd pour la loi normale.

La simulation de données consiste typiquement à tirer un échantillon ( X 1 , . . ., Xn ) de variables indépendantes et identiquement distribuées (en abrégé i.i.d.) suivant la même loi qu’une variable aléatoire X . Les fonctions utilisées sur R pour simuler des données suivant une loi s’écrivent de la même manière que pour la densité, la fonction de répartition ou le quantile, mais cette fois avec le préfixe r .

Le tableau suivant récapitule les suffixes et les fonctions d, p, q et r des lois les plus utilisées en statistique :

4.1 Génération d’un échantillon gaussien

⋆ retourne un vecteur de 100 valeurs tirées d’une distribution gaussienne de moyenne 0 et d’écart-type 1.

 rnorm(100,mean=0,sd=1)
 rnorm(100,0,1)
 rnorm(100)

⋆ obtenir une distribution de valeurs bimodale.

 rnorm(1000,mean=c(0,3),sd=c(0.5,1))

4.2 Génération d’un échantillon suivant une loi uniforme | Date : 14-11-2021.

Col1	suffixe	masse ou densité	répartition	quantile	génération de nombres aléatoires
Binomiale	`binom`\|`dbinom`\|`pbinom`\|`qbinom`\|`rbinom`
Binomiale négative	`nbinom`\|`dnbinom`\|`pnbinom`\|`qnbinom`\|`rnbinom`
Géométrique	`geom`\|`dgeom`\|`pgeom`\|`qgeom`\|`rgeom`
Hypergéométrique	`hyper`\|`dhyper`\|`phyper`\|`qhyper`\|`rhyper`
Poisson	`pois`\|`dpois`\|`ppois`\|`qpois`\|`rpois`
Uniforme	`unif`\|`dunif`\|`punif`\|`qunif`\|`runif`
Exponentielle	`exp`\|`dexp`\|`pexp`\|`qexp`\|`rexp`
Normale	`norm`\|`dnorm`\|`pnorm`\|`qnorm`\|`rnorm`
Log-normale	`lnorm`\|`dlnorm`\|`plnorm`\|`qlnorm`\|`rlnorm`
Beta	`beta`\|`dbeta`\|`pbeta`\|`qbeta`\|`rbeta`
Gamma	`gamma`\|`dgamma`\|`pgamma`\|`qgamma`\|`rgamma`
Cauchy	`cauchy`\|`dcauchy`\|`pcauchy`\|`qcauchy`\|`rcauchy`
Student	`t`\|`dt`\|`pt`\|`qt`\|`rt`
Fisher	`f`\|`df`\|`pf`\|`qf`\|`rf`
Weibull	`weibull`\|`dweibull`\|`pweibull`\|`pweibull`\|`rweibull`
χ2	`chisq`\|`dchisq`\|`pchisq`\|`qchisq`\|`rchisq`
Logistique	`logis`\|`dlogis`\|`plogis`\|`qlogis`\|`rlogis`

4.2 Génération d’un échantillon suivant une loi uniforme

⋆ retourne 100 valeurs issues d’une distribution uniforme entre 0 et 5.

 runif(100,0,5)
 runif(100,min=0,max=5)

⋆ obtenir une distribution de valeurs bimodale.

 runif(100,c(0,1),c(5,7))

4.3 Exemples

Donner toutes les valeurs de la loi binomiale de paramètres n = 10 et p = 1 / 3 ( B (10 , 1 / 3)).

 options(digits = 7) # par défaut
dbinom(0:10, 10, 1/3)
[1] 1.734153e-02 8.670765e-02 1.950922e-01 2.601229e-01 2.276076e-01 1.365645e-01
[7] 5.690190e-02 1.625768e-02 3.048316e-03 3.387018e-04 1.693509e-05

 options(digits = 4)
dbinom(0:10, 10, 1/3)
[1] 1.734e-02 8.671e-02 1.951e-01 2.601e-01 2.276e-01 1.366e-01 5.690e-02 1.626e-02
[9] 3.048e-03 3.387e-04 1.694e-05

Quelle est la probabilité d’obtenir 1 avec une loi binomiale B (10 , 1 / 3) ?

 dbinom(1, 10, 1/3)
[1] 0.08671

4.4 Exercices | Date : 14-11-2021.

Quelle est la probabilité d’obtenir plus de 45 et moins de 55 avec une loi binomiale B (100 , 1 / 2) ?

# Calculer P(45 < X < 55) pour X ∼ B (100 , 1 / 2)

sum(dbinom(46:54, 100, 0.5))
[1] 0.6318

Quelle est la probabilité d’obtenir au plus 1 avec une loi binomiale pour B (10 , 1 / 3) ?

 pbinom(1, 10, 1/3)
[1] 0.104

 dbinom(0, 10,1/3)
[1] 0.01734
dbinom(1, 10,1/3)
[1] 0.08671
dbinom(0, 10,1/3) + dbinom(1, 10,1/3)
[1] 0.104

N.B. : Noter que d donne les valeurs P( X = k ) et que p donne les valeurs P( X ⩽ x ).

4.4 Exercices

4.4.1 Exercice 1

⋆ Quelle est la probabilité de dépasser strictement 4 pour une loi de Poisson de paramètre 2.7 ( P (2 . 7)) ?

⋆ Quelle est la probabilité de dépasser 1.96 pour une loi normale centrée réduite ( N (0 , 1)) ?

⋆ Quelle est la valeur x telle que P (X ⩽ x) = 0.975 pour une loi normale (quantile) ?

⋆ Quel est le quantile 1 % pour une loi T de Student à 5 ddl ?

⋆ Donner un échantillon aléatoire simple de 10 valeurs d’une loi de Poisson de paramètre égale à 2.7.

⋆ Donner un échantillon aléatoire simple de 10 valeurs d’une loi normale réduite.

⋆ Donner un échantillon aléatoire simple de 10 valeurs d’une loi du Khi2 à 2 ddl.

⋆ Donner un échantillon aléatoire simple de 10 valeurs d’une loi binomiale n = 100 et p = 1 / 2.

4.4.2 Exercice 2

Soit la variable aléatoire X qui suit la loi normale standard.

a) Calculer les probabilités suivantes :

i) P( X < 1 . 45).

ii) P( X > 1 . 45).

iii) P( |X| < 2 . 05).

4.4 Exercices | Date : 14-11-2021.

b) Trouver la valeur de u dans les cas suivants :

i) P( X < u ) = 0 . 6.

ii) P( X > u ) = 0 . 6.

Soit la variable aléatoire X qui suit la loi normale de paramètres µ = 30 et σ = 1. Calculer P( X < 1 . 45) et P( X < 31 . 45).
Soit la variable aléatoire X qui suit la loi normale de paramètres µ = 0 et σ = 10. Calculer P( X < 1 . 45) et P( X < 14 . 45).

 #–––––––––––––––––––––––––––––––––––
# Solution Exercise 2
#–––––––––––––––––––––––––––––––––––

 #–––––––––––- question 1a –––––––––––––––––pnorm(1.45)
pnorm(1.45, lower.tail=FALSE)
pnorm(2.05)- pnorm(-2.05)
pnorm(2.05)- pnorm(2.05, lower.tail=FALSE)
#–––––––––––- question 1b –––––––––––––––––qnorm(0.60)
qnorm(0.60, lower.tail=FALSE)
#–––––––––––- question 2 –––––––––––––––––pnorm(1.45, mean=30)
pnorm(31.45, mean=30)
#–––––––––––- question 3 –––––––––––––––––pnorm(1.45, mean=0,sd=10)
pnorm(14.5, mean=0,sd=10)

4.4.3 Exercice 3

On jette dix fois de suite une pièce équilibrée et on s’intéresse à la variable aléatoire X correspondant au nombre de Pile obtenus sur les 10 lancers.

Quelle est la loi de X ?
Déterminer P( X = 7).
Représenter la loi de la variable X (argument type="h" de la fonction plot).
Soit Y le nombre de Pile obtenus après 10 lancers lorsque la pièce est truquée avec probabilité d’obtenir Pile égale à p = 0 . 75. Représenter sur deux graphiques différents mais dans une même fenêtre graphique les lois des variables X et Y .
Calculer la valeur de la fonction de répartition de la variable Y au point 7.
Représenter la fonction de répartition de Y (plot avec l’argument type="s").
Lire sur le graphique la médiane de Y .
Retrouver la valeur de cette médiane ainsi que les valeurs des quartiles.

4.4 Exercices | Date : 14-11-2021.

 #–––––––––––––––––––––––––––––––––––
# Solution Exercise 3
#–––––––––––––––––––––––––––––––––––
# Question 1 : X suit une loi binomiale B(10,0.5)
# Question 2
dbinom(7,size=10,prob=0.5) ou directement dbinom(7,10,0.5)
# Question 3
k=0:10
p=dbinom(k,10,0.5)
plot(k,p,type="h")
# Question 4
k=0:10
px=dbinom(k,10,0.5)
py=dbinom(k,10,0.75)
par(mfrow=c(1,2))
plot(k,px,type="h",ylab="Loi de X")
plot(k,py,type="h",ylab="Loi de Y") # X est stochastiquement inferieure a Y par(mfrow=c(1,1)
# Question 5
pbinom(7,10,0.75)
# Question 6
k=0:10
F=pbinom(k,10,0.75)
plot(k,F,type="s",ylab="cdf de Y")
# Question 7
abline(h=0.5,lty=3) # mediane de Y = 8
# Question 8
qbinom(0.5,10,0.75)
qbinom(c(0.25,0.5,0.75),10,0.75)