POL6022 - Méthodes quantitatives avancées

Théories de la probabilité

Dans ce laboratoire, nous (1) créons des espaces échantillonnaux et des sous-ensembles d’événements, (2) calculons des probabilités additives pour des événements indépendants et non indépendants et (3) combinons et permutons des observations avec et sans remise.

Fonctions pertinentes: sample(), seq(), rep(), combn(), permn()

1. Espaces échantillonnaux et événements

1.1 Créer un espace échantillonnal

Un espace échantillonnal (également représenté par la lettre grecque oméga) englobe tous les événements qui…

omega_piece  <-c("Face","Pile") # Création d'un vecteur appelé omega_piece

omega_piece # Affichage des éléments de notre vecteur omega_piece

## [1] "Face" "Pile"

# Création d'un vecteur appelé omega_de

# Méthode 1: Saisie manuelle des valeurs
omega_de <-c(1,2,3,4,5,6) #

# Méthode 2: Utilisation de la fonction seq()
omega_de <- seq(1,6) 

 # Affichage des éléments de notre vecteur omega_de
omega_de

## [1] 1 2 3 4 5 6

1.2 Extraction d’un espace événementiel

Nous voulons maintenant extraire des espaces événementiels (représentés par des lettres romaines majuscules, comme A) de ces espaces échantillonnaux. Ce sont les événements de l’espace échantillonnal qui se sont effectivement produits.

Nous pouvons utiliser la fonction sample() pour échantillonner à partir d’un vecteur, en spécifiant la probabilité de chaque résultat. Pour ces exemples, nous utilisons un échantillonnage avec remise car chaque essai est indépendant des autres. Concrètement, obtenir Face lors d’un premier pile ou face ne retire pas Face de l’espace échantillonnal pour le deuxième lancer.

# Définition du "seed" pour la reproductibilité
set.seed(123)

# Échantillonnage d'une observation de notre vecteur omega_piece (avec remise)
A <- sample(omega_piece , size=1, replace=TRUE, prob=c(0.5,0.5))

# Affichage de nos résultats
A

## [1] "Pile"

# Échantillonnage de 5 observations de notre vecteur omega_piece (avec remise)
B <- sample(omega_piece , size=5, replace=TRUE, prob=c(0.5,0.5))

# Affichage de nos résultats
B

## [1] "Face" "Pile" "Face" "Face" "Pile"

# Définition du "seed" pour la reproductibilité
set.seed(123)

# Affichage de notre vecteur omega_piece pour voir l'ordre des arguments
omega_piece

## [1] "Face" "Pile"

# Échantillonnage de 5 observations de notre vecteur omega_piece (avec remise)
C <- sample(omega_piece , size=5, replace=TRUE, prob=c(0.9,0.1))

# Affichage de nos résultats
C

## [1] "Face" "Face" "Face" "Face" "Pile"

# Définition de la "seed" pour la reproductibilité
set.seed(123)

# La fonction rep() permet de répéter un élément donné plusieurs fois
rep(0.16,6)

## [1] 0.16 0.16 0.16 0.16 0.16 0.16

# Échantillonnage d'une observation de notre vecteur omega_de (avec remise)
A <- sample(omega_de , size=1, replace=TRUE, prob=rep(0.16,6))

# Affichage de nos résultats
A

## [1] 3

# Échantillonnage de 5 observations de notre vecteur omega_de (avec remise)
B <- sample(omega_de , size=5, replace=TRUE, prob=rep(0.16,6))

# Affichage de nos résultats
B

## [1] 6 4 1 1 2

# N'oubliez pas de définir votre "seed" à 123!
set.seed(123)

# Réponse

set.seed(123)

sample(omega_de, size=5, replace=TRUE, prob=c(rep(0.12,5),0.4))

## [1] 6 2 3 1 1

2. Combinaisons sans remise

Nous pourrions vouloir mesurer toutes les façons dont un ensemble d’éléments peut être distribué, si l’ordre n’a pas d’importance. Disons que les éléments qui nous intéressent sont des lettres de l’alphabet.

Nous avons vu en cours que la formule pour les combinaisons sans remise était la suivante :

\[ \left({{n}\atop k}\right) = \frac{n!}{k! * (n-k)!} \]

Où n est le nombre d’objets que nous voulons distribuer et k est la taille du groupe.

Si nous nous concentrons uniquement sur les lettres de l’alphabet a, b et c, de combien de façons pouvons-nous arranger des dyades (c’est-à-dire des paires de deux) à partir de ces éléments si leur ordre n’a pas d’importance ?

\[ \left({{3}\atop 2}\right) = \frac{3!}{2! * (3-2)!} = \frac{3 * 2 * 1}{(2 * 1) * 1} = \frac{6}{2} = 3 \]

Voici comment nous ferions cela en R (et visualiserions toutes les combinaisons) :

# Génération de toutes les combinaisons possibles pour les trois premières lettres de 
# l'alphabet (paires de 2 éléments)
combn(letters[1:3], 2)

##      [,1] [,2] [,3]
## [1,] "a"  "a"  "b" 
## [2,] "b"  "c"  "c"

# Affichage du nombre de combinaisons (c'est-à-dire le nombre de colonnes dans le tableau)
ncol(combn(letters[1:3], 2))

## [1] 3

Nous pouvons en déduire que vous auriez \(P(a,b) = \frac{1}{3} = 0,\bar{3}\), puisqu’il y a trois façons possibles d’arranger ces lettres.

ncol(combn(1:10, 3))

## [1] 120

3. Permutations

Nous pourrions vouloir mesurer toutes les façons dont un ensemble d’éléments peut être distribué, si l’ordre compte. Disons que les éléments qui nous intéressent sont les jours de l’année (n = 365).

3.1 Avec remise

Nous avons vu en cours que la formule pour les permutations avec remise était la suivante :

\[ _nP_r = n^r \]

Où n est le nombre d’objets que nous voulons distribuer et r est la taille du groupe. Imaginons que je veuille connaître toutes les permutations possibles de dates de naissance entre deux personnes (avec remise, puisque les deux individus peuvent être nés le même jour).

\[ _{365}P_2 = 365^2 = 133 225 \]

Voici comment calculer cela en R, ainsi que voir toutes les permutations possibles :

# Chargement de la bibliothèque gtools
library(gtools)

# Création d'un vecteur avec tous les jours de l'année 
jours <- format(seq(as.Date("2022-01-01"), as.Date("2022-12-31"),
                   by="+1 day"), format="%d-%m")

# Affichage des 5 premières permutations
head(permutations(n=365,r=2,v=jours,repeats.allowed=TRUE),5)

##      [,1]    [,2]   
## [1,] "01-01" "01-01"
## [2,] "01-01" "01-02"
## [3,] "01-01" "01-03"
## [4,] "01-01" "01-04"
## [5,] "01-01" "01-05"

# Affichage du nombre de permutations (c'est-à-dire le nombre de lignes dans le tableau)
nrow(permutations(n=365,r=2,v=jours,repeats.allowed=TRUE))

## [1] 133225

3.2 Sans remise

Introduisons maintenant un nouveau concept : celui des permutations sans remise. C’est un cas où l’ordre compte, et où vous ne pouvez pas assigner le même élément deux fois. Une fois qu’un élément est assigné dans le groupe, il n’est plus disponible pour être assigné à nouveau.

\[ _nP_r = \frac{n!}{(n-r)!} \]

Si nous reprenons le même exemple qu’avant, avec toutes les permutations possibles de dates de naissance pour deux individus, mais cette fois ils ne peuvent pas partager le même anniversaire (permutation sans remise), voici à quoi ressemble la formule :

\[ _{365}P_2 = \frac{365!}{(365-2)!} = \frac{365!}{363!} = 365 * 364 = 132 860 \]

Voici comment faire cela en R :

# Installation et chargement de la bibliothèque combinat

# install.packages("combinat")
library(combinat)

# Création d'un vecteur avec tous les jours de l'année 
jours <- format(seq(as.Date("2022-01-01"), as.Date("2022-12-31"),
                   by="+1 day"), format="%d-%m")

# Affichage du nombre de permutations (c'est-à-dire le nombre de lignes dans le tableau)
nrow(permutations(n=365,r=2,v=jours,repeats.allowed=FALSE))

## [1] 132860

Une dernière remarque : lorsque le nombre d’objets est égal à la taille du groupe, la formule peut être simplifiée en :

\[ \text{si n = r, alors } _nP_r = n! \]

Par exemple :

\[ _3P_3 = \frac{3!}{(3-3)!} = \frac{3!}{0!} = \frac{3!}{1} = 3!\]

Rappel : la factorielle de 0 est 1.

Application : Le problème des anniversaires

Quelles sont les chances que deux personnes dans une classe donnée aient le même anniversaire ? Ou plutôt, puisque c’est plus facile à mesurer, quelle est la probabilité que personne dans une classe donnée ne partage le même anniversaire ?

\[ P(\text{personne ne partage d'anniversaire}) = \frac{\text{Permutations de k anniversaires sans remise}}{\text{Permutations de k anniversaires avec remise}}\]

La formule ressemble à ceci :

\[ P(\text{personne ne partage d'anniversaire}) = \frac{365!}{365^k * (365-k)!}\]

Comme exemple, la probabilité que personne ne partage un anniversaire dans une classe de 3 personnes est :

\[ P(\text{personne ne partage d'anniversaire}) = \frac{365!}{365^3 * (365-3)!}\]

\[ P(\text{personne ne partage d'anniversaire}) = \frac{365 * 364 * 363 * 362 * 361 *\text{ ... }* 2 * 1}{365^3 * (362 * 361 * 360 *\text{ ... }* 2 * 1)!}\]

\[ P(\text{personne ne partage d'anniversaire}) = \frac{365 * 364 * 363 * \cancel{362} * \cancel{361} *\text{ ... }* \cancel{2} * \cancel{1}}{365^3 * (\cancel{362} * \cancel{361} * \cancel{360} *\text{ ... }* \cancel{2} * \cancel{1})!}\]

\[ P(\text{personne ne partage d'anniversaire}) = \frac{365 * 364 * 363}{365^3}\]

\[ P(\text{personne ne partage d'anniversaire}) = 0,9917958\]

Par conséquent, la probabilité qu’au moins deux personnes partagent un anniversaire dans un groupe de 3 personnes est :

\[ P(\text{au moins deux personnes partagent un anniversaire}) = 1 - 0,9917958 = 0,0082042 \]

# Fonction birthday() de Imai
birthday <-function(k) {logdenom <- k*log(365) + lfactorial(365 - k)

lognumer <-lfactorial(365)

pr <- 1 - exp(lognumer - logdenom)

return(pr)

}

# Test du résultat pour un groupe de 3 personnes
birthday(3)

## [1] 0.008204166

# Indice : un exemple général sur l'utilisation de la fonction which() :
which(1:10 > 3) 

## [1]  4  5  6  7  8  9 10

# Ce code affiche toutes les observations dans notre vecteur 1:10 qui sont
# supérieures à 3

# Réponse

which(birthday(1:100) > .5)[1]

## [1] 23