[PYTHON] Ableitung des EM-Algorithmus und Berechnungsbeispiel für den Münzwurf
Ein Beispiel für die Lösung eines Problems beim Münzwurf mit dem EM-Algorithmus wird in mathematischen Formeln und Python-Code gezeigt. Anschließend wird die Ableitung des EM-Algorithmus selbst gezeigt.
Beispiel: Schätzen der Wahrscheinlichkeit, dass ein Münzwurf auftritt
** Problem **: Ich habe zwei Münzen A und B. Beim Münzwurf ist bekannt, dass A leichter herauszukommen ist als B. Als Ergebnis der Wiederholung des Versuchs, den Münzwurf 20 Mal mit einer der beiden Münzen 10 Mal durchzuführen, war die Häufigkeit, mit der die Tabelle erschien, wie folgt. Schätzen Sie die Wahrscheinlichkeit, dass die Vorderseite jeder Münze von A und B erscheint.
| Häufigkeit, mit der die Tabelle angezeigt wurde | 15 | 4 | 6 | 13 | 3 | 5 | 4 | 4 | 7 | 2 |
|---|
** Antwort **: Die Wahrscheinlichkeit, dass jede Tabelle von A und B erscheint, ist $ \ theta_A, \ theta_B $ und aus $ N $ Versuchen das Ergebnis des Münzwurfs von $ M $ in der $ i $ -Versuchung, der Tabelle Sei $ x_i $ die Häufigkeit, mit der die Ausgabe erfolgt, und $ z_i \ in \ {A, B \} $ die verwendeten Münzen.
Der folgende Algorithmus, der als EM-Algorithmus bezeichnet wird, ermittelt den geschätzten Wert von $ \ theta = (\ theta_A, \ theta_B) $ $ \ hat {\ theta} $.
- Bestimmen Sie den Anfangswert des geschätzten Werts von $ \ theta $ $ \ theta ^ {(0)} $
- Berechnen Sie die hintere Verteilung $ p (z_i | x_i, \ theta ^ {(t)}) $ von $ z_i $ in jedem Versuch $ i $ unter Verwendung der Schätzung $ \ theta ^ {(t)} $
\sum_i \sum_{z_i} p(z_i|x_i,\theta^{(t)}) \log p(x_i|z_i,\theta) Um zu maximieren\theta Berechnen und\theta^{(t+1)} Zu- Setzen Sie $ \ theta ^ {(t + 1)} $ in $ \ theta ^ {(t)} $, bis $ \ theta ^ {(t)} $ konvergiert, kehren Sie zu Schritt 1 zurück und wiederholen Sie die Berechnung.
Die spezifische Formel für jeden Schritt und der Code in Python sind unten aufgeführt.
Schritt 1
Sei $ \ 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])
** Schritt 2. **
Die hintere Verteilung von $ z_i $ bei $ \ theta ^ {(t)} $ ist
\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}
Die Gleichung in der zweiten Zeile verwendet den Satz von Bayes, und die Gleichung in der dritten Zeile enthält keine Informationen darüber, welche Münze verwendet wurde, sodass die vorherige Verteilung verwendet wird.
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
** Schritt 3 **
Berechnen Sie anhand des Berechnungsergebnisses in Schritt 1 das Maximum $ \ theta $ nach der Gradientenmethode.
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
** Schritt 4 **
Wiederholen, bis konvergiert.
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
Ergebnis
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])
Bei diesem Problem war die Anzahl der Daten gering, sodass wir sie nicht so genau abschätzen konnten.
Ableitung des EM-Algorithmus
Sei $ x = \ {x_0, \ cdots, x_N \} $ die Stichprobe der in N Versuchen erhaltenen Wahrscheinlichkeitsvariablen, $ z $ die nicht beobachtbare Wahrscheinlichkeitsvariable und $ \ theta $ der Parameter, den Sie schätzen möchten. Die Protokollwahrscheinlichkeitsfunktion $ l (\ theta) $ für $ \ theta $ kann wie folgt geschrieben werden: (Im Folgenden betrachten wir $ \ sum_z $ unter der Annahme, dass $ z $ gemäß diesem Problem eine diskrete Variable ist. Wenn jedoch $ z $ eine kontinuierliche Variable ist, ist es dasselbe, wenn wir die Integration berücksichtigen.)
l(\theta) := \log p(x|\theta) = \log \sum_z p(x,z|\theta) \ . (1)
In der eigentlichen Berechnung
l(\theta) = \sum_i \log \sum_{z_i} p(x_i,z_i|\theta)\ ,
Es ist einfacher zu berechnen, aber es wird schwierig sein, die Formel zu sehen. Deshalb werden wir hier mit Formel (1) fortfahren. Betrachten Sie nun eine beliebige Wahrscheinlichkeitsverteilung $ Q (z) $ für $ z $ und verwenden Sie Jensens Ungleichung.
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)}\ ,
Ist festgelegt. In der obigen Gleichung hält die Ungleichung die gleiche Zahl, vorausgesetzt, die Konstante ist $ C $.
\frac{p(x,z|\theta)}{Q(z)} = C \ ,
Nur wenn ist. $ Q (z) $, das dies erfüllt, ist unter Berücksichtigung der Standardisierungskonstante
\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}
ist. Das heißt, wenn Sie ein bestimmtes $ \ theta ^ {(t)} $ angeben,
\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)})}\ ,
Ist wahr, egal was $ \ theta ^ {(*)} \ neq \ theta ^ {(t)} $ gegeben ist
\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)
Ist wahr. Daher ist die Untergrenze der Ungleichung, dh
\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}
Dann können wir sehen, dass die folgende Ungleichung gilt.
\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}
Hier haben wir die Beziehung (2) für die erste Ungleichung und die Beziehung (3) für die zweite Ungleichung verwendet.
Sie haben jetzt den EM-Algorithmus abgeleitet! Um das Obige zusammenzufassen, ist bei gegebenem Parameter $ \ theta ^ {(t)} $ das $ \ theta ^ {(t + 1)} $, das durch die folgende zweistufige Berechnung erhalten wird, $ . Es stellt sich heraus, dass $ l (\ theta ^ {(t)}) \ leq l (\ theta ^ {(t + 1)}) $ immer für die Wahrscheinlichkeitsfunktion $ l (\ theta) $ von Theta $ gilt. Es war.
- Verwenden Sie für die versteckte Variable $ z $ $ \ theta ^ {(t)} $, um die hintere Verteilung $ p (z | x, \ theta ^ {(t)}) $ zu berechnen
\sum_z p(z|x,\theta^{(t)}) \log p(x,z|\theta) Um zu maximieren\theta Berechnen und\theta^{(t+1)} Zu
Mit anderen Worten, wenn Sie $ \ theta ^ {(0)} $ entsprechend entscheiden und den obigen Vorgang wiederholt wiederholen, erhalten Sie den lokal optimalen geschätzten Wert $ \ hat {\ theta} $. Dies ist der EM-Algorithmus.
Referenzmaterial
- CS229 Lecture notes, by Andrew NG
- Computational Statistics in Python, Expectation Maximization (EM) Algorithm
Recommended Posts