function [bhat,s,fval,H,V,scores]=estimate_logit_FE(y,x,nn,bstart,imet,print)
% bstart is a single constant followed by betas
% 


    n=length(y);
    bhat=bstart(2:length(bstart));
    
    if imet==0

        da=5;
        al=zeros(n,1)+bstart(1);

        fun=@(b)concentrate_FE(y,x,nn,b,al,imet,da);
        opti=optimoptions(@fminunc,'MaxFunctionEvaluations',10000,'MaxIterations',10000,'Display','off');
        [bhat,fval] = fminunc(fun, bhat,opti);
        if print
            diary on
            fprintf('\n  Fminunc =%9.4f \n',fval)
            diary off
        end

    %     [bhat,fval] = powell(fun,bhat,100,1.0d-6);
    %  
    %     if print
    %         diary on
    %         fprintf('\n  Powell  =%9.4f \n',fval)
    %         diary off
    %     end
    %     
    %     [bhat,fval] = fminunc(fun, bhat,opti);
    %     if print
    %         diary on
    %         fprintf('\n  Fminunc =%9.4f n\',fval)
    %         diary off
    %     end
    %     
        s=0*bhat-99.999;
        H=0;
        V=0;
        scores=0;
        
    else
        
        k=size(x,2);
        no=size(nn,1);
        ii=ones(no,1);
        indx=ones(n,1);
        f=zeros(no,1);
        for i=1:no
            ss=sum(y(nn(i,1):nn(i,2)));
            if (ss==0)||(ss==(nn(i,2)-nn(i,1)+1))
                f(i)=0;
                ii(i)=0;
                indx(nn(i,1):nn(i,2))=0;
            end
        end
        nntemp=nn(ii==1,:);
        xx=x(indx==1,:);
        yy=y(indx==1,:);
        no1=size(nntemp,1);
        nn1=zeros(size(nntemp));
        nn2(1,1)=1;
        for i=1:no1-1
            nn2(i,2)=nn2(i,1)+(nntemp(i,2)-nntemp(i,1));
            nn2(i+1,1)=nn2(i,2)+1;
        end
        nn2(no1,2)=nn2(no1,1)+(nntemp(no1,2)-nntemp(no1,1));
        
        ss=zeros(no1,1);
        for i=1:no1
            ss(i)=sum(y(nn(i,1):nn(i,2)));
        end
        
        al=zeros(no1,1)+bstart(1);
        
        ahat=zeros(length(yy),1);
        not_converged=true;
        iter=0;
        tol=1.0d-7;
        while not_converged
            bstart=bhat;
            iter=iter+1;
            for i=1:no1
                ahat(nn2(i,1):nn2(i,2))=al(i);
            end
            xb=xx*bhat;
            phat=Lambda(xb+ahat);
%             L=yy'*log(phat)+(1-yy)'*log(1-phat)

            for i=1:no1
                phati=phat(nn2(i,1):nn2(i,2));
                LLi=yy(nn2(i,1):nn2(i,2))'*log(phati)+(1-yy(nn2(i,1):nn2(i,2)))'*log(1-phati);
                da=sum(yy(nn2(i,1):nn2(i,2))-phati)/sum(phati.*(1-phati));
                al_new=al(i)+da;
                go_on=true;
                while go_on
                    phati=Lambda(xb(nn2(i,1):nn2(i,2))+al_new);
                    LLi_new=yy(nn2(i,1):nn2(i,2))'*log(phati)+(1-yy(nn2(i,1):nn2(i,2)))'*log(1-phati);
                    if LLi_new>=LLi
                        go_on=false;
                        al(i)=al_new;
                    else
                        da=da/2;
                        al_new=al(i)+da;
                    end
                end
            end
            for i=1:no1
                ahat(nn2(i,1):nn2(i,2))=al(i);
            end
            xb=xx*bhat;
            phat=Lambda(xb+ahat);
            LL=yy'*log(phat)+(1-yy)'*log(1-phat);
            db1=xx'*(yy-phat);
            db2=(repmat(phat.*(1-phat),1,k).*xx)'*xx;
            db=db2\db1;
 
            go_on=true;
            while go_on
                bnew=bhat+db;
                phat=Lambda(xx*bnew+ahat);
                LL_new=yy'*log(phat)+(1-yy)'*log(1-phat);
                if LL_new>=LL
                    go_on=false;
                    bhat=bnew;
                else
                    db=db/2;
                end
            end
            
            if iter>50 || max(abs(bstart-bhat))<tol
                not_converged=false;
            end
        end
        
        fval=LL;
        s=0*bhat-99.999;
        H=0;
        V=0;
        scores=0;
    end
end
            
                
            
        
        
  

