Exemple d’une suite de lancers à pile ou face

But : prédire le résultat du prochain lancer

On se donne N=10 observations. Ce sera notre vecteur de données D dont nous devons tenir compte pour notre prédiction :

D <- rbinom(10,1,0.5)
N_1 <- sum(D)
D

##  [1] 0 0 1 0 1 1 1 0 1 0

N_1

## [1] 5

Nous disons de façon arbitraire qu’un résultat “pile” est noté 0 et un résultat “face” est noté 1. Avec cette convention nous poursuivons nos raisonnements.

Likelihood : La vraisemblance \(L(\theta)\) est donnée par la loi binomiale \[L(\theta) = p(D|\theta) \propto \theta^{N_1} \cdot (1-\theta)^{N_0},\] où \(N_1=\sum_{i=1}^N I(x_i=1)\) est le nombre de “face”, et de façon analogue, \(N_0=\sum_{i=1}^N I(x_i=0)\) est le nombre de “pile”, et \(I(\cdot)\) est telle que \(I(VRAI)=1\) et \(I(FAUX)=0\).
Ces deux comptes \(N_0\) et \(N_1\) sont les statistiques suffisantes (the sufficient statistics) du vecteur des données \(D\), dans le sens où elles nous suffisent pour inférer le paramètre \(\theta\).

Prior : La loi a priori du paramètre \(\theta\) doit avoir l’intervalle \([0,1]\) comme support (dans la mesure où le paramètre \(\theta\) signifie une probabilité). En outre cette loi doit avoir la même forme que la vraisemblance, de façon à ce qu’elle en soit le prior conjugué. En entend par là que la densité de la loi de \(\theta\) doit avoir la même forme que celle de la loi binomiale. C’est pourquoi nous choisissons que le prior est de loi béta : \(\theta \sim Beta(\theta,1,1)\).
Ici, le choix des paramètres \(a=1\) et \(b=1\) de la loi générale \(Beta(\theta,a,b)\) correspond à un prior non informatif.

Nous vérifions qu’en effet la loi binomiale et la loi Beta ont des densités de même forme, à savoir \[\theta^v \cdot (1-\theta)^w.\] La différence réside dans le fait que pour la loi binomiale, \(\theta\) est le paramètre et \(v\) et \(w\) sont les données statistiques suffisantes, tandis que pour la loi Beta, c’est l’inverse.

Posterior : Nous rappelons que la densité de la loi \(Beta(\theta|a,b)\) est

\[ \begin{eqnarray*} & & \frac{1}{\mbox{Beta(a,b)}} \cdot \theta^{a-1} \cdot (1-\theta)^{b-1} \\ & \propto \mbox{ } & \theta^{a-1} \cdot (1-\theta)^{b-1} \end{eqnarray*} \]

et que la densité de la loi binomiale \(\mathcal{B}(k|\theta,n)\), pour \(k\leq n\), est

\[ \begin{eqnarray*} & & \binom{n}{k} \cdot \theta^k \cdot (1-\theta)^{n-k} \\ & \propto \mbox{ } & \theta^k \cdot (1-\theta)^{n-k} \end{eqnarray*} \]

En faisant un abus d’écriture pour désigner la densité de la loi Beta comme la loi elle-même : \[p(\theta|D) \propto \mathcal{B}(N_1|\theta,N_0+N_1) \cdot Beta(\theta|a,b) =Beta(\theta|a+N_1,b+N_0)\] Pour cette raison, on parle des paramètres \(a\) et \(b\) comme de pseudo-comptes, i.e. des quantités qui s’ajoutent aux comptes \(N_1\), respectivement \(N_0\).

Dans notre cas et puisque nous avons choisi \(a=b=1\), la loi a posteriori - ledit posterior - du paramètre \(\theta\), qui, rappelons-le, signifie une probabilité, à savoir celle de tirer “face” en un lancer, est donnée par :

\[ p(\theta|D) \propto Beta(\theta|N_1+1,N_0+1). \]

Comme estimateur de ce paramètre, nous considérons le maximum de la densité de sa répartition (the MAP estimator) :

\[ \widehat{\theta} = \mbox{argmax}_\theta \space p(\theta|D) \]

Affirmation : Le mode d’une v.a. de loi \(Beta(a,b)\) est \(\mbox{mode}=\frac{a-1}{a+b-2}\).

Conséquence (où \(N=N_0+N_1\)) :

\[ \widehat{\theta}_{MAP} = \frac{a+N_1-1}{a+b+N-2} = \frac{N_1}{N} \quad\mbox{(la fréquence empirique)} \]

Dans notre cas, \(a=b=1\), le nombre total de lancers est \(N=10\) et le nombre de “face” est :

## N_1 =  5

Par conséquent notre estimation du paramètre \(\theta\) est

cat("theta_hat = ",(1+N_1-1)/(1+1+10-2))

## theta_hat =  0.5

Posterior predictive : Nous avons inféré à partir des données \(D\) la valeur du paramètre inconnu \(\theta\). Nous voulons à présent faire la meilleure prédiction pour le lancer suivant, encore inconnu. Par commodité nous définissons la v.a. \(X:=\) résultat du 11e lancer. La formule de marginalisation nous permet alors d’écrire :

\[ p(X=1|D) = \int_0^1 p(X=1|\theta,D) \cdot p(\theta|D) \]

On reconnaît en effet la forme générale d’une formule de marginalisation

\[ \mu(A) = \int_0^1 \mu(A|B) \cdot \mu(B)\,dB, \]

où la mesure \(\mu(\cdot)\) est définie par la probabilité conditionnelle : \(\mu(\cdot):=p(\cdot|D)\).

Dans notre cas,

\[ \begin{eqnarray*} p(X=1|D) & = & \int_0^1 \theta \cdot Beta(\theta|a+N_1,b+N_0)\,d\theta \\ & = & \mathbb{E}\left[\theta|D\right] \\ & = & \frac{a+N_1}{a+N_1+b+N_0} \\ & = & \frac{a+N_1}{a+b+N} \end{eqnarray*} \]

Numériquement :

## [1] 0.5

Or la v.a. \(X\) qui prédit le résultat du 11e lancer est bien évidemment de loi de Bernoulli. Ce que nous venons de calculer, c’est le paramètre de cette loi. De là, la distribution de \(X\) est donnée par :

\[ \begin{eqnarray*} p(X|D) & = & \mathcal{Bern}(X \,|\, \mathbb{E}\left[\theta|D\right]) \\ & = & \mathcal{Bern}\left(X \,\bigg|\frac{a+N_1}{a+b+N}\right) \\ & = & \mathcal{Bern}\left(X \,\bigg|\frac{1+N_1}{12}\right) \end{eqnarray*} \] Nous pouvons à présent générer un 11e lancer en tenant compte des précédents résultats :

parametre <- (1+N_1)/12
parametre

## [1] 0.5

cat("x = ",rbinom(1,1,parametre))

## x =  0

Représentations graphiques

a<-20
b<-0.5
distrib_gamma <- function(t){
  y <- dgamma(t,a,b)
  return(y)
}
t <- seq(0,1,0.01)
plot(t,distrib_gamma(t),type="l")

a<-10
b<-10
distrib_beta <- function(t){
  y <- dbeta(t,a,b)
  return(y)
}
t <- seq(0,1,0.01)
plot(t,distrib_beta(t),type="l")
lines(t,dnorm((t-0.5),0,0.115),col="red")

Un classificateur naïf bayesien Bernoulli multivarié

On se donne une matrice des observations, composée de deux variables (colonnes) explicatives \(X_1\) et \(X_2\) ainsi que de la variable expliquée \(Y\) qui figure l’association d’un couple d’observations \((x_1,x_2)\) à une classe. On classe les observations (lignes de la matrices) en deux catégories. Les associations connues correspondent chacune à une ligne de la matrice.
But : Pour un nouveau couple d’observations non encore associé à une classe, prédire l’association.

Nous générons notre matrice \(M\) des observations :

Y <- rbinom(10,1,0.5)
X_1 <- rbinom(10,1,0.5)
X_2 <- rbinom(10,1,0.5)
cat("X_1 = ",X_1,"\n")

## X_1 =  1 1 0 0 1 0 1 0 0 1

cat("X_2 = ",X_2,"\n")

## X_2 =  1 1 0 1 1 0 1 0 0 0

cat("Y = ",Y)

## Y =  0 1 1 0 1 0 0 1 0 0

M <- matrix(nrow=10,ncol=3)
M[,1] <- X_1
M[,2] <- X_2
M[,3] <- Y
M

##       [,1] [,2] [,3]
##  [1,]    1    1    0
##  [2,]    1    1    1
##  [3,]    0    0    1
##  [4,]    0    1    0
##  [5,]    1    1    1
##  [6,]    0    0    0
##  [7,]    1    1    0
##  [8,]    0    0    1
##  [9,]    0    0    0
## [10,]    1    0    0

Modèle probabiliste dit naïf bayesien : Conditionnellement à une classe, chaque \(X_i\) est de loi de Bernoulli de paramètre dépendant de la classe, et \(Y\) est de loi de Bernoulli de paramètre \(\phi\).

Fixons les notations : Soit \(\underline{X}:=(X_1,X_2)\) le vecteur aléatoire des variables explicatives. Soit \(\underline{\Theta}\) la matrice \((\theta)_{c,j}\) des paramètres des 4 lois de Bernoulli possibles, suivant la catégorie \(c \in \{0,1\}\) et suivant l’indice de la variable explicative \(j \in \{1,2\}\). Soit \(\underline{\theta}_{\,.\,c} = (\theta_{1,c},\theta_{2,c})\) le vecteur colonne de cette matrice pour l’une des catégories \(c\). A l’aide de ces notations, le modèle s’exprime comme :

\[ p(\underline{X}\,|\,Y=c,\, \underline{\theta}_{\,.\,c}) = \prod_{j=1}^2 p(X_j\,|\ y=c,\ \theta_{j,c}) \]

Le modèle est dit naïf parce qu’il suppose que, conditionnellement à la catégorie \(c\), les variables explicatives (the features) sont indépendantes : c’est ce que dit la formule ci-dessus.

Pour prendre davantage de hauteur conceptuelle, remarquons que ce modèle est une factorisation particulière de la probabilité conjointe (joint probability) des \(X_i\) et \(Y\), i.e. une factorisation de \(p(X_1,X_2,Y)\), factorisation donnée par la règle de la chaîne (chain rule) et qu’on représente bien à l’aide d’un graphique directionnel acyclique (DAG).

Model fitting : On calcule le maximum de la vraisemblance afin d’estimer les coefficients de la matrice des paramètres \(\Theta\).

La probabilité pour une seule ligne de notre matrice des observations \(M\) est donnée par (formule de Bayes) :

\[ p(\underline{X},\, Y \,|\, \Theta,\phi) = p(Y \,|\, \underline{\phi}) \cdot \prod_j p(X_{j} \,|\, Y,\,\underline{\theta}_{\,j\,.})\;, \]

où \(\underline{\theta}_{\,j\,.}\) est le vecteur ligne de la matrice \(\Theta\) des coefficients associés à la loi conditionnelle \(p(\cdot\,|\,Y)\) de la v.a. vectorielle \(\underline{X}\).

Considérons la fonction indicatrice \(I(Y=c)\) définie là où vit notre nouvelle observation, i.e. l’espace des valeurs de la v.a. vectorielle \(\underline{X}\), et qui teste si telle réalisation \((x_1,x_2)\) de \(\underline{X}\) est classée dans la classe \(c\) ou non. En particulier, chaque ligne de la matrice \(M\) des observations peut être un argument de cette fonction indicatrice. A l’aide de celle-ci nous pouvons exprimer autrement la relation ci-dessus :

\[ p(\underline{X},\, Y \,|\, \Theta,\phi) = \prod_c \phi_c^{I(Y=c)} \cdot \prod_j \, \prod_c p(X_j\,|\,\theta_{jc})^{I(Y=c)} \]

Cette expression donne la probabilité conjointe (joint probability) des v.a. \(X_1\), \(X_2\) et \(Y\) en fonction des paramètres. Si nous nous donnons des réalisations concrètes \((x_1,x_2,y)\) de ces trois v.a. et que nous les insérons dans l’expression ci-dessus, nous obtenons une probabilité qui est fonction des paramètres, i.e. une fonction de vraisemblance. Cette fonction de vraisemblance est dans ce cas associée à un seul triplet de réalisations, formant une unique observation. Une telle observation pourrait être une ligne de la matrice \(M\).

Afin d’insérer les valeurs concrètes des \(x_j\) valant 0 ou 1, il faut encore préciser, concernant l’expression \(p(X_j\,|\,\theta_{jc})^{I(Y=c)}\), à laquelle des deux valeurs 0 ou 1 est associé le paramètre \(\theta_{jc}\). Nous décidons donc qu’il est associé à la valeur 1, si bien que

\[ p(X_j\,|\,\theta_{jc})^{I(Y=c)} = \theta_{jc}^{I(x_j=1)} \cdot (1-\theta_{jc})^{I(x_j=0)} \]

i.e. \(p(X_j\,|\,\theta_{jc})^{I(Y=c)}\) vaut \(\theta_{jc}\) si \(x_j\) vaut 1 et \(1-\theta_{jc}\) si \(x_j\) vaut 0. Il s’agit en effet d’une loi de Bernoulli.

A ce stade nous avons exprimé la vraisemblance pour une seule observation ou réalisation \((x_1,x_2,y)\) au moyen de fonctions caractéristiques, une pour chaque élément du triplet : les fonctions \(I(x_1=1)\), \(I(x_2=1)\) et \(I(y=c)\), pour \(-\)disons\(-\) \(c=1\).
Remarquons, puisque nos v.a. sont toutes de loi de Bernoulli, que \(I(x_1=0)=1-I(x_1=1)\).

La vraisemblance pour une seule observation \((x_1,x_2,y)\) (i.e. une seule ligne de la matrice \(M\) des observations) est ainsi :

\[ p((x_1,x_2,y)\,|\, \Theta,\phi) = \prod_c \phi_c^{I(y=c)} \cdot \prod_j \, \prod_c \theta_{jc}^{I(x_j=1)} \cdot (1-\theta_{jc})^{I(x_j=0)} \]

La vraisemblance pour toutes nos observations, i.e. pour toutes les lignes \((x_1^i,x_2^i,y^i) \quad i=1,2,\ldots,\mbox{<nombre de lignes>}\) de la matrice \(M\), et non plus seulement pour une seule ligne, est donnée par :

\[ L(\Theta,\phi) = \prod_i p((x_1^i,x_2^i,y^i)\,|\, \Theta,\phi) = \prod_c \phi_c^{I(y^i=c)} \cdot \prod_i \, \prod_j \, \prod_c \theta_{jc}^{I(x_j^i=1)} \cdot (1-\theta_{jc})^{I(x_j^i=0)} \]

Afin de déterminer le maximum de cette vraisemblance (MLE), nous en considérons le logarithme :

\[ \log L(\Theta,\phi) = \sum_i \log \left( p((x_1^i,x_2^i,y^i)\,|\, \Theta,\phi) \right) = \\ \sum_c I(y^i=c) \log(\phi_c) + \sum_i \, \sum_j \, \sum_c I(x_j^i=1) \log(\theta_{jc}) + I(x_j^i=0) \log((1-\theta_{jc})) \] La somme courant sur l’indice i des lignes de la matrice des observations se laisse exprimer à l’aide de comptes, en définissant
\(N_{y=c}:=\sum_i I(y^i=c)\)
ainsi que \(N_{x_j=1}:=\sum_i I(x_j^i=1), \quad j=1,2\).

De là notre log-vraisemblance devient :

\[ \log L(\Theta,\phi) = \sum_c N_{y=c} \log(\phi_c) + \sum_j \, \sum_c N_{x_j=1} \log(\theta_{jc}) + N_{x_j=0} \log((1-\theta_{jc})) \]

Nous observons que cette expression se décompose en une somme de termes, un concernant \(\phi\) et \(2*2*2\) termes concernant les \(\theta_{jc}\). Cette expression se laisse par conséquent optimiser séparément terme par terme. D’autre part nous savons que le maximum de vraisemblance (MLE) pour le modèle Beta-binomial (et d’ailleurs également pour le modèle Dirichlet-multinomial) est la moyenne empirique.

Nous nous rappelons que \(c\in\{0,1\}\)). L’estimateur du maximum de vraisemblance pour le paramètre de la loi de \(Y\) est \(\widehat{\phi}_1 = N_{y=1}/N\), étant bien entendu qu’alors \(\widehat{\phi}_0 = 1-\widehat{\phi}_1\).

De façon analogue, l’estimateur du paramètre \(\theta_{1j}\) est \(\widehat{\theta}_{1j} = N_{x_j=1}/N_{y=1}\), i.e. le nombre de lignes de la matrice des observations associées à la classe \(c=1\) telles que la colonne \(X_j\) prend la valeur 1 divisé par le nombre total des observations associées à la classe \(c=1\).

Quant à l’estimateur du paramètre \(\theta_{0j}\), on a que \(\widehat{\theta}_{0j} = 1-\widehat{\theta}_{1j}\).

Nous sommes à présent en mesure d’implémenter le processus d’apprentissage des paramètres (the model fitting procedure) :

N_classe1 <- vector("numeric",length=2) # C'est le vecteur (N_classe1,1 ; N_classe1,2) des comptes des lignes associées à la classe c=1, pour chacune des variables #explicatives X_j, avec dans notre cas j=1,2.

for(i in 1:10){ # indice i courant sur les lignes de M
  c <- M[i,3] # class label pour la ième ligne de la matrice des observations
  for(j in 1:length(N_classe1)){  # indice j courant sur les variables explicatives X
    if(M[i,j]==1){
      N_classe1[j] <- N_classe1[j]+c      
    }
  }
}
phi <- sum(N_classe1)/10

# Le compte N_xj des lignes où X_j vaut 1, indépendamment de la classe :
N_x1 <- sum(M[,1])
N_x2 <- sum(M[,2])

# Matrice des coefficients :
Theta <- matrix(nrow=2,ncol=2)
# Convention : Nous convenons que la 1ère colonne est pour la classe c=0, la seconde pour la classe c=1, et que la 1ère ligne est pour la variable X_1, la seconde pour X_2.

#pour la classe c=1 (c'est la deuxième colonne de Theta) :
#... et pour X_1 :
Theta[1,2] <- N_classe1[1]/sum(N_classe1)
#... et pour X_2 :
Theta[2,2] <- N_classe1[2]/sum(N_classe1)

#pour la classe c=0 (c'est la première colonne de Theta) :
#... et pour X_1 :
Theta[1,1] <- 1-Theta[1,2]
#... et pour X_2 :
Theta[2,1] <- 1-Theta[2,2]

cat("N_classe1 = ",N_classe1,"\n")

## N_classe1 =  2 2

cat("phi = ",phi,"\n")

## phi =  0.4

print("Theta = ")

## [1] "Theta = "

Theta

##      [,1] [,2]
## [1,]  0.5  0.5
## [2,]  0.5  0.5

La matrice \(\Theta\) est la matrice \((\theta)_{j\,c} \;\), i.e. telle que la première colonne concerne la classe \(c=0\), la seconde la classe \(c=1\), et que la première ligne concerne la variable explicative \(X_1\), première colonne de la matrice des observations \(M\), et la seconde ligne concerne la variable explicative \(X_2\), deuxième colonne de \(M\).

Notre modèle est à présent entraîné (fitted) !

Ce que nous voulons maintenant faire, c’est utiliser notre modèle pour faire une prévision. Qu’entendons-nous par là exactement ? Nous nous donnons une observation non déjà classée, i.e. un couple \(\underline{x}=(x_1,x_2)\), et nous voulons prédire la valeur de \(y\) correspondante. Le but est de calculer :

\[ p(Y=c \,|\, \underline{X}=\underline{x},D) \propto p(Y=c \,|\, D) \,\cdot \prod_j p(X_j=x_j \,|\, y=c,D) \]

Ici, \(D\) figure les données de la matrice des observations \(M\). Quant aux v.a. \(\underline{X}\) et \(Y\), elles se rapportent à la nouvelle observation. Pour comprendre cette relation, il faut comparer avec la relation générale :

\[ P(A|B)= P(A) \cdot P(B|A) \cdot \mbox{constante}, \]

où \(\mbox{constante} = 1/P(B)\). Poser \(A:=\{Y=c\}\), \(B:=\{\underline{X}=\underline{x}\}\) et \(P(\cdot):=p(\cdot|D)\).

Maintenant que nous avons compris cette formule, nous essayons d’en calculer le membre de droite.

C’est alors que nous remarquons que \(p(Y=c \,|\, D)\), i.e. \(P(y=c)\), dépend du paramètre \(\phi\) et que cette dépendance n’apparaît pas dans cette expression, raison pour laquelle nous la faisons apparaître explicitement en remplaçant cette expression selon la formule de marginalisation :

\[ P(y=c) = \int P(y=c \,|\, \phi) \cdot P(\phi)\,d\phi\,,\quad \mbox{où } P(\cdot):=p(\cdot|D). \] On constate que, si on choisit que la loi du prior \(P(\phi)\) est une loi Beta non informative (\(a=b=1\)), cette expression est égale à la moyenne empirique que nous avons déjà calculée \(\widehat{\phi}=N_c/N\).

Quant à l’expression \(p(X_j=x_j \,|\, y=c,D)\), c’est-à-dire \(P(X_j=x_j \,|\, y=c)\) nous procédons de façon analogue, de façon à faire apparaître explicitement la dépendance de \(X_j\) par rapport à son prior \(\theta_{j\,c}\) :

\[ P(X_j \,|\, y=c) = \int Bern((X_j=x_j \,|\, y=c, \, \theta_{j\,c}) \cdot P(\theta_{j\,c}) \; d\theta_{j\,c} \]

Cette expression valant, en choisissant un prior de loi Beta non informative,

\[ (\widehat{\theta}_{j\,c})^{I(X_j=1)} \cdot (1-\widehat{\theta}_{j\,c})^{I(X_j=0)} \]

Par conséquent nous obtenons pour notre posterior predictive :

\[ \begin{eqnarray*} p(Y=c \,|\, \underline{X}=\underline{x},D) & \propto & p(Y=c \,|\, D) \,\cdot \prod_j p(X_j=x_j \,|\, y=c,D) \\ & = & P(Y=c) \,\cdot \prod_j P(X_j=x_j \,|\, y=c) \\ & = & \widehat{\phi} \,\cdot \prod_j (\widehat{\theta}_{j\,c})^{I(X_j=1)} \cdot (1-\widehat{\theta}_{j\,c})^{I(X_j=0)}\,, \end{eqnarray*} \]

où \(\widehat{\theta}_{j\,c}=N_{j\,c}/N_c\quad\) et \(\quad \widehat{\phi}=N_c/N\).

Soit \(v\) la nouvelle observation non encore classée dont nous voulons prédire la classe :

v <- vector("numeric",2)
v <- c(1,0)
v

## [1] 1 0

Nous calculons le posterior predictive, pour \(c=1\) :

posterior_predictive <- phi *   Theta[1,2]^v[1] * (1-Theta[1,2])^(1-v[1]) *
        Theta[2,2]^v[2] * (1-Theta[2,2])^(1-v[2])
posterior_predictive

## [1] 0.1

Finalement et puisque \(Y\) est de loi de Bernoulli, nous produisons notre prédiction en simulant une v.a. de même loi et ayant pour paramètre le posterior predictive que nous venons de calculer :

rbinom(1,1,posterior_predictive)

## [1] 0

Ceci achève l’étude du modèle naïf Bayesien Bernoulli multinomial.
Nous pouvons par la suite porter notre attention sur la possibilité d’avoir une v.a. de loi normale et une v.a. de loi multinomiale, ainsi que 3 classes ou plus.

Le classificateur naïf Bayesien adapté à une matrice des observations \(M\) plus générale

Nous nous donnons une matrice des observations \(M\) structurée de la façon suivante : Chaque ligne correspond à une observation, chaque colonne à part celle tout à droite correspond à une variable explicative \(X_i\), la colonne tout à droite correspondant à la classe \(Y\). Nous fixons les grandeurs qui décrivent notre matrice : soit \(N\) le nombre des observations, i.e. le nombre de lignes de la matrice, soient \(a\) et \(b\) tels que \(a + b = K\) le nombre des variables explicatives, la matrice ayant donc au total \(K+1\) colonnes. Les colonnes sont rangées de gauche à droite de façon à ce que celles occupant de la première position à la \({a}^\mbox{ème}\) position correspondent aux v.a. de lois multinoulli et que les colonnes suivantes (à part bien sûr la dernière) correspondent aux v.a. de lois normales. Il y a donc \(a\) v.a. explicatives de loi multinoulli et \(b\) v.a. explicatives de loi normale. Quant à la v.a. expliquée \(Y\), elle est de loi multinoulli.

Posterior predictive :

Ici, \(D\) figure les données de la matrice des observations \(M\). Quant aux v.a. \(\underline{X}\) et \(Y\), elles se rapportent à la nouvelle observation.

Concernant notre nouvelle observation, nous savons que \(X_j=x_j\;,\quad j=1,\cdots,K\), ce que nous notons vectoriellement par \(\underline{X}=\underline{x}\). Nous voulons déterminer la probabilité que \(Y=c\), pour \(c \in \{1,...,C\}\) l’une des \(C\) classes possibles.

\[ p(Y=c \,|\, \underline{X}=\underline{x},D) \propto p(Y=c \,|\, D) \,\cdot \prod_j p(X_j=x_j \,|\, y=c,D) \]

On pose que \(P(\cdot):=p(\cdot|D)\). Nous avons alors que :

\[ p(Y=c \,|\, D) = P(y=c) = \int P(y=c \,|\, \phi_c) \cdot P(\phi_c)\,d\phi_c \,. \]

\(Y\) étant de loi multinoulli, cette expression est celle de l’espérance du paramètre associé à la classe \(c\), i.e. \(\mathbb{E}[\phi_c]\), parce que \(P(y=c \,|\, \phi_c)=\phi_c\).

Soit \(\underline{\phi}:=(\phi_c)_{c \in C}\) le vecteur des paramètres de la loi multinoulli de \(Y\). Si l’on choisit que la loi du prior \(P(\underline{\phi})\) soit une loi de Dirichlet non informative (tous les \(\alpha_k=1\)), cette expression est égale à la moyenne empirique \(\widehat{\phi_c}=N_c/N\).

Quant au produit \(\prod_j p(X_j=x_j \,|\, y=c,D)\), il convient de distinguer entre les \(a\) v.a. de loi multinoulli, que nous notons par l’exposant “A”, et les \(b\) v.a. de loi normale que nous notons par l’exposant “B”. C’est pourquoi nous écrivons :

\[ \prod_j p(X_j=x_j \,|\, y=c,D) = \prod_{j=1}^a p(X^A_j=x^A_j \,|\, y=c,D) \cdot \prod_{j=a+1}^{a+b} p(X^B_j=x^B_j \,|\, y=c,D)\quad , \] Ici nous avons scindé le vecteur aléatoire selon \(\underline{X}=(\underline{X}^A,\underline{X}^B)\) ainsi que le vecteur de la dernière observation \(\underline{x}=(\underline{x}^A,\underline{x}^B)\).

Portons à présent notre attention sur l’expression \(p(X^A_j=x^A_j \,|\, y=c,D)\).

Puisque \(X^A_j\) est une v.a. de loi multinoulli, son paramètre - son prior - est un vecteur dont la dimension égale le nombre des valeurs discrètes que \(X^A_j\) peut prendre. Nous le notons \(\underline{\theta}\) mais bien entendu ce paramètre dépend de \(j\) et de la classe \(c\). Pour davantage de lisibilité nous laissons néanmoins de côté ces indices \(\underline{\theta}_{j\,c}\), dans la mesure où et \(c\) et \(j\) sont des constantes dans l’expression ci-dessus. Mais il faut alors garder en tête que \(\underline{\theta}\) est le paramètre de la loi multinoulli que suit la \(j^\mbox{ème}\) v.a. explicative \(X^A_j\) lorsque la classe (i.e. \(Y\)) vaut \(c\). C’est le vecteur \(\underline{\theta}=(\theta_1, \cdots,\theta_{a_j})\), où \(a_j\) est le nombre de catégories de la loi multinoulli que suit \(X^A_j\).

Or l’expression \(p(X^A_j=x^A_j \,|\, y=c,D)\) n’exprime pas la dépendance de \(X^A_j\) par rapport à \(\underline{\theta}\). Ce que nous voulons, c’est exprimer explicitement la dépendance de \(X^A_j\) par rapport au prior \(\underline{\theta}\).

On définit pour cela la mesure de probabilité \(\tilde{P}(\cdot):=p(\cdot\,|\, y=c,D)\). Alors

L’intégrale intérieure est une intégrale multiple prise sur toutes les composantes du vecteur \(\underline{\theta}\) à l’exception de la composante \(\theta_{a_k}\), c’est ce que qu’indique \(d{\underline{\theta}}_{-a_k}\), où le signe “moins” est là qui signifie qu’on a intégré sur toutes les composantes à l’exception de la \(a_k^\mbox{ème}\).

Puisque \(x_j^A\) est une valeur prise dans l’ensemble des valeurs possibles \(\{1,2,\cdots , a_j\}\), nous choisissons que, ci-dessus, \(a_k\) soit tel que \(x_j^A=a_k\). Dans ce cas,

C’est l’espérance, étant donné \(Y=c\), du paramètre associé à la valeur discrète \(x_j^A\in\{1,2,\cdots,a_j\}\) de la v.a. explicative multinoulli \(X_j^A\). Si l’on choisit que la loi du prior \(\underline{\theta}\) soit une loi de Dirichlet non informative, cette expression est égale à la moyenne empirique \(\widehat{\theta}_{a_k}=N_{X^A_j=x^A_k}/N_{y=c}\) ,

i.e. le nombre de lignes de la matrice des observations associées à la classe \(c\) telles que la colonne \(X^A_j\) prend la valeur \(x^A_j \in \{1,2,\cdots,a_j\}\) divisé par le nombre total des observations associées à la classe \(c\).

Nous en concluons que le produit \(\prod_{j=1}^a p(X^A_j=x^A_j \,|\, y=c,D)\) est égal à \(\prod_{j=1}^a \mathbb{E}[\theta_{x^A_j}]\) c’est-à-dire à \(\prod_{j=1}^a \widehat{\theta}_{x^A_j}\).

Portons à présent notre attention sur l’expression \(p(X^B_j=x^B_j \,|\, y=c,D)\).

Puisque \(X^B_j\) est une v.a. de loi normale, son paramètre - son prior - est un vecteur de dimension 2. Nous le notons \(\underline{\theta}\) mais bien entendu ce paramètre dépend de \(j\) et de la classe \(c\). Pour davantage de lisibilité nous laissons néanmoins de côté ces indices \(\underline{\theta}_{j\,c}\), dans la mesure où et \(c\) et \(j\) sont des constantes dans l’expression ci-dessus. Mais il faut alors garder en tête que \(\underline{\theta}\) est le paramètre de la loi normale que suit la \(j^\mbox{ème}\) v.a. explicative \(X^B_j\) lorsque la classe (i.e. \(Y\)) vaut \(c\). C’est le vecteur \(\underline{\theta}=(\mu, \sigma^2)\).

Or l’expression \(p(X^B_j=x^B_j \,|\, y=c,D)\) n’exprime pas la dépendance de \(X^B_j\) par rapport à \(\underline{\theta}\). Ce que nous voulons, c’est exprimer explicitement la dépendance de \(X^B_j\) par rapport au prior \(\underline{\theta}\).

On définit pour cela la mesure de probabilité \(p_c(\cdot):=p(\cdot\,|\, y=c)\). Alors

\[ \begin{eqnarray} p(X^B_j=x^B_j \,|\, y=c,D) &=& p_c(X^B_j=x^B_j\,|D) \\ &=& \int p_c(X^B_j=x^B_j \,|\,\underline{\theta},D) \cdot p_c(\underline{\theta},D) \; d\underline{\theta} \\ &=& \int \int p_c(X^B_j=x^B_j \,|\,\mu,\sigma^2,D) \cdot p_c(\mu,\sigma^2|D) \; d\mu \, d\sigma \\ \end{eqnarray} \]

\(p_c(\mu,\sigma^2|D)\) est la distribution jointe a posteriori (joint posterior). Nous l’écrivons comme le produit d’un prior avec la vraisemblance : \(p_c(\mu,\sigma^2|D) \propto p_c(\mu,\sigma^2) \cdot p_c(D \,|\, \mu,\sigma^2)\). Comment choisir la distribution du prior joint (joint prior) \(p_c(\mu,\sigma^2)\) ?

Nous prenons l’option de construire une distribution jointe qui soit conjuguée (à la vraisemblance).

Nous savons que dans le cas où \(\mu\) est inconnu et \(\sigma^2\) est connue, un possible prior sur \(\mu|\sigma^2\) conjugué à la vraisemblance est de loi normale \(\mathcal{N}(\mu_0,\sigma^2/\kappa_0)\), avec la liberté de choix concernant \(\mu_0\) et \(\kappa_0\).

Nous savons aussi que dans le cas où c’est \(\sigma^2\) qui est inconnue et \(\mu\) connu, un possible prior sur \(\sigma^2\) conjugué à la vraisemblance est de loi scaled inverse-\(\chi^2(\nu,s^2)\). C’est la loi inverse-gamma \(\mbox{Inv-}\Gamma(\nu/2,\,s^2\cdot\nu/2)\), alors que la densité de la loi gamma est, rappelons-le \(\Gamma(\theta|\alpha,\beta):=\frac{\beta^\alpha}{\Gamma(\alpha)}\theta^{\alpha-1}\,e^{-\beta\theta}\), et que celle de la loi inverse-gamma est \(\mbox{Inv-}\Gamma(\theta|\alpha,\beta)=\frac{\beta^\alpha}{\Gamma(\alpha)}\theta^{-\alpha-1}\,e^{-\beta/\theta}\). La loi inverse-gamma est la loi de \(1/\theta\) lorsque \(\theta\) est de loi gamma.

Nous écrivons la densité jointe a priori (joint prior) comme \(p_c(\mu,\sigma^2) \propto p_c(\mu|\sigma^2) \cdot p_c(\sigma^2)\). En choisissant

\(\mu|\sigma^2 \sim \mathcal{N}(\mu_0,\sigma^2/\kappa_0)\)

\(\sigma^2 \sim \mbox{Inv-}\Gamma(\nu_0/2,\,\sigma_0^2\cdot\nu_0/2)\),

on obtient pour la densité jointe a priori :

\[ \begin{eqnarray} p_c(\mu_0,\sigma^2) &\propto& (\sigma^2/\kappa_0)^{-1/2} \exp\left(\frac{(\theta-\mu_0)^2}{2\sigma^2/\kappa_0}\right) \cdot \frac{(\sigma_0^2\cdot\nu_0/2)^{\nu_0/2}}{\Gamma(\nu_0/2)}\theta^{-\nu_0/2-1}\,e^{-(\sigma_0^2\cdot\nu_0/2)/\theta} \\ &\propto& \sigma^{-1}(\sigma^2)^{-(\nu_0/2+1)} \exp\left( \frac{-1}{2\sigma^2}\cdot [\nu_0\sigma_0^2 + \kappa_0(\mu_0-\mu)^2] \right) \end{eqnarray} \]

Afin de pouvoir y faire référence par la suite, nous donnons un nom à cette distribution. Nous l’appelons \(\mathcal{N}-\mbox{Inv-}\chi^2\) et nous disons qu’elle correspond ci-dessus aux paramètres \(\mu_0\) et \(\sigma_0^2/\kappa_0\) associés au prior conditionnel \(p_c(\mu \,|\,\sigma^2)\) ainsi qu’aux paramètres \(\nu_0\) et \(\sigma_0^2\) associés au prior \(p_c(\sigma^2)\).

La densité jointe a priori est donc de loi \(\mathcal{N}-\mbox{Inv-}\chi^2(\mu_0,\sigma_0^2/\kappa_0 \,;\, \nu_0,\sigma_0^2)\).

Nous obtenons la densité jointe a posteriori en multipliant ce résultat par la vraisemblance :

\(p_c(\mu,\sigma^2|D) \propto p_c(\mu,\sigma^2) \cdot p_c(D \,|\, \mu,\sigma^2)\)

La vraisemblance dont il est question ici est :

\[ \begin{eqnarray*} p_c(D \,|\, \mu,\sigma^2) &=& \Pi_{i=1}^N \, p_c(y_i \,|\,\mu,\sigma) \\ &\propto& \Pi_{i=1}^N \, \sigma^{-1} \exp\left( -\frac{1}{2\sigma^2}(y_i-\mu)^2 \right) \\ &=& \sigma^{-1} \exp\left( -\frac{1}{2\sigma^2}(\sum_{i=1}^N y_i-\mu)^2 \right) \end{eqnarray*} \]

La densité jointe a posteriori est donc :

\[ \begin{eqnarray*} p_c(\mu,\sigma^2|D) &\propto& p_c(\mu,\sigma^2) \cdot p_c(D \,|\, \mu,\sigma^2) \\ &\propto& \sigma^{-1}(\sigma^2)^{-(\nu_0/2+1)} \exp\left( -\frac{1}{2\sigma^2}\cdot [\nu_0\sigma_0^2 + \kappa_0(\mu_0-\mu)^2] \right) \cdot \sigma^{-1} \exp\left( -\frac{1}{2\sigma^2}(\sum_{i=1}^N y_i-\mu)^2 \right) \\ &=& \sigma^{-1-\nu_0+2-1} \, \exp\left( -\frac{1}{2\sigma^2}\cdot \left[ \nu_0\sigma_0^2 + \kappa_0(\mu_0-\mu)^2 + (\sum_{i=1}^N y_i-\mu)^2 \right] \right) \end{eqnarray*} \]

Un calcul standard montre que la densité jointe a posteriori (joint posterior) est de loi

\[ \mathcal{N}-\mbox{Inv-}\chi^2(\mu_n,\sigma_n^2/\kappa_n \,;\, \nu_n,\sigma_n^2)\,, \]

avec

\[ \begin{eqnarray*} \mu_n &=& \frac{\kappa_0}{\kappa_0+n}\mu_0 + \frac{n}{\kappa_0+n}\overline{y} \\ \kappa_n &=& \kappa_0+n \\ \nu_n &=& \nu_0+n \\ \nu_n \sigma_n^2 &=& \nu_0 \sigma_0^2 + (n-1) s^2 + \frac{\kappa_0 n}{\kappa_0+n}(\overline{y}-\mu_0)^2. \\ \end{eqnarray*} \] Ici, on a posé que

\[ \begin{eqnarray*} \overline{y} &:=& \frac{1}{N}\sum_{i=1}^N y_i & & \mbox{la moyenne empirique}\\ s^2 &:=& \frac{1}{N-1}\left( \sum_{i=1}^N y_i^2-N\overline{y}^2 \right) & & \mbox{la variance empirique} \end{eqnarray*} \] Les paramètres de la distribution jointe a posteriori combinent l’information contenue dans la distribution a priori avec celles contenue dans les données empiriques. Notamment \(\mu_n\) est une moyenne pondérée de la moyenne à priori et de la moyenne empirique.

Maintenant que nous disposons d’une densité pour le joint posterior \(p_c(\mu,\sigma^2|D)\), nous reprenons nos calculs visant à obtenir le posterior predictive. Mais tout d’abord concluons à propos de l’expression \(p(X^B_j=x^B_j \,|\, y=c,D)\) :

\[ \begin{eqnarray} p(X^B_j=x^B_j \,|\, y=c,D) &=& p_c(X^B_j=x^B_j\,|D) \\ &=& \int p_c(X^B_j=x^B_j \,|\,\underline{\theta},D) \cdot p_c(\underline{\theta},D) \; d\underline{\theta} \\ &=& \int \int p_c(X^B_j=x^B_j \,|\,\mu,\sigma^2,D) \cdot p_c(\mu,\sigma^2|D) \; d\mu \, d\sigma \\ &\propto& \int \int p_c(X^B_j=x^B_j \,|\,\mu,\sigma^2,D) \cdot \mathcal{N}-\mbox{Inv-}\chi^2(\mu_n,\sigma_n^2/\kappa_n \,;\, \nu_n,\sigma_n^2) \; d\mu \, d\sigma \\ &\propto& \int \int p_c(X^B_j=x^B_j \,|\,\mu,\sigma^2,D) \, \cdot \\ & & \sigma^{-1}(\sigma^2)^{-(\nu_0/2+1)} \exp\left( \frac{-1}{2\sigma^2}\cdot [\nu_0\sigma_0^2 + \kappa_0(\mu_0-\mu)^2] \right) \; d\mu \, d\sigma \\ \end{eqnarray} \] Nous pouvons à présent mettre bout à bout nos résultats partiels. Le posterior predictive est ainsi donné par :

\[ \begin{eqnarray*} p(Y=c \,|\, \underline{X}=\underline{x},D) &\propto& p(Y=c \,|\, D) \,\cdot \prod_j p(X_j=x_j \,|\, y=c,D) \\ &\propto& p(Y=c \,|\, D) \,\cdot \prod_{j=1}^a p(X^A_j=x^A_j \,|\, y=c,D) \cdot \prod_{j=a+1}^{a+b} p(X^B_j=x^B_j \,|\, y=c,D) \\ &\propto& N_c/N \,\cdot \prod_{j=1}^a \widehat{\theta}_{x^A_j} \quad\cdot \prod_{j=a+1}^{a+b} \int \int p_c(X^B_j=x^B_j \,|\,\mu,\sigma^2,D) \, \cdot \\ & & \sigma^{-1}(\sigma^2)^{-(\nu_0/2+1)} \exp\left( \frac{-1}{2\sigma^2}\cdot [\nu_0\sigma_0^2 + \kappa_0(\mu_0-\mu)^2] \right) \; d\mu \, d\sigma \end{eqnarray*} \]