def EM(X, k, w, m):
    delta_max=0
    sigma = np.random.rand(2,2,2)
    g=1/2*np.ones((X.shape[0], k))
    g0=np.random.rand(X.shape[0], k)
    while True:
        delta_max=0
        for i in range(X.shape[0]):
            for j in range(k):
                g0[i,j] = np.copy(g[i,j])
                g[i,j] = w[j]* fi(X[i], m[j], sigma[j])
                sum=0
                for s in range(k):
                    sum += w[s]*fi(X[i], m[s], sigma[s])
                g[i,j]/=sum
                delta_max=np.maximum(delta_max, np.abs(g[i,j]-g0[i,j]))

        #print(delta_max)
        for j in range(k):

            m[j] = 1/(X.shape[0]*w[j])#*np.sum(g[:,j].reshape(200,1)*X, axis=0) ##
            L=0
            H=0
            for i in range(X.shape[0]):
                H+=g[i,j]*X[i]
                a= (X[i]-m[j])
                a=a.reshape(X.shape[1],1)

                L+=g[i,j]*np.dot(a,a.transpose())
                #print(np.dot(a,a.transpose()))
            #print(L)
            sigma[j] = (1/(X.shape[0]*w[j]))*L
            m[j]*=H
            w[j]=np.mean(g[:,j])
            #print (w[j])
        #print(g[:,1])
        print(delta_max)
        if delta_max < 0.0001:
            return (w, m, sigma)