Dérivation de l'algorithme EM et exemple de calcul pour le tirage au sort

Un exemple de résolution d'un problème de tirage au sort à l'aide de l'algorithme EM est présenté dans les formules mathématiques et le code Python, puis la dérivation de l'algorithme EM lui-même est montrée.

Exemple: estimation de la probabilité qu'un tirage au sort apparaisse

** Problème **: j'ai deux pièces A et B. Lors du tirage au sort, on sait que A est plus facile à sortir que B. Suite à la répétition de l'essai consistant à effectuer un tirage au sort 20 fois en utilisant l'une ou l'autre des pièces 10 fois, le nombre de fois où la table est apparue était le suivant. Estimez la probabilité que le recto de chaque pièce de A et B apparaisse.

Nombre d'apparitions de la table 15 4 6 13 3 5 4 4 7 2

** Réponse **: La probabilité que chaque table de A et B apparaisse est $ \ theta_A, \ theta_B $, et sur $ N $ essais, le résultat du tirage au sort $ M $ dans l'essai $ i $, la table Soit $ x_i $ le nombre de fois que le est émis, et $ z_i \ in \ {A, B \} $ les pièces utilisées.

L'algorithme suivant, appelé algorithme EM, trouve la valeur estimée de $ \ theta = (\ theta_A, \ theta_B) $ $ \ hat {\ theta} $.

  1. Déterminez la valeur initiale de la valeur estimée de $ \ theta $ $ \ theta ^ {(0)} $
  2. Calculer la distribution postérieure $ p (z_i | x_i, \ theta ^ {(t)}) $ de $ z_i $ dans chaque essai $ i $ en utilisant l'estimation $ \ theta ^ {(t)} $
  3. \sum_i \sum_{z_i} p(z_i|x_i,\theta^{(t)}) \log p(x_i|z_i,\theta)Pour maximiser\thetaCalculer et\theta^{(t+1)}À
  4. Mettez $ \ theta ^ {(t + 1)} $ dans $ \ theta ^ {(t)} $ jusqu'à ce que $ \ theta ^ {(t)} $ converge, revenez à l'étape 1 et répétez le calcul.

La formule spécifique pour chaque étape et le code en Python sont indiqués ci-dessous.

étape 1

Soit $ \ theta_A ^ {(0)} = 0,5, \ theta_B ^ {(0)} = 0,4 $.

import numpy as np
import scipy.stats as st

np.random.seed(0)

# Problem setting
N = 10
M = 20
theta_true = np.array([0.8, 0.2])
z_true = np.random.randint(2, size=N)
x = np.zeros(N, dtype=int)
for i in range(0, N):
    x[i] = st.binom.rvs(M, theta_true[z_true[i]])
print(x)

# Initial estimates
t = 0
t_max = 50
thetas = np.zeros([2, t_max])
thetas[:,t] = np.array([0.5, 0.4])

** Étape 2. **

Étant donné $ \ theta ^ {(t)} $, la distribution postérieure de $ z_i $ est

\begin{eqnarray}
p(z_i|x_i,\theta^{(t)})&=& p(z_i|x_i,\theta^{(t)})\\
 &=&  \frac{p(x_i|z_i, \theta^{(t)}) p(z_i)}{\sum_{z_i} P(x_i|z_i,\theta^{(t)})p(z_i)}\\
 &=&  \frac{p(x_i|z_i, \theta^{(t)})}{\sum_{z_i} P(x_i|z_i,\theta^{(t)})}\\
 &=&  \frac{B(x_i|M,\theta_{z_i}^{(t)})}{\sum_{z_i} B(x_i|M,\theta_{z_i}^{(t)})} \ .
\end{eqnarray}

L'équation de la deuxième ligne utilise le théorème de Bayes, et l'équation de la troisième ligne ne contient pas d'informations sur la pièce de monnaie utilisée, donc la distribution précédente est utilisée.p(z_i=A)=p(z_i=B)Calculé comme. Aussi,B(x|n,p)Représente une distribution binomiale. Distribution binaireB(x|n,p)Est un événementU,VDes événementsUEst la probabilitépL'essai sélectionné parnQuand je suis allé indépendammentxFoisUReprésente la probabilité sélectionnée.

def get_post_z_i(z_i, x_i, M, theta):
    norm_c = 0 # Normalization term in denominator
    for j in range(0,2):
        p_j = st.binom.pmf(x_i, M, theta[j])
        norm_c = norm_c + p_j
        if z_i == j:
            ll = p_j
    return ll / norm_c

** Étape 3 **

En utilisant le résultat du calcul de l'étape 1, calculez le maximum $ \ theta $ par la méthode du gradient.

def neg_expect(theta_next, theta_cur, x):
    N = x.size
    e = 0
    for i in range(0, N): # for trial i
        for j in range(0, 2): # used coin j
            post_z = get_post_z_i(j, x[i], M, theta_cur)
            prob_x = st.binom.logpmf(x[i], M, theta_next[j])
            e = e + post_z * prob_x
    return -e

# Sample calculation
bnds = [(0,0.99), (0,0.99)]
res = minimize(neg_expect, thetas[:,t], args=(thetas[:,t], x),
          bounds=bnds, method='SLSQP', options={'maxiter': 1000})
res.x

** Étape 4 **

Répétez jusqu'à converger.

t = 0
while t < t_max-1:
    if t % 10 == 0:
        print(str(t) + "/" + str(t_max))
    res = minimize(neg_expect, thetas[:,t], args=(thetas[:,t], x),
         bounds=bnds, method='SLSQP', options={'maxiter': 1000})
    thetas[:,t+1] = res.x
    t = t + 1

résultat

import matplotlib.pyplot as plt
import seaborn as sns
# result
plt.plot(np.arange(0,t_max), thetas[0,:])
plt.plot(np.arange(0,t_max), thetas[1,:])
# true value
plt.plot([0, t_max], np.ones(2) * theta_true[0])
plt.plot([0, t_max], np.ones(2) * theta_true[1])
plt.ylim([0, 1])

result.png

Dans ce problème, le nombre de données était petit, nous ne pouvions donc pas l'estimer avec autant de précision.

Dérivation de l'algorithme EM

Soit $ x = \ {x_0, \ cdots, x_N \} $ l'échantillon des variables de probabilité obtenues dans N essais, $ z $ la variable de probabilité non observable et $ \ theta $ le paramètre que vous voulez estimer. La fonction de vraisemblance du journal $ l (\ theta) $ pour $ \ theta $ peut être écrite comme suit: (Dans ce qui suit, nous considérerons $ \ sum_z $ en supposant que $ z $ est une variable discrète selon ce problème, mais si $ z $ est une variable continue, il en va de même si nous considérons l'intégration.)

l(\theta) := \log p(x|\theta) = \log \sum_z p(x,z|\theta) \ . (1)

Dans le calcul réel

l(\theta) = \sum_i \log \sum_{z_i} p(x_i,z_i|\theta)\ ,

Il est plus facile à calculer, mais il sera difficile de voir la formule, nous allons donc procéder ici avec la formule (1). Considérons maintenant une distribution de probabilité arbitraire $ Q (z) $ pour $ z $ et utilisez l'inégalité de Jensen.

l(\theta) = \log \sum_z Q(z) \frac{p(x,z|\theta)}{Q(z)} \geq \sum_z Q(z) \log \frac{p(x,z|\theta)}{Q(z)}\ ,

Est établi. Dans l'équation ci-dessus, l'inégalité détient le nombre égal, en supposant que la constante est $ C $.

\frac{p(x,z|\theta)}{Q(z)} = C \ ,

Seulement quand. $ Q (z) $ qui satisfait cela est, en considérant la constante de standardisation

\begin{eqnarray}
Q(z)&=\frac{p(x,z|\theta) }{\sum_zp(x,z|\theta)} \\
&=\frac{p(x,z|\theta)}{p(x|\theta)}\\
&=p(z|x,\theta)  \ ,
\end{eqnarray}

est. Autrement dit, si vous donnez un certain $ \ theta ^ {(t)} $,

\log \sum_z p(z|x,\theta^{(t)}) \frac{p(x,z|\theta^{(t)})}{p(z|x,\theta^{(t)})} = \sum_z p(z|x,\theta^{(t)}) \log \frac{p(x,z|\theta^{(t)})}{p(z|x,\theta^{(t)})}\ ,

Est vrai, peu importe ce que $ \ theta ^ {(*)} \ neq \ theta ^ {(t)} $ est donné

\log \sum_z p(z|x,\theta^{(t)}) \frac{p(x,z|\theta^{(*)})}{p(z|x,\theta^{(t)})} \geq \sum_z p(z|x,\theta^{(t)}) \log \frac{p(x,z|\theta^{(*)})}{p(z|x,\theta^{(t)})}\ ,  \ \ (2)

Est vrai. Par conséquent, la limite inférieure de l'inégalité, c'est-à-dire\log p(x,z|\theta^{(\*)})/p(z|x,\theta^{(t)})dezMaximisez les attentes concernant(Maximization of Expextation)Aimer faire\theta^{(*)}À\theta^{(t+1)}ça ira;

\begin{eqnarray}
\theta^{(t+1)} &=& \mathrm{arg \ max}_{\theta^{(*)}}  \sum_z p(z|x,\theta^{(t)}) \log \frac{p(x,z|\theta^{(*)})}{p(z|x,\theta^{(t)})} \  \ (3)\\
&=& \mathrm{arg \ max}_{\theta^{(*)}}  \sum_z p(z|x,\theta^{(t)}) \log p(x,z|\theta^{(*)}) \ .  \  \ (4)
\end{eqnarray}

Ensuite, nous pouvons voir que l'inégalité suivante est vraie.

\begin{eqnarray}
l(\theta^{(t+1)})
&\geq &\sum_z p(z|x,\theta^{(t)})\log  \frac{p(x,z|\theta^{(t+1)})}{p(z|x,\theta^{(t)})} \\
&\geq& \sum_z p(z|x,\theta^{(t)})\log  \frac{p(x,z|\theta^{(t)})}{p(z|x,\theta^{(t)})} \\
&=& l(\theta^{(t)}) \ .
\end{eqnarray}

Ici, nous avons utilisé la relation (2) pour la première inégalité et la relation (3) pour la seconde inégalité.

Vous avez maintenant dérivé l'algorithme EM! Pour résumer ce qui précède, étant donné un paramètre $ \ theta ^ {(t)} $, le $ \ theta ^ {(t + 1)} $ obtenu par le calcul en deux étapes suivant est $ . Il s'avère que $ l (\ theta ^ {(t)}) \ leq l (\ theta ^ {(t + 1)}) $ est toujours valable pour la fonction de vraisemblance $ l (\ theta) $ de theta $. C'était.

  1. Pour la variable cachée $ z $, utilisez $ \ theta ^ {(t)} $ pour calculer la distribution postérieure $ p (z | x, \ theta ^ {(t)}) $
  2. \sum_z p(z|x,\theta^{(t)}) \log p(x,z|\theta)Pour maximiser\thetaCalculer et\theta^{(t+1)}À

En d'autres termes, si vous décidez $ \ theta ^ {(0)} $ de manière appropriée et répétez la procédure ci-dessus à plusieurs reprises, vous obtiendrez la valeur estimée localement optimale $ \ hat {\ theta} $. C'est l'algorithme EM.

Matériel de référence

Recommended Posts

Dérivation de l'algorithme EM et exemple de calcul pour le tirage au sort
Préparé pour le calcul de la date et l'automatisation de mon bot
Dérivation et implémentation d'équations de mise à jour pour la décomposition de facteurs tensoriels non négatifs
Explication et implémentation de l'algorithme ESIM
Calcul de la classe auto-fabriquée et de la classe existante
[Calcul scientifique / technique par Python] Dérivation de solutions analytiques pour équations quadratiques et cubiques, formules, sympy
Exemple de code python pour la distribution exponentielle et l'estimation la plus probable (MLE)
Explication et implémentation de l'algorithme Decomposable Attention
EM de distribution gaussienne mixte
Algorithme EM modèle mixte gaussien [apprentissage automatique statistique]
Distribution gaussienne mixte et logsumexp
Dérivation de l'algorithme EM et exemple de calcul pour le tirage au sort
PRML Chapitre 9 Implémentation Python de distribution gaussienne mixte
Estimation de la distribution gaussienne mixte par la méthode de Bayes
(Apprentissage automatique) J'ai essayé de comprendre attentivement l'algorithme EM dans la distribution gaussienne mixte avec l'implémentation.