Index: tol/bmath/bstat/logit.cpp
===================================================================
--- tol/bmath/bstat/logit.cpp	(revision 41d0e3cb6ed5dcce53913fc55cc68d79bf9fc9c2)
+++ tol/bmath/bstat/logit.cpp	(revision 88060096c3dc2ab3d6decfbcfc131401a5215421)
@@ -63,5 +63,5 @@
   BTimer tm(BText("Logit model (")+N+"x"+n+")");
 
-  BMatrix<BDat> dif, f(N,1), H_(n,n);
+  BMatrix<BDat> dif, f(N,1), H_(n,n), Xb(N,1);
   e.Alloc(N,1);
   p.Alloc(N,1);
@@ -74,4 +74,5 @@
   clock_t tm_1, tm_0 = BTimer::Clocks();
   BReal tm_dif;
+  BDat Lp, L1p;
 
   for(iter=0; iter<BDat::MaxIter(); iter++)
@@ -81,5 +82,5 @@
     if(iter==0)
     {
-      p.SetAllValuesTo(0);
+      Xb.SetAllValuesTo(0);
     }
     else
@@ -87,6 +88,6 @@
     //Std(BText("\nTRACE Logit B=\n")+B.Name());
     //Std(BText("\nTRACE Logit X=\n")+X.Name());
-      p = X*B;
-    //Std(BText("\nTRACE Logit p=\n")+p.Name());
+      Xb = X*B;
+    //Std(BText("\nTRACE Logit Xb=\n")+Xb.Name());
     }
     oldLlh = lnLikelyhood;
@@ -94,6 +95,13 @@
     for(i=0; i<N; i++)
     {
-      p(i,0) = 1.0/(1+Exp(-p(i,0)));
-      f(i,0) = p(i,0)*(1-p(i,0));
+      if(Xb(i,0)>=0)
+      {
+        p(i,0) = 1.0/(Exp(-Xb(i,0))+1.0);
+      }
+      else
+      {
+        p(i,0) = Exp(Xb(i,0))/(Exp(Xb(i,0))+1.0);
+      }
+      f(i,0) = p(i,0)*(1.0-p(i,0));
       e(i,0) = y(i,0)-p(i,0);
       for(j=0; j<n; j++)
@@ -105,19 +113,23 @@
         }
       }
-	    if(Abs(e(i,0))<=BDat::Tolerance()) { lnL(i,0) = 0.0;	     }
-      else if(Abs(e(i,0))==1	) 
-      {
-        #if defined(_MSC_VER) || defined(UNIX)
-        //lnL(i,0) = -INFINITY;
-        lnL(i,0) = BDat::NegInf();
-        #else
-         lnL(i,0) = -1.0/0.0;
-        #endif
-      }
-      else
-      {
-	      lnL(i,0) = y(i,0)*Log(p(i,0))+(1-y(i,0))*Log(1-p(i,0));
-      }
-      // Std(BText("\ny=")+y(i,0)+"\tp="+p(i,0)+"\te="+e(i,0)+"\tlnL="+lnL(i,0));
+	    if(Abs(e(i,0))<=BDat::Tolerance()) { lnL(i,0) = 0.0; }
+      else 
+      {
+        if(Xb(i,0)>=0)
+        {
+	        Lp  = Log(p(i,0));
+          L1p = -Xb(i,0)-Log(1+Exp(-Xb(i,0)));
+        }
+        else 
+        {
+	        Lp  =  Xb(i,0)-Log(1+Exp( Xb(i,0)));
+          L1p = Log(1-p(i,0));
+        }
+        lnL(i,0) = y(i,0)*Lp + (1-y(i,0))*L1p;
+      }
+      if(!lnL(i,0).IsFinite())
+      {
+        Std(BText("\n  Non finite log-likelihood i=")+i+"\ty="+y(i,0)+"\tXb="+Xb(i,0)+"\tp="+p(i,0)+"\tlog(p)="+Lp+"\tlog(1-p)="+L1p+"\tlnL="+lnL(i,0));
+      }
       lnLikelyhood += lnL(i,0);
     }
@@ -141,9 +153,10 @@
 
     if(advance.IsUnknown()) { break; }
-    B -= dif;
     if(norm	 < DEpsilon()	    ) { break; }
     if(maxAbsDif < BDat::Tolerance()) { break; }
     if(advance	 < DEpsilon()	    ) { break; }
     if(Abs(oldLlh-lnLikelyhood) < DEpsilon()) { break; }
+    if(iter==BDat::MaxIter()-1) { break; }
+    B -= dif;
   }
   return BTRUE;
@@ -194,4 +207,5 @@
   BReal tm_dif;
   lnLikelyhood = BDat::NegInf();
+  BDat _p, _1p;
   for(iter=0; iter<BDat::MaxIter(); iter++)
   {
@@ -203,8 +217,10 @@
     for(i=0; i<N; i++)
     {
-      p(i,0) = BNormalDist::Dist01(XB(i,0));
+      _1p = BNormalDist::Dist01(-XB(i,0));
+       _p = BNormalDist::Dist01( XB(i,0));
+      p(i,0) = _p;
       f(i,0) = Exp(-0.5*(XB(i,0)*XB(i,0)))/Sqrt(2.0*BDat::Pi());//BNormalDist::Dens01(XB(i,0));
-      l(i,0) = (y(i,0)==0)?-f(i,0)/(1-p(i,0)):f(i,0)/p(i,0);
-      e(i,0) = y(i,0)-p(i,0);
+      l(i,0) = (y(i,0)==0)?-f(i,0)/(_1p):f(i,0)/_p;
+      e(i,0) = y(i,0)-_p;
       for(j=0; j<n; j++)
       {
@@ -216,14 +232,10 @@
       }
 	    if(Abs(e(i,0))<=BDat::Tolerance()) { lnL(i,0) = 0.0;	     }
-      else if(Abs(e(i,0))==1		    ) 
-      {
-        lnL(i,0) = BDat::NegInf();
-      }
       else
       {
-	      lnL(i,0) = y(i,0)*Log(p(i,0))+(1-y(i,0))*Log(1-p(i,0));
+	      lnL(i,0) = y(i,0)*Log(_p)+(1-y(i,0))*Log(_1p);
       }
 //    Std(BText("\ny=")+y(i,0)+"\tXB="+XB(i,0)+
-//		"\tp="+p(i,0)+"\tf="+f(i,0)+"\tl="+l(i,0)+
+//		"\tp="+_p+"\tf="+f(i,0)+"\tl="+l(i,0)+
 //		"\te="+e(i,0)+"\tlnL="+lnL(i,0));
       lnLikelyhood += lnL(i,0);
@@ -247,9 +259,10 @@
 
     if(advance.IsUnknown()) { break; }
-    B -= dif;
     if(norm	     < chop 	          ) { break; }
     if(maxAbsDif < BDat::Tolerance()) { break; }
     if(advance	 < chop     	      ) { break; }
     if(Abs(oldLlh-lnLikelyhood) < DEpsilon()) { break; }   
+    if(iter==BDat::MaxIter()-1) { break; }
+    B -= dif;
   }
   return BTRUE;