# Another EM-algorithm problem

I have the following problem:

I have a random vector $y$ which has length $l$. The first $z$ bits come from a Bernoulli random variable with parameter $\theta_1$ and the next $l-z$ come from a Bernoulli random variable with parameter $\theta_2$.

I need to find the values of $\theta_1$ and $\theta_2$. I don't need $z$ but expect to find it anyway.

My idea of the solution:

E-step: given some values of $\theta_1$ and $\theta_2$, find $E[z]$.

This seems to be a double sum: $\frac{\sum_{z=0}^l z*\theta_1^{\sum_{i=0}^z y_i}*(1-\theta_1)^{z - \sum_{i=0}^z y_i}*\theta_2^{\sum_{i=z+1}^l y_i}*(1-\theta_2)^{l-z - \sum_{i=z+1}^l y_i}}{\sum_{z=0}^l \theta_1^{\sum_{i=0}^z y_i}*(1-\theta_1)^{z - \sum_{i=0}^z y_i}*\theta_2^{\sum_{i=z+1}^l y_i}*(1-\theta_2)^{l-z - \sum_{i=z+1}^l y_i}}$

This monster of a formula is essentially just multiplying every possible value of $z$ by the likelihood of the sample given this $z$ divided by the sum of likelihoods themselves. Seems to be the expectation of $z$. Just two nested loops.

M-step: given $z$, find $\theta_1$ and $\theta_2$.

Well, theoretically we would need to compute an $argmax$, but here, with just one parameter for every piece of the vector, we have:

$\theta_1 = {\sum_{i=0}^z y_i}/z$

$\theta_2 = {\sum_{i=z}^l y_i}/(l-z)$

I have some python code, which does exactly that:

import numpy as np
import random
import matplotlib.pyplot as plt
import math
%matplotlib inline

def e_step( data, theta1, theta2 ):#quadratic time, but okay for l=550
l = len(data)
z = 0
final_likelihood = 0;
#in the next loop I am computing the expectation of Z
#by multiplying every possible z by it's likelihood
for zi in range(0,l):
likelihood = 1
for i in range(0, zi):#on these elements we suspect to have p=\theta_1
if data[i]==1:
likelihood = likelihood*theta1
else:
likelihood = likelihood*(1-theta1)
for i in range(zi+1, l):#on these elements we suspect to have p=\theta_2
if data[i]==1:
likelihood = likelihood*theta2
else:
likelihood = likelihood*(1-theta2)
final_likelihood = final_likelihood+likelihood
z = z + zi*likelihood
z1 = z/final_likelihood #\sum z_i * p(z_i)
return math.floor(z1)

def m_step( data, z ):
l = len(data)
sum1=0;
for i in range(0,z):
sum1=sum1+data[i]#maximum likelihood \theta_1 is just the mean value (Bernoulli)
if z<=1:
theta1=0
else:
theta1=sum1/(z-1);
sum2=0
for i in range(z+1,l):#maximum likelihood \theta_2 is just the mean value (Bernoulli)
sum2=sum2+data[i]
if (l-z)<=1:
theta2=0
else:
theta2=sum2/(l-z-1);

return theta1, theta2

def em_algorithm( data, n_iter=50):
theta1 = 0.5
theta2 = 0.5
z = e_step(data, theta1, theta2 )
theta1_array=[]
theta2_array=[]
z_array = []
plot_num = n_iter

for i in range(0, n_iter):
(theta1, theta2) = m_step( data, z)
z = e_step(data, theta1, theta2 )
likelihood=0
theta1_array.append(theta1*100)#we scale \theta by 100 to fit \theta and z on one plot
theta2_array.append(theta2*100)
z_array.append(z)

x = np.arange(plot_num)
l1, = plt.plot( x, theta1_array, label='theta1' )
l2, = plt.plot( x, theta2_array, label='theta2' )
l3, = plt.plot( x, z_array , label='z')

plt.legend(handles=[l1, l2, l3])
plt.show()

return theta1, theta2, z

#Test 1
z=50
theta1=0.1
theta2=0.9
l=100
a = np.random.binomial(1, theta1, size=z)
b = np.random.binomial(1, theta2, size=l-z)
data = np.concatenate((a,b))
random.randint( 0, len(data))
em_algorithm( data, 10 )

#Test 2
z=500
theta1=0.49
theta2=0.51
l=550
a = np.random.binomial(1, theta1, size=z)
b = np.random.binomial(1, theta2, size=l-z)
data = np.concatenate((a,b))
random.randint( 0, len(data))
a = em_algorithm( data, 50 )


Well, the algorithm doesn't work. It does converge to some random numbers in 2-3 steps.

Could someone have a look at either the formulae or the code?

• This isn't the EM algorithm: the E step of the EM algorithm computes the expectation of the loglikelihood, not of the missing data. In a lot of simple examples the expectation of the loglikelihood is what you get by plugging in the expectation of the missing data, but not in general -- and I think not here Mar 26, 2021 at 22:10

Seems you have a float-int misuse. Try to change you e_step by returning z as float. Then in m_step convert z to int on demand (in loops). See here the patch:

def e_step( data, theta1, theta2 ):#quadratic time, but okay for l=550
l = len(data)
z = 0
final_likelihood = 0;
#in the next loop I am computing the expectation of Z
#by multiplying every possible z by it's likelihood
for zi in range(0,l):
likelihood = 1
for i in range(0, zi):#on these elements we suspect to have p=\theta_1
if data[i]==1:
likelihood = likelihood*theta1
else:
likelihood = likelihood*(1-theta1)
for i in range(zi+1, l):#on these elements we suspect to have p=\theta_2
if data[i]==1:
likelihood = likelihood*theta2
else:
likelihood = likelihood*(1-theta2)

final_likelihood = final_likelihood+likelihood
z = z + zi*likelihood
z1 = z/final_likelihood #\sum z_i * p(z_i)
#     return math.floor(z1)
return z1

def m_step( data, z ):
l = len(data)
sum1=0;
print(z)
for i in range(0,int(z)):
sum1=sum1+data[i]#maximum likelihood \theta_1 is just the mean value (Bernoulli)
if z<=1:
theta1=0
else:
theta1=sum1/(z-1);
sum2=0
for i in range(int(z)+1,l):#maximum likelihood \theta_2 is just the mean value (Bernoulli)
sum2=sum2+data[i]
if (l-z)<=1:
theta2=0
else:
theta2=sum2/(l-z-1);

return theta1, theta2


So, I got results for the 1st test:

(0.082442583084907359, 0.94985216350816926, 49.518615627077224)


2nd test:

(0.49059338206625835, 0.50516900450659064, 265.45141608291095)


I've tested on python notebook with python 2.7.

Your initial version of code didn't work with me due to error

<ipython-input-1-dda33dc9ff4c> in m_step(data, z)
31 l = len(data)
32 sum1=0;
---> 33 for i in range(0,z):
34 sum1=sum1+data[i]#maximum likelihood \theta_1 is just the mean value (Bernoulli)
35 if z<=1:
TypeError: range() integer end argument expected, got float.

• Tried your code, but didn't manage ... It still doesn't converge. Jun 7, 2017 at 19:42
• Ihave result for the 2nd test:(0.49059338206625835, 0.50516900450659064, 265.45141608291095) Jun 8, 2017 at 5:18
• Hmm, my launch of the code converges fast. I've added the info about my test into main answer, hope you'll find the key. Jun 8, 2017 at 5:26
• Hm. Maybe I misunderstand it. It does converge, just as in your case. But! 265 isn't 500, is it? Although thetas seem reasonable. Jun 8, 2017 at 9:14
• I think it is because theta1=0.51, theta2=0.49 are similiar in values. As well as you have just l=550 (it is not a lot of statistics for making theta values distinguishable enough). Jun 8, 2017 at 18:53