[PYTHON] 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} $.
- Déterminez la valeur initiale de la valeur estimée de $ \ theta $ $ \ theta ^ {(0)} $
- 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)} $
\sum_i \sum_{z_i} p(z_i|x_i,\theta^{(t)}) \log p(x_i|z_i,\theta) Pour maximiser\theta Calculer et\theta^{(t+1)} À- 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.
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])
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
\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.
- Pour la variable cachée $ z $, utilisez $ \ theta ^ {(t)} $ pour calculer la distribution postérieure $ p (z | x, \ theta ^ {(t)}) $
\sum_z p(z|x,\theta^{(t)}) \log p(x,z|\theta) Pour maximiser\theta Calculer 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
- CS229 Lecture notes, by Andrew NG
- Computational Statistics in Python, Expectation Maximization (EM) Algorithm
Recommended Posts