close Warning: Can't synchronize with repository "(default)" (/var/svn/mms does not appear to be a Subversion repository.). Look in the Trac log for more information.

Ticket #63: nameblock_logit.tol

File nameblock_logit.tol, 14.7 KB (added by Pedro Gea, 16 years ago)
Line 
1//////////////////////////////////////////////////////////////////////////////
2// FILE    :
3// PURPOSE :
4//////////////////////////////////////////////////////////////////////////////
5
6
7//_.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-.
8NameBlock Logit =
9//_.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-.
10[[
11   Real _.ctInfo           = 0;
12   Real _.dEpsilon         = DiffDist;
13   Real _.maxIter          = MaxIter;
14   Real _.tolerance        = 10^(-8);
15   Real _.tolerance.Rec    = 10^(-4);
16   Real _.probability.Init = 1/2;
17   Real _.num.Step         = 8;
18
19   Matrix _Error(Matrix y, Matrix p)
20   {
21     y-p //p = Probability(X, B)
22   };
23   
24   Matrix _Gradient(Matrix X, Matrix error)
25   {
26     Tra(Tra(error)*X) //error = y - Probability(X, B)
27   };
28
29   Matrix _Hessian(VMatrix vX, Matrix p)
30   {
31     VMatrix W = Mat2VMat(p$*RSum(-p,1));
32     VMatrix V = Eye(Rows(p), Rows(p), 0, W);
33     VMat2Mat(-Tra(V*vX)*vX)     
34   };
35
36   Matrix _Probability(Matrix X, Matrix B, Matrix LogitGroupProb)
37   {
38     Matrix p = RPow(RSum(Exp(-X*B-LogitGroupProb), 1), -1);
39     p
40   };
41
42   Matrix _MLnLikelyhood(Matrix y, Matrix p)
43   {
44     y$*Log(p)+RSum(-y,1)$*Log(RSum(-p,1))
45   };
46
47   Matrix _MLikelyhood(Matrix mLnLikelyhood)
48   {
49     Exp(mLnLikelyhood)
50   };
51
52   Set Estimate.MaxLikelyhood
53   (
54     Matrix y,        // Matriz y de variable salida binaria
55     Matrix XIni,     // Matriz x de variables entrada
56     Matrix B0Ini,    // Matriz de parametros iniciales
57     Anything ctInfo, // Un real 0 sin constante,
58                      // Un real 1 estima constante,
59                      // Una matriz 
60     Real dEpsilon,   // Diferencia de paso
61     Real maxIter,    // Numero maximo de iteraciones
62     Real tolerance   // Tolerancia al error
63   )
64   {
65     Text grCtInfo = Grammar(ctInfo);
66     
67     Matrix X = Case
68     (
69/*
70       And(grCtInfo=="Real", EQ(ctInfo, 0)), XIni,
71       And(grCtInfo=="Real", EQ(ctInfo, 1)),
72       {
73         XIni|Rand(Rows(XIni), 1, 1, 1)
74       },
75*/
76       grCtInfo=="Matrix", XIni,
77       1, XIni
78     );
79
80     Matrix IdCt = Case
81     (
82//       And(grCtInfo=="Real", EQ(ctInfo, 0)), Rand(Rows(XIni), 1, 0, 0),
83//       And(grCtInfo=="Real", EQ(ctInfo, 1)), Rand(Rows(XIni), 1, 0, 0),
84       grCtInfo=="Matrix", ctInfo,
85       1, Rand(Rows(XIni), 1, 0, 0)
86     );
87       
88     Matrix B0 = Case
89     (
90       grCtInfo=="Matrix", B0Ini,
91       1, B0Ini
92     );
93     
94     Real n = Columns(X); // Number of variables
95     Real N = Rows(X);    // Number of data
96
97     Text iniMsg =
98      "Model Logit Init ("+IntText(N)+"x"+IntText(n)+")"+Time;
99     Real WriteLn(iniMsg); 
100
101     VMatrix vX       = Mat2VMat(X);
102
103     Matrix zeroMat   = Rand(N, 1, 0, 0);
104     Matrix oneMat    = Rand(N, 1, 1, 1);
105     Matrix negInfMat = Rand(N, 1, -1/0, -1/0); 
106     Matrix tolMat    = Rand(N, 1, tolerance, tolerance);   
107     
108     Real completeTime    = 0;
109     Real difTimeMin      = 0;
110     Real exit            = 0;
111     Matrix B             = Copy(B0);
112     Real iter            = 0;
113     Real oldLnLikelyhood = 0;
114   
115     Set cycle = Empty;
116     
117     Real While(Not(exit),
118     {
119       Real time0 = Copy(Time);
120
121//Real CMsg::Trace::show(1, "Printing oldLnLikelyhood="<<oldLnLikelyhood);
122       Matrix p     = _Probability(X, B, IdCt);
123
124//Real CMsg::Trace::show(1, "Printing p="<<p);
125       Matrix error = _Error(y, p);
126
127       Matrix G     = _Gradient(X, error);
128       Real time0.1 = Copy(Time);
129       Matrix H     = _Hessian(vX, p);
130       Real time0.2 = Copy(Time); 
131       Matrix dif   = MinimumResidualsSolve(H, G);
132 
133       Real norm      = MatFrobeniusNorm(G);
134       Real advance   = MatFrobeniusNorm(dif);
135       Real maxAbsDif = MatMax(Abs(dif));
136       
137
138       Matrix mLnLikelyhood = _MLnLikelyhood(y, p);
139       Matrix mLLCorrect    = IfMat(EQ(Abs(error), oneMat), negInfMat,
140       IfMat(LT(Abs(error), tolMat), zeroMat, mLnLikelyhood));
141
142       Real lnLikelyhood     = MatSum(mLLCorrect);
143//Real CMsg::Trace::show(1, "Printing dif"<<dif);
144//Real CMsg::Trace::show(1, "Printing mLnLikelyhood"<<mLnLikelyhood);
145//Real CMsg::Trace::show(1, "Printing mLLCorrect"<<mLLCorrect);
146//Real CMsg::Trace::show(1, "Printing error"<<error);
147
148//Real CMsg::Trace::show(1, "Printing lnLikelyhood="<<lnLikelyhood);
149
150       Real difTimeMin := (time0.2-time0.1)+difTimeMin;
151       
152       Real difTime = Copy(Time)-time0;
153       Real completeTime := completeTime+difTime;
154       Real exitTreatment = Case
155       ( 
156         IsUnknown(advance),
157         {
158           Real CMsg::Trace::show(1, "Advance is unknown");
159           Real exit:=1
160         },
161         GE(iter, maxIter),
162         {
163           Real CMsg::Trace::show(1, "MaxIter reached");
164           Real exit:=1
165         },
166         LT(norm, dEpsilon),
167         {
168           Real CMsg::Trace::show(1, "Norm lower than "<<dEpsilon);
169           Real exit:=1
170         },
171         LT(maxAbsDif, tolerance),
172         {
173           Real CMsg::Trace::show(1, "MaxAbsDif lower than "<<tolerance);
174           Real exit:=1
175         },
176         LT(Abs(oldLnLikelyhood - lnLikelyhood), dEpsilon),
177         {
178           Real CMsg::Trace::show(1,
179           "LnLikelihood is stable: lnLikelyhood    ="<<lnLikelyhood+
180                                  " oldLnLikelyhood ="<<oldLnLikelyhood);
181           Real exit:=1
182         },         
183         1,
184         {
185           Text exitMsg = " Logit model iteration("+
186            FormatReal(iter, "%0"+IntText(Floor(Log10(maxIter)+1))+".lf")+
187           ")"+"        "+
188           " LogLikelyhood = "+FormatReal(lnLikelyhood, "%.4E")+"       "+
189           " MaxAbsDif = "+FormatReal(maxAbsDif, "%.4E")+"      "+
190           " Gradient Norm = "+FormatReal(norm, "%.4E")+"       "+
191           " Time = "+FormatReal(difTime, "%.4lf");
192           Real CMsg::Trace::show(1, exitMsg);
193           Real iter := iter+1;   
194           Real oldLnLikelyhood := Copy(lnLikelyhood);
195           Matrix B:=(Matrix B-dif);
196           Real 0           
197         }
198       );       
199       Real exitCycle = If(EQ(exitTreatment, 1),
200       {
201         Set cycle := SetOfAnything
202         (B, p, G, H, dif, norm, advance,
203          maxAbsDif, error, lnLikelyhood, oldLnLikelyhood, ctInfo, X, XIni, y,
204          mLLCorrect);
205         1
206       }, 0)
207     });
208     Real time2     = Copy(Time);
209/*
210     Matrix p       = _Probability(X, B);
211
212     Matrix error   = _Error(y, p);
213
214     Matrix G       = _Gradient(X, error);
215     Matrix H       = _Hessian(vX, p);
216     Matrix dif     = MinimumResidualsSolve(H, G);
217
218     Real norm            = MatFrobeniusNorm(G);
219     Real advance         = MatFrobeniusNorm(dif);
220     Real maxAbsDif       = MatMax(Abs(dif));
221     Matrix mLnLikelyhood = _MLnLikelyhood(y, p);
222     Real lnLikelyhood    = MatSum(mLnLikelyhood);
223*/
224     Real difTime2 = Copy(Time)-time2; 
225     Real totalTime = completeTime+difTime2;
226     Text endMsg =
227      "Model Logit Ended. Time:"<<totalTime+" SplitTime:"<<difTimeMin+NL+NL;
228     Real WriteLn(endMsg); 
229
230/*
231     SetOfAnything
232     (B, p, G, H, dif, norm, advance,
233      maxAbsDif, error, lnLikelyhood, oldLnLikelyhood)
234*/
235     cycle
236   };
237
238   Set Estimate.MaxLikelyhood.Default(Matrix y, Matrix X)
239   {
240     Matrix B0Ini = Rand(Columns(X), 1, 0, 0);
241     Estimate.MaxLikelyhood
242     (
243       y,      // Matriz y de variable salida binaria
244       X,      // Matriz x de variables entrada
245       B0Ini,
246       _.ctInfo,
247       _.dEpsilon, // Diferencia de paso
248       _.maxIter,  // Numero maximo de iteraciones
249       _.tolerance // Tolerancia al error
250     )
251   };
252   Set Estimate.MaxLikelyhood.Constant(Matrix y, Matrix X, Matrix constant)
253   {
254     Matrix B0Ini = Rand(Columns(X), 1, 0, 0);
255     Estimate.MaxLikelyhood
256     (
257       y,      // Matriz y de variable salida binaria
258       X,      // Matriz x de variables entrada
259       B0Ini,
260       constant,
261       _.dEpsilon, // Diferencia de paso
262       _.maxIter,  // Numero maximo de iteraciones
263       _.tolerance // Tolerancia al error
264     )
265   };
266   Set Estimate.MaxLikelyhood.ProbB0
267   (
268     Matrix y,
269     Matrix X,
270     Matrix B0Ini,
271     Real prob
272   )
273   {
274     Matrix probMat  = Rand(Rows(y), 1, prob, prob);
275     Matrix constant = Log(probMat$/RSum(-probMat, 1));
276
277     Estimate.MaxLikelyhood
278     (
279       y,      // Matriz y de variable salida binaria
280       X,      // Matriz x de variables entrada
281       B0Ini,
282       constant,
283       _.dEpsilon, // Diferencia de paso
284       _.maxIter,  // Numero maximo de iteraciones
285       _.tolerance // Tolerancia al error
286     )
287   };
288
289   Set Estimate.MaxLikelyhood.ProbRec
290   (
291     Matrix y,
292     Matrix X,
293     Matrix B0Cur,
294     Real step,
295     Real tolerance,
296     Real probCur,
297     Real probEnd
298   )
299   {
300     Text iniMsg = NL+NL+"Model Logit (p="<<probCur+") "+Time;
301     Real WriteLn(iniMsg); 
302
303     Set logitResult = Estimate.MaxLikelyhood.ProbB0(y, X, B0Cur, probCur);
304     If(LT(Abs(probCur - probEnd), tolerance), logitResult,
305     {
306       Real prob = probCur*step;
307       Matrix B0 = logitResult::B;
308       Estimate.MaxLikelyhood.ProbRec
309       (y, X, B0, step, tolerance, prob, probEnd)
310     })
311   };
312
313   Set Estimate.MaxLikelyhood.Prob
314   (
315     Matrix y,
316     Matrix X,
317     Real prob
318   )
319   {
320     Real step    = (prob/_.probability.Init)^(1/_.num.Step);
321     Matrix B0Ini = Rand(Columns(X), 1, 0, 0);
322     Estimate.MaxLikelyhood.ProbRec
323     (y, X, B0Ini, step, _.tolerance.Rec, _.probability.Init, prob)
324   };
325   
326   Set Diagnosis
327   (
328     Matrix y,
329     Matrix X,
330     Matrix B,
331     Matrix error,
332     Real lnLikelyhood,
333     Set names,
334     Matrix H  //Hessian
335   )
336   {
337     Text iniMsg =
338      "Model Logit Diagnosis. Init Time:"+Time;
339     Real WriteLn(iniMsg);
340
341     Real N = Rows(X);    //
342     Real n = Columns(X); // Parameter number
343
344     Matrix FIM = -H; // Fisher Information Matrix
345     Matrix COV = SVDInverse(FIM); // Varianze Covarianza Parameter Matrix
346     
347     Set Parameters = For(1, Columns(X), Set(Real k)
348     {
349       Text name  = names[k];
350       Real value = MatDat(B, k, 1);
351       Real var   = MatDat(COV, k, k);
352       Real std   = SqRt(var);
353       Real tStudent = value/std;
354       Real refProb  = 2*(1-DistT(Abs(tStudent), N-n-1));
355       ParameterInf
356       (
357         name,     //Name
358         0,        //Factor
359         0,        //Order
360         value,    //Value
361         std,      //StDs
362         tStudent, //TStudent
363         refProb   //RefuseProb
364       )
365     });
366 
367     Real VarTot    = MatVar(y);
368     Real VarError  = MatVar(error);
369     Real R2 = 1-VarError/VarTot;
370
371     Real MaxProb   = MatSum(y)/N;
372     Real lnLikelyhoodIntercep =
373                      MatSum(_MLnLikelyhood(y, Rand(N, 1, MaxProb, MaxProb)));
374
375     Real Nagelkerke.R2    = 1-(Exp((2/N)*(lnLikelyhoodIntercep-lnLikelyhood)));
376     Real Nagelkerke.R2Max = 1-Exp((2/N)*lnLikelyhoodIntercep);
377     Real Nagelkerke.R2MaxRescaled = Nagelkerke.R2/Nagelkerke.R2Max;   
378     Text endMsg =
379      "Model Logit Diagnosis. End Time:"+Time+NL;
380     Real WriteLn(endMsg);
381 
382     SetOfAnything
383     (
384       Parameters,
385       FIM,
386       COV,
387       MaxProb,
388       R2,
389       Nagelkerke.R2,
390       Nagelkerke.R2Max,
391       Nagelkerke.R2MaxRescaled,
392       lnLikelyhood,
393       lnLikelyhoodIntercep
394     )
395   };
396
397  Set PreTesting(Matrix Y, Matrix X, Set varNames)
398  {
399    Text iniMsg =
400     "Model Logit PreTesting. Init Time:"+Time;
401    Real WriteLn(iniMsg);
402 
403    Real WriteLn("  Checking column stability..."+Time);
404    Matrix unkX    = IsUnknown(X);
405    Matrix posInfX = IsPosInf(X);
406    Matrix negInfX = IsNegInf(X);
407 
408    Matrix unkY    = IsUnknown(Y);
409    Matrix posInfY = IsPosInf(Y);
410    Matrix negInfY = IsNegInf(Y);
411 
412    Real isUnkX    = MatSum(unkX);
413    Real isPosInfX = MatSum(posInfX);
414    Real isNegInfX = MatSum(negInfX);
415 
416    Real isUnkY    = MatSum(unkY);
417    Real isPosInfY = MatSum(posInfY);
418    Real isNegInfY = MatSum(negInfY);
419 
420    Real n = Rows(Y);
421    Real balanced = MatSum(Y)/n;
422 
423    Real valid =
424     Not(Or(isUnkX, isPosInfX, isNegInfX, isUnkY, isPosInfY, isNegInfY));
425
426    Set checkValid = If(EQ(valid, 1), Empty,
427    {
428      Set data = For(1, Card(varNames), Set(Real k)
429      {
430        Text name = varNames[k];
431        Real kUnkX    = MatSum(SubCol(unkX, [[k]])); 
432        Real kPosInfX = MatSum(SubCol(posInfX, [[k]]));
433        Real kNegInfX = MatSum(SubCol(negInfX, [[k]]));
434        SetOfAnything(name, kUnkX, kPosInfX, kNegInfX)
435      });
436      Set header = SetOfText("VarName", "Unk", "PosInf", "NegInf");
437      SetOfSet(header)<<
438      SetOfSet(SetOfAnything("Y", isUnkY, isPosInfY, isNegInfY))<<
439      data
440    });
441    Real WriteLn("  Adjusting X matrix..."+Time);
442
443    Matrix YPre =
444     If(isUnkY, IfMat(unkY, VMat2Mat(Eye(Rows(Y), 1, 0, 0)), Y), Y);
445    Matrix XPre =
446     If(isUnkX, IfMat(unkX, VMat2Mat(Eye(Rows(X), Columns(X), 0, 0)), X), X);   
447     
448    VMatrix vXPre = Mat2VMat(XPre);
449    Set index    = Range(1, Columns(X), 1);
450//    Real WriteLn("  Correlation Y|X matrix..."+Time);
451//    Matrix corVarX = Cor(Tra(YPre|XPre));
452
453    Real WriteLn("  X information column..."+Time); 
454    Set ColInfo = For(1, Card(varNames), Set(Real k)
455    {
456      Text name  = varNames[k];
457      Matrix col = SubCol(XPre, [[k]]);
458      Real min   = MatMin(col);
459      Real max   = MatMax(col);
460      Real stds  = MatStDs(col);
461      Real avr   = MatAvr(col);
462      Matrix freq   = Frequency(col, 100, min, max);
463      Matrix freq01 = Frequency(col$*Y, 100, min, max);
464      Matrix ratioDisc = freq01$/freq;
465      Real minValue    = MatDat(freq01, 1, 2);
466      Real maxValue    = MatDat(freq01, 100, 2);
467      Real maxRatio    = MatDat(ratioDisc, 100, 2)/balanced;
468      Real dicotomicRatio = (minValue+maxValue)/n ;
469
470      VMatrix y = SubCol(vXPre, [[k]]);
471      VMatrix x = SubCol(vXPre, index-[[k]]);
472
473      Set linReg   = LinReg::Get.GeneralInformation(y, x);
474      Real R2MultiColinearity = linReg::R2;
475 
476      SetOfAnything
477      (
478        name,
479        min,
480        max,
481        stds,
482        avr,
483        freq,
484        freq01,
485        ratioDisc,
486        maxRatio,
487        maxValue,
488        dicotomicRatio,
489        R2MultiColinearity
490      )
491    });
492
493    Text endMsg =
494     "Model Logit PreTesting. End Time:"+Time+NL;
495    Real WriteLn(endMsg);
496
497    SetOfAnything
498    (unkX, posInfX, negInfX, unkY, posInfY, negInfY, valid, checkValid,
499     balanced, /*corVarX,*/
500     ColInfo, YPre, XPre)
501  }
502]];
503//_.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-._.-.
504