3
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ 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 $\endgroup$ Mar 26, 2021 at 22:10

1 Answer 1

0
$\begingroup$

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.
$\endgroup$
6
  • $\begingroup$ Tried your code, but didn't manage ... It still doesn't converge. $\endgroup$ Jun 7, 2017 at 19:42
  • $\begingroup$ Ihave result for the 2nd test:(0.49059338206625835, 0.50516900450659064, 265.45141608291095) $\endgroup$ Jun 8, 2017 at 5:18
  • $\begingroup$ 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. $\endgroup$ Jun 8, 2017 at 5:26
  • $\begingroup$ Hm. Maybe I misunderstand it. It does converge, just as in your case. But! 265 isn't 500, is it? Although thetas seem reasonable. $\endgroup$ Jun 8, 2017 at 9:14
  • $\begingroup$ 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). $\endgroup$ Jun 8, 2017 at 18:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.