Programmation avec R

Programming, Probability, Simulation · lab

Browse all programmation documents

| Date : 14-11-2021.

Programmation avec R

Intervenant : Yousri Henchiri.

Lois de probabilit s et simulation d chantillons

Lobjectif dans ce TP est de se familiariser avec la manipulation des lois de probabilit

classiques ainsi quavec leur simulation sur R et illustration de quelques jeux de hasard.

1 G n rer d chantillons al atoires partir dun vecteur

(cid:70) 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)

(cid:70) fournit un chantillon al atoire de 5 valeurs prises dans le vecteur donn , mais avec re-

mise. (replace = TRUE signiant remettre = VRAI) :

sample(vect,5,replace=TRUE)

(cid:70) fournit un chantillon al atoire de length(vect) valeurs prises dans vect (cest- -dire une

permutation).

sample(vect)

(cid:70) fournit un chantillon al atoire de 5 valeurs parmi les 10 premiers entiers (i.e. parmi

seq(1,10)).

sample(10,5)

(cid:70) 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 largument

prob. La somme des poids na 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 dun 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 lobjet RienNeVaPlus, signiant

quil 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 nissons un objet appel pi ce contenant les deux r sultats possibles lors dun

lancer de pi ce, puis nous lutilisons 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 puisquelle 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 dun d

Nous commen ons par la simulation du lancer dun 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 darriv 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 dun objet il sut dentrer 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, cest- -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) signie simplement que lon ne veut faire quun

seul tirage. Pour faire plusieurs lancers de d sans avoir retaper la commande en entier il

sut dutiliser la touche pour invoquer une instruction ant rieure.

sample(d , 1)

[1] 1

sample(d , 1)

[1] 2

sample(d , 1)

[1] 6

Advertisement

On peut galement faire plusieurs tirages dun 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, cest- -dire un syst me

au comportement compl tement d terministe. Il y a quelque chose dun 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 lon a obtenu exactement le m me tirage. Cest pourquoi on parle de

tirage pseudo-al atoires. Si on change la condition initiale on obtient des tirages di 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 dun tirage pseudo-al atoire ! !

Ce que lon souhaite cest 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 :

(cid:70) La face 1 cent-mille-cent-vingt-neuf (100129) fois.

(cid:70) La face 2 quatre-vingt-dix-neuf-mille-neuf-cent-vingt-quatre (99924) fois.

(cid:70) La face 3 quatre-vingt-dix-neuf-mille-sept-cent-trente-quatre (99734) fois.

(cid:70) La face 4 cent-mille-trois-cent-quatre-vingt-neuf (100389) fois.

(cid:70) La face 5 cent-mille-deux-cent-sept (100207) fois.

(cid:70) 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 quun long discours on peut ga-

lement r sumer les r sultats de notre exp rience avec un diagramme en b tons :

barplot(obs, main = "R sultat de lexp rience", xlab = "Face du d ", ylab ="Nombre

de tirages obtenus")

Nous pouvons donc constater que le d nest 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 :

(cid:70) un pr xe : d pour la densit , p pour la fonction de r partition et q pour la fonction

quantile.

(cid:70) un suxe sp ciant la loi laquelle on sint resse : binom pour la loi binomiale, norm

pour la loi normale, etc.

Les param tres des lois sont pr cis s lint rieur de la fonction :

(cid:70) size et prob pour la loi binomiale ;

(cid:70) mean et sd pour la loi normale.

La simulation de donn es consiste typiquement tirer un chantillon (X1, . . . , Xn) de va-

riables ind pendantes et identiquement distribu es (en abr g i.i.d.) suivant la m me loi quune

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 xe r.

Le tableau suivant r capitule les suxes et les fonctions d, p, q et r des lois les plus utilis es

en statistique :

4.1 G n ration dun chantillon gaussien

(cid:70) retourne un vecteur de 100 valeurs tir es dune distribution gaussienne de moyenne 0 et

d cart-type 1.

Advertisement

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

rnorm(100,0,1)

rnorm(100)

(cid:70) obtenir une distribution de valeurs bimodale.

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

4.2 G n ration dun chantillon suivant une loi uniforme

| Date : 14-11-2021.

suxe

masse ou densit

r partition

quantile

Binomiale

Binomiale n gative

G om trique

Hyperg om trique

Poisson

Uniforme

Exponentielle

Normale

Log-normale

Beta

Gamma

Cauchy

Student

Fisher

Weibull

2

Logistique

binom

nbinom

geom

hyper

pois

unif

exp

norm

lnorm

beta

gamma

cauchy

t

f

weibull

chisq

logis

dbinom

dnbinom

dgeom

dhyper

dpois

dunif

dexp

dnorm

dlnorm

dbeta

dgamma

dcauchy

dt

df

dweibull

dchisq

dlogis

pbinom

pnbinom

pgeom

phyper

ppois

punif

pexp

pnorm

plnorm

pbeta

pgamma

pcauchy

pt

pf

pweibull

pchisq

plogis

qbinom

Advertisement

qnbinom

qgeom

qhyper

qpois

qunif

qexp

qnorm

qlnorm

qbeta

qgamma

qcauchy

qt

qf

pweibull

qchisq

qlogis

g n ration de nombres

al atoires

rbinom

rnbinom

rgeom

rhyper

rpois

runif

rexp

rnorm

rlnorm

rbeta

rgamma

rcauchy

rt

rf

rweibull

rchisq

rlogis

4.2 G n ration dun chantillon suivant une loi uniforme

(cid:70) retourne 100 valeurs issues dune distribution uniforme entre 0 et 5.

runif(100,0,5)

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

(cid:70) 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 dobtenir 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 dobtenir 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 dobtenir 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 (cid:54) x).

4.4 Exercices

4.4.1 Exercice 1

(cid:70) Quelle est la probabilit de d passer strictement 4 pour une loi de Poisson de param tre

2.7 (P(2.7)) ?

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

(cid:70) Quelle est la valeur x telle que P (X (cid:54) x) = 0.975 pour une loi normale (quantile) ?

(cid:70) Quel est le quantile 1 % pour une loi T de Student 5 ddl ?

(cid:70) Donner un chantillon al atoire simple de 10 valeurs dune loi de Poisson de param tre

gale 2.7.

Advertisement

(cid:70) Donner un chantillon al atoire simple de 10 valeurs dune loi normale r duite.

(cid:70) Donner un chantillon al atoire simple de 10 valeurs dune loi du Khi2 2 ddl.

(cid:70) Donner un chantillon al atoire simple de 10 valeurs dune loi binomiale n = 100 et

p = 1/2.

4.4.2 Exercice 2

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

2. 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).

3. 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 sint resse la variable al atoire X

correspondant au nombre de Pile obtenus sur les 10 lancers.

1. Quelle est la loi de X ?

2. D terminer P(X = 7).

3. Repr senter la loi de la variable X (argument type="h" de la fonction plot).

4. Soit Y le nombre de Pile obtenus apr s 10 lancers lorsque la pi ce est truqu e avec

probabilit dobtenir Pile gale p = 0.75. Repr senter sur deux graphiques di rents

mais dans une m me fen tre graphique les lois des variables X et Y .

5. Calculer la valeur de la fonction de r partition de la variable Y au point 7.

6. Repr senter la fonction de r partition de Y (plot avec largument type="s").

7. Lire sur le graphique la m diane de Y .

8. 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)