root/tolp/trunk/tol/bmath/bstat/logit.cpp @ 203

Revision 203, 12.6 kB (checked in by vdebuen, 5 years ago)

At first iteration of Logit is B=0 and so p=X*B=0

  • Property svn:eol-style set to native
Line 
1/* logit.cpp: logistic model.
2              GNU/TOL Language.
3
4   Copyright (C) 2003 - Bayes Decision, SL (Spain [EU])
5
6   This program is free software; you can redistribute it and/or modify
7   it under the terms of the GNU General Public License as published by
8   the Free Software Foundation; either version 2, or (at your option)
9   any later version.
10
11   This program is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14   GNU General Public License for more details.
15
16   You should have received a copy of the GNU General Public License
17   along with this program; if not, write to the Free Software
18   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
19   USA.
20 */
21
22#if defined(_MSC_VER)
23#include <win_tolinc.h>
24#endif
25
26#include <tol/tol_matrix.h>
27#include <tol/tol_blogit.h>
28#include <tol/tol_blinalg.h>
29#include <tol/tol_bmfstpro.h>
30#include <tol/tol_bstat.h>
31#include <tol/tol_bltrimat.h>
32#include <tol/tol_btimer.h>
33
34BTraceInit("logit.cpp");
35
36
37//--------------------------------------------------------------------
38BBool Logit(const BMatrix    <BDat>& y,
39            const BMatrix    <BDat>& X,
40                  BMatrix    <BDat>& B,
41                  BMatrix    <BDat>& e,
42                  BMatrix    <BDat>& p,
43                  BMatrix    <BDat>& G,
44                  BSymMatrix <BDat>& H,
45                  BMatrix    <BDat>& lnL,
46                  BDat&              lnLikelyhood,
47                  BDat&              likelyhood)
48
49/*! Estimates the boolean model
50 *
51 *              y = F(X * B) + e;
52 *
53 *                dim(B) = n   : number of variables
54 *                dim(y) = N   : number of data
55 *                dim(X) = N,n
56 *                dim(e) = N
57 *
58 *  where y is 0 or 1 and F is the logistic distribution.
59 */
60//--------------------------------------------------------------------
61{
62  BInt i, j, k, iter, n=X.Columns(), N=X.Rows();
63  BTimer tm(BText("Logit model (")+N+"x"+n+")");
64
65  BMatrix<BDat> dif, f(N,1), H_(n,n);
66  e.Alloc(N,1);
67  p.Alloc(N,1);
68  lnL.Alloc(N,1);
69  B.Alloc(n,1);
70  G.Alloc(n,1);
71  H.Alloc(n);
72  BDat norm, advance, maxAbsDif, oldLlh;
73  B.SetAllValuesTo(0);
74  clock_t tm_1, tm_0 = BTimer::Clocks();
75  BReal tm_dif;
76
77  for(iter=0; iter<BDat::MaxIter(); iter++)
78  {
79    H_.SetAllValuesTo(0);
80    G.SetAllValuesTo(0);
81    if(iter==0)
82    {
83      p.SetAllValuesTo(0);
84    }
85    else
86    {
87    //Std(BText("\nTRACE Logit B=\n")+B.Name());
88    //Std(BText("\nTRACE Logit X=\n")+X.Name());
89      p = X*B;
90    //Std(BText("\nTRACE Logit p=\n")+p.Name());
91    }
92    oldLlh = lnLikelyhood;
93    lnLikelyhood=0;
94    for(i=0; i<N; i++)
95    {
96      p(i,0) = 1.0/(1+Exp(-p(i,0)));
97      f(i,0) = p(i,0)*(1-p(i,0));
98      e(i,0) = y(i,0)-p(i,0);
99      for(j=0; j<n; j++)
100      {
101              G(j,0) += X(i,j)*e(i,0);
102              for(k=0; k<n; k++)
103        {
104                H_(j,k) -= f(i,0)*X(i,j)*X(i,k);
105        }
106      }
107            if(Abs(e(i,0))<=BDat::Tolerance()) { lnL(i,0) = 0.0;             }
108      else if(Abs(e(i,0))==1    ) 
109      {
110        #if defined(_MSC_VER) || defined(UNIX)
111        //lnL(i,0) = -INFINITY;
112        lnL(i,0) = BDat::NegInf();
113        #else
114         lnL(i,0) = -1.0/0.0;
115        #endif
116      }
117      else
118      {
119              lnL(i,0) = y(i,0)*Log(p(i,0))+(1-y(i,0))*Log(1-p(i,0));
120      }
121      // Std(BText("\ny=")+y(i,0)+"\tp="+p(i,0)+"\te="+e(i,0)+"\tlnL="+lnL(i,0));
122      lnLikelyhood += lnL(i,0);
123    }
124    H = H_;
125    //dif           = MinimumResidualsSolve(H_,G);
126      //MinimumResidualsSolve(H_,G, BDat(DEpsilon()));
127    dif = gsl_MinimumResidualsSolve(H_,G);
128    norm    = G.FrobeniusNorm();
129    advance = dif.FrobeniusNorm();
130    maxAbsDif = dif.MaxAbsNorm();
131    tm_1 = BTimer::Clocks();
132    tm_dif = BTimer::ClockToSecond(tm_1-tm_0);
133    tm_0 = tm_1;
134    //likelyhood = Exp(lnLikelyhood);
135    Std(BText("\n  Logit model iteration(")+iter+") "+
136              "  LogLikelyhood = "+lnLikelyhood+
137        //"      Likelyhood = "+likelyhood+
138              "  maxAbsDif = "+maxAbsDif+
139              "  Gradient Norm = "+norm+
140          I2("\tTime : "," Tiempo : ") + tm_dif+I2(" seconds"," segundos"));
141
142    if(advance.IsUnknown()) { break; }
143    B -= dif;
144    if(norm      < DEpsilon()       ) { break; }
145    if(maxAbsDif < BDat::Tolerance()) { break; }
146    if(advance   < DEpsilon()       ) { break; }
147    if(Abs(oldLlh-lnLikelyhood) < DEpsilon()) { break; }
148  }
149  return BTRUE;
150}
151
152
153
154//--------------------------------------------------------------------
155BBool Probit(const BMatrix    <BDat>& y,
156                   const BMatrix    <BDat>& X,
157                               BMatrix    <BDat>& B,
158                               BMatrix    <BDat>& e,
159                               BMatrix    <BDat>& p,
160                               BMatrix    <BDat>& G,
161                               BSymMatrix <BDat>& H,
162                               BMatrix    <BDat>& lnL,
163                               BDat&                lnLikelyhood,
164                               BDat&                likelyhood)
165
166/*! Estimates the boolean model
167 *
168 *              y = F(X * B) + e;
169 *
170 *                dim(B) = n   : number of variables
171 *                dim(y) = N   : number of data
172 *                dim(X) = N,n
173 *                dim(e) = N
174 *
175 *  where y is 0 or 1 and F is the logistic distribution.
176 */
177//--------------------------------------------------------------------
178{
179  BInt i, j, k, iter, n=X.Columns(), N=X.Rows();
180  BTimer tm(BText("Probit model (")+N+"x"+n+")");
181  double chop = Sqrt(DEpsilon()/N);
182
183  BMatrix<BDat> dif, f(N,1), l(N,1), H_(n,n), XB;
184  e.Alloc(N,1);
185  p.Alloc(N,1);
186  lnL.Alloc(N,1);
187  B.Alloc(n,1);
188  G.Alloc(n,1);
189  H.Alloc(n);
190  p.Alloc(N,1);
191  BDat norm, advance, maxAbsDif, oldLlh;
192  B.SetAllValuesTo(0);
193  clock_t tm_1, tm_0 = BTimer::Clocks();
194  BReal tm_dif;
195  lnLikelyhood = BDat::NegInf();
196  for(iter=0; iter<BDat::MaxIter(); iter++)
197  {
198    H_.SetAllValuesTo(0);
199    G.SetAllValuesTo(0);
200    XB = X*B;
201    oldLlh = lnLikelyhood;
202    lnLikelyhood=0;
203    for(i=0; i<N; i++)
204    {
205      p(i,0) = BNormalDist::Dist01(XB(i,0));
206      f(i,0) = Exp(-0.5*(XB(i,0)*XB(i,0)))/Sqrt(2.0*BDat::Pi());//BNormalDist::Dens01(XB(i,0));
207      l(i,0) = (y(i,0)==0)?-f(i,0)/(1-p(i,0)):f(i,0)/p(i,0);
208      e(i,0) = y(i,0)-p(i,0);
209      for(j=0; j<n; j++)
210      {
211              G(j,0) += X(i,j)*l(i,0);
212              for(k=0; k<n; k++)
213              {
214                H_(j,k) -= l(i,0)*(l(i,0)+XB(i,0))*X(i,j)*X(i,k);
215              }
216      }
217            if(Abs(e(i,0))<=BDat::Tolerance()) { lnL(i,0) = 0.0;             }
218      else if(Abs(e(i,0))==1                ) 
219      {
220        lnL(i,0) = BDat::NegInf();
221      }
222      else
223      {
224              lnL(i,0) = y(i,0)*Log(p(i,0))+(1-y(i,0))*Log(1-p(i,0));
225      }
226//    Std(BText("\ny=")+y(i,0)+"\tXB="+XB(i,0)+
227//              "\tp="+p(i,0)+"\tf="+f(i,0)+"\tl="+l(i,0)+
228//              "\te="+e(i,0)+"\tlnL="+lnL(i,0));
229      lnLikelyhood += lnL(i,0);
230    }
231    H = H_;
232  //dif     = MinimumResidualsSolve(H_,G);
233    dif       = gsl_MinimumResidualsSolve(H_,G);
234    norm      = G.FrobeniusNorm();
235    advance   = dif.FrobeniusNorm();
236    maxAbsDif = dif.MaxAbsNorm();
237    tm_1      = BTimer::Clocks();
238    tm_dif    = BTimer::ClockToSecond(tm_1-tm_0);
239    tm_0      = tm_1;
240//  likelyhood = Exp(lnLikelyhood);
241    Std(BText("\n  Probit model iteration(")+iter+") "+
242              "  LogLikelyhood = "+lnLikelyhood+
243//            "  Likelyhood = "+likelyhood+
244              "  maxAbsDif = "+maxAbsDif+
245              "  Gradient Norm = "+norm+
246             I2("\tTime : "," Tiempo : ") + tm_dif+I2(" seconds"," segundos"));
247
248    if(advance.IsUnknown()) { break; }
249    B -= dif;
250    if(norm          < chop               ) { break; }
251    if(maxAbsDif < BDat::Tolerance()) { break; }
252    if(advance   < chop               ) { break; }
253    if(Abs(oldLlh-lnLikelyhood) < DEpsilon()) { break; }   
254  }
255  return BTRUE;
256}
257
258
259
260//--------------------------------------------------------------------
261BInt BooleanModelIteration(const BArray     <BDat>& y,
262                           const BMatrix    <BDat>& X,
263                                 BArray     <BDat>& B,
264                                 BArray     <BDat>& e,
265                                 BSymMatrix <BDat>& cov,
266                                 BArray     <BDat>& p,
267                                 BArray     <BDat>& yPrev,
268                                 BDat&              s,
269                                 BDat&              R2,
270                                 BDat&              F,
271                                 BInt&              v1,
272                                 BInt&              v2,
273                                 BDat&              alfa,
274                                 BDat&              resSqSum,
275                                 BDat               tolerance,
276                                 BProbDist*         dist)
277//--------------------------------------------------------------------
278{
279  BArray<BDat> T,b,err;
280  BInt m = X.Columns();
281  BInt n = X.Rows();
282  BInt i, j;
283  BMatrix<BDat> C = X;
284  BArray<BDat>  u = y;
285  BArray<BDat>  v = y;
286  BDat g, absb;
287
288  BDat sOld = s;
289  BDat prop;
290
291  for(i=0; i<n; i++)
292  {
293    v(i) = Sqrt(p(i)*(1.0-p(i)));
294    u(i) = (y(i)-p(i))/v(i);
295    for(j=0; j<m; j++)
296    {
297      C(i,j) *= v(i);
298    }
299  }
300  C*=(-1);
301  if(!HouseholderTransformation(T,C))
302  {
303    return(-1);
304  }
305
306        LeastSqrHouseholder(T,C,u,b,err);
307
308  for(i=0; i<m; i++) { B(i) += b(i); }
309
310  e = y;
311  prop = s = 0;
312  for(i=0; i<n; i++)
313  {
314    g = 0;
315    for(j=0; j<m; j++) { g += B(j) * X(i,j); }
316    p(i) = dist->Dist(g);
317    yPrev(i) = Round(p(i));
318    if(yPrev(i)!=y(i)) { prop+=1; }
319    e(i) -= p(i);
320    s += e(i) ^ 2;
321  }
322
323  prop /= n;
324  s /= n;
325  s.Sqrt();
326
327  BLowTrMatrix<BDat> Rt(m);
328
329  for(i=0; i<m; i++)
330  {
331    for(j=0; j<i; j++) { Rt(i,j)=C(j,i); }
332    Rt(i,i) = T(i);
333  }
334
335  BLowTrMatrix<BDat> Rti;
336  Rt.Inverse(Rti);
337  cov = MtMSqr(Rti);
338
339  cov *= s*s;
340  for(i=0; i<n; i++)
341  {
342    for(j=0; j<m; j++)
343    {
344      cov(i,j) /= v(i);
345    }
346  }
347  BArray<BDat> yp = y;
348  for(i=0; i<n; i++) { yp(i) += e(i); }
349  R2 = R2Coefficient(y,yp);
350  v1 = m;
351  v2 = n-m-1;
352  F  = R2*BDat(v2) / ((1-R2)*BDat(v1));
353  alfa = BFSnedecorDist(v1,v2).Dist(F);
354  resSqSum = Moment(e,2);
355
356
357  BDat ds = sOld-s;
358  BDat d  = FrobeniusNorm(b);
359/*
360  Std(BText("\n######################################################\n"));
361  Std(BText("\nB=\n")+Name(B));
362//Std(BText("\ne=\n")+Name(e));
363//Std(BText("\nP=\n")+Name(p));
364  Std(BText("\nprop=\n")+prop.Name());
365  Std(BText("\ns=\n")+s.Name());
366  Std(BText("\nd=\n")+d.Name());
367  Std(BText("\ntolerance=\n")+tolerance.Name());
368  Std(BText("\n######################################################\n"));
369*/
370  return((d<=tolerance) && (ds>=0) && (ds<=tolerance));
371}
372
373
374
375//--------------------------------------------------------------------
376BBool BooleanModel(const BArray     <BDat>& y,
377                   const BMatrix    <BDat>& X,
378                         BArray     <BDat>& B,
379                         BArray     <BDat>& e,
380                         BSymMatrix <BDat>& cov,
381                         BArray     <BDat>& p,
382                         BArray     <BDat>& yPrev,
383                         BDat&              s,
384                         BDat&              R2,
385                         BDat&              F,
386                         BInt&              v1,
387                         BInt&              v2,
388                         BDat&              alfa,
389                         BDat&              resSqSum,
390                         BDat&              totSqSum,
391                         BDat               tolerance,
392                         BProbDist*         dist)
393   
394/*! Estimates the boolean model
395 *
396 *              y = F(X * A) + e;
397 *
398 *                dim(B) = m   : number of variables
399 *                dim(y) = n   : number of data
400 *                dim(X) = n,m
401 *                dim(e) = n
402 *
403 *              where y is 0 or 1.
404 */
405//--------------------------------------------------------------------
406{
407  BInt i, j, result;
408  BInt n = X.Rows();
409  BInt m = X.Columns();
410  p.ReallocBuffer(n);
411  B.ReallocBuffer(m);
412  BInt iter = 0;
413  totSqSum = Varianze(y);
414  s = Sqrt(totSqSum);
415
416  for(i=0; i<n; i++)
417  {
418    p(i) = 0.5;
419    if((y(i)!=0.0)&&(y(i)!=1.0))
420    {
421      Error(I2("There non 0 or 1 values for output variable in a Logistic "
422               "Model.",
423               "Existen valores distintos de 0 ó 1 en la variable respuesta "
424               "de un Modelo Logístico."));
425      return(BFALSE);
426    }
427  }
428  for(j=0; j<m; j++)
429  {
430    B(j) = 0.0;
431  }
432
433  BTimer tmr(I2("Logistic estimation","Estimación logística"));
434  Std(I2("\n Number of variables ", "\nNúmero de variables ") + n);
435  Std(I2("\n Number of data ", "\nNúmero de datos ") + m);
436  Std(I2("\n Standard deviation "," Desviación estandard ") + s);
437  yPrev.ReallocBuffer(n);
438  do
439  {
440    iter++;
441    result=BooleanModelIteration(y,X,B,e,cov,p,yPrev,
442                                 s,R2,F,v1,v2,alfa,resSqSum,tolerance,dist);
443    Std(I2("\nIteration ","\nIteración ")+iter+
444        I2(" Standard deviation "," Desviación estandard ") + s);
445  }
446  while(!result && (iter<=BDat::MaxIter()));
447
448  return(result==1);
449}
450
451
452
453
454/*
455//--------------------------------------------------------------------
456   BBool Probit(const BArray     <BDat>& y,
457                const BMatrix    <BDat>& X,
458                      BArray     <BDat>& B,
459                      BArray     <BDat>& e,
460                      BSymMatrix <BDat>& cov,
461                      BArray     <BDat>& p,
462                      BArray     <BDat>& yPrev,
463                                  BDat & s,
464                                  BDat & R2,
465                                  BDat & F,
466                                  BInt & v1,
467                                  BInt & v2,
468                                  BDat & alfa,
469                                  BDat & resSqSum,
470                                  BDat & totSqSum,
471                                  BDat   tolerance)
472
473// PURPOSE: Estimates the boolean model
474//
475//              y = F(X * A) + e;
476//
477//                dim(B) = m   : number of variables
478//                dim(y) = n   : number of data
479//                dim(X) = n,m
480//                dim(e) = n
481//
482//              where y is 0 or 1.
483//
484//--------------------------------------------------------------------
485{
486  static BNormalDist dist;
487  return(BooleanModel(y,X,B,e,cov,p,yPrev,s,
488                      R2,F,v1,v2,alfa,resSqSum,totSqSum,tolerance,&dist));
489}
490*/
491
Note: See TracBrowser for help on using the browser.