function [beta,s,f,H,V,scores]=estimate_logit(y,x)


sx=std(x);
ii=find(sx==0);
theta=4*x\y;
theta(ii)=theta(ii)-2;
npar=length(theta);
sdf=ones(1,npar);
iter=0;

iter=iter+1;
while ((mean(abs(sdf))>1.0d-4) && (iter<250))
        beta=theta;
        xb=x*beta;
        xhess=zeros(npar,npar);
        exb=exp(trunbo(xb,30));
        p=exb./(1+exb);
        f=y'*log(p)+(1-y)'*log(1-p);
        sdf=x'*(y-p);
        xhess=-(repmat((p.*(1-p)),1,npar).*x)'*x;

        theta=theta+0.5*inv(-xhess+eye(npar)*1.e-12)*sdf;
        theta=trunbo(theta,20);
 end

 beta=theta;
 xb=x*beta;
 xhess=zeros(npar,npar);
 exb=exp(xb);
 p=exb./(1+exb);
 f=y'*log(p)+(1-y)'*log(1-p);

 df=x.*repmat((y-p),1,npar);
 H=-(repmat((p.*(1-p)),1,npar).*x)'*x;

Hinv=inv(-H+eye(npar)*1.0d-16);

scores=df;
V=scores'*scores;
s=sqrt(diag(Hinv'*V*Hinv));


end




    