function xp=Mod_ML(y,alpha)
% INPUT
% ======================
% y ........a column vector containing the frequencies  
% OUTPUT
% ======================
% xp .......solution of the problem max \{ \sum y_i \log(p_i) : p_i \geq
% \alpha/n\} where n=sum(y)

n=sum(y)
f=@(lam)sum(max(y/lam,alpha/n))-1;
lb=min(y)/2
ub=2*n*max(y)/alpha
lam=bisection(f,lb,ub,1e-10)
xp=max(y/lam,alpha/n)

function z=bisection(f,lb,ub,eps)

if (f(lb)*f(ub)>0)
    error('f(lb)*f(ub)>0')
end

iter=0;
while (ub-lb>eps)
    z=(lb+ub)/2;
    iter=iter+1;
    if(f(lb)*f(z)>0)
        lb=z;
    else
        ub=z;
    end
    fprintf('iter_numer = %3d current_sol = %2.6f \n',iter,z);
end