oldlik = -inf;
newlik = -1e100;

eye_d = eye(d);
eye_k = eye(k);

[lambda_psi_inv, lambda_t, psi_inv] = lambda_psi_inverse(lambda, psi, eye_d, eye_k);

while (newlik-oldlik)/abs(newlik) > converge_ratio
    step = step + 1;

    beta = lambda'*lambda_psi_inv;


    [lambda, psi] = mstep(beta, s, lambda, psi, eye_k);
    [lambda_psi_inv, lambda_t, psi_inv] = lambda_psi_inverse(lambda, psi, eye_d, eye_k);

    oldlik = newlik;
    newlik = likelihood(lambda, lambda_t, psi, psi_inv, lambda_psi_inv, s, eye_d);
    likhist(step) = newlik;
    if newlik < oldlik
        error('new likelihood worse than old likelihood');
    end
    disp([step, newlik]);
end