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?
 A: 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.

