Changeset 946360a in tol


Ignore:
Timestamp:
Jul 2, 2012, 10:15:36 AM (14 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
4fc9c7e
Parents:
c92ed76
Message:

Refs #1579

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@4903 be4846ae-f051-483a-831f-b64e7d08527f

Location:
tol/btol/bmodel
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • tol/btol/bmodel/estim.cpp

    rc92ed76 r946360a  
    188188    M_->IntegrateMissingValues();
    189189    M_->Statistics();
    190     //M_->Diagnostics();
     190    M_->Diagnostics();
    191191  }
    192192  return(ok_);
    … …  
    14571457
    14581458  res = TolLapack::dpotrf(CblasLower,cov,L);
     1459  if(!res) { return(BDat::NegInf()); }
    14591460  res = TolBlas::dtrsm(CblasLeft,CblasLower,CblasNoTrans,CblasNonUnit,alpha,
    14601461                           L,dif,eps);
     1462  if(!res) { return(BDat::NegInf()); }
    14611463  const BDat* eps_i = eps.Data().Buffer();
    14621464  BDat eps_sumSqr = 0;
    … …  
    15621564  int s = T*2;
    15631565  contens_.Alloc(sampleSize,1);
     1566  contens_.GetData().Replicate(BDat::NegInf(),sampleSize);
    15641567  int k;
    15651568  int res;
    … …  
    15701573    _covACF_Bartlett(cov,acf,n,T,s);
    15711574    res = TolLapack::dpotrf(CblasLower,cov,L);
     1575    if(res) { return; }
    15721576  };
    15731577  for(k=0; k<sampleSize; k++)

  • tol/btol/bmodel/modcalc.cpp

    rc92ed76 r946360a  
    20312031  }
    20322032
    2033   //-------------------------------------------------------
    2034   // Autocorrelation Matrix and Box-Pierce-Ljung Statistics
    2035   //-------------------------------------------------------
    2036   aCorNum_ = N/4;
    2037   BData acor(aCorNum_), pacor(aCorNum_);
    2038   {
    2039 //  BTimer tm("Autocorrelation Matrix and Box-Pierce-Ljung Statistics");
    2040     if(periodicity_>0)
    2041     {
    2042       if(aCorNum_< 3*periodicity_) { aCorNum_ =  3*periodicity_; }
    2043       if(aCorNum_>10*periodicity_) { aCorNum_ = 10*periodicity_; }
    2044     }
    2045     else
    2046     {
    2047       if(aCorNum_>100) { aCorNum_ = 100; }
    2048     }
    2049     if(aCorNum_>N/2) { aCorNum_ = N/2; }
    2050     resACor_ = BMat(4,aCorNum_);
    2051     BInt  p = arParam_;
    2052   //BInt  q = maParam_;
    2053     AutoCor(residuous_,acor,aCorNum_);
    2054     BDat Q = 0;
    2055     for(i=0; i<aCorNum_; i++)
    2056     {
    2057       resACor_(0,i)=i+1;
    2058       resACor_(1,i)=acor(i);
    2059       resACor_(2,i)=sigma;
    2060       Q += ((resACor_(1,i)^2)/BDat(N-i-1));
    2061       resACor_(3,i)=Q*BDat(N*(N+2));
    2062     }
    2063   }
    2064 
    2065   //-------------------------------------------------------
    2066   // Partial Autocorrelation Matrix
    2067   //-------------------------------------------------------
    2068   {
    2069 //  BTimer tm("Partial Autocorrelation Matrix");
    2070     resPACor_=BMat(3,aCorNum_);
    2071     PartAutoCor(acor,pacor);
    2072 
    2073     for(i=0; i<aCorNum_; i++)
    2074     {
    2075       resPACor_(0,i)=i+1;
    2076       resPACor_(1,i)=pacor(i);;
    2077       resPACor_(2,i)=sigma;
    2078     }
    2079 
    2080     resACor_  = resACor_ .T();
    2081     resPACor_ = resPACor_.T();
    2082   }
    2083 }
     2033  Std(BText("\nTRACE BModel::Statistics() N_=")+N_);
     2034  aCorNum_ = N_/4;
     2035  Std(BText("\nTRACE BModel::Statistics() aCorNum_=")+aCorNum_);
     2036  if(periodicity_>0)
     2037  {
     2038    if(aCorNum_< 2*periodicity_+4) { aCorNum_ =  2*periodicity_+4; }
     2039    Std(BText("\nTRACE BModel::Statistics() aCorNum_=")+aCorNum_);
     2040  }
     2041  if(aCorNum_<N_/2) { aCorNum_ = N_/2; }
     2042  Std(BText("\nTRACE BModel::Statistics() aCorNum_=")+aCorNum_);
     2043
     2044  if(aCorNum_>0)
     2045  {
     2046    resACor_.Alloc(aCorNum_,4);
     2047    resPACor_.Alloc(aCorNum_,3);
     2048    BData acor, pacor;
     2049    //-------------------------------------------------------
     2050    // Autocorrelation Matrix and Box-Pierce-Ljung Statistics
     2051    //-------------------------------------------------------
     2052    {
     2053    //BTimer tm("Autocorrelation Matrix and Box-Pierce-Ljung Statistics");
     2054      BInt  p = arParam_;
     2055    //BInt  q = maParam_;
     2056    //Std(BText("\nA_.Data()=\n")+Name(A_.Data())+"\n");
     2057      AutoCor(A_.Data(),acor,aCorNum_);
     2058    //Std(BText("\nacor=\n")+Name(acor)+"\n");
     2059      BDat Q = 0;
     2060      for(i=0; i<aCorNum_; i++)
     2061      {
     2062        resACor_(i,0)=i+1;
     2063        resACor_(i,1)=acor(i);
     2064        resACor_(i,2)=sigma;
     2065        Q += ((acor(i)^2)/BDat(N-i-1));
     2066        resACor_(i,3)=Q*BDat(N*(N+2));
     2067      }
     2068    }
     2069
     2070    //-------------------------------------------------------
     2071    // Partial Autocorrelation Matrix
     2072    //-------------------------------------------------------
     2073    {
     2074    //BTimer tm("Partial Autocorrelation Matrix");
     2075      PartAutoCor(acor,pacor);
     2076      for(i=0; i<aCorNum_; i++)
     2077      {
     2078        resPACor_(i,0)=i+1;
     2079        resPACor_(i,1)=pacor(i);;
     2080        resPACor_(i,2)=sigma;
     2081      }
     2082
     2083    }
     2084  }
     2085}

  • tol/btol/bmodel/moddiag.cpp

    rc92ed76 r946360a  
    5959//--------------------------------------------------------------------
    6060{
    61     if(necessary)
    62     {
    63         if(testAccept_(n) <      testRefuse_(n))
    64         {
    65             if     (x<=testAccept_(n)) { diagQualify_(n) = BDIAGGOOD;   }
    66             else if(x<=testRefuse_(n)) { diagQualify_(n) = BDIAGACCEPT; }
    67             else                       { diagQualify_(n) = BDIAGREFUSE; }
    68         }
    69         else
    70         {
    71             if     (x>=testAccept_(n)) { diagQualify_(n) = BDIAGGOOD;   }
    72             else if(x>=testRefuse_(n)) { diagQualify_(n) = BDIAGACCEPT; }
    73             else                       { diagQualify_(n) = BDIAGREFUSE; }
    74         }
    75         diagPunct_(n) = x;
     61  if(necessary)
     62  {
     63    if(testAccept_(n) <   testRefuse_(n))
     64    {
     65      if     (x<=testAccept_(n)) { diagQualify_(n) = BDIAGGOOD;  }
     66      else if(x<=testRefuse_(n)) { diagQualify_(n) = BDIAGACCEPT; }
     67      else                       { diagQualify_(n) = BDIAGREFUSE; }
    7668    }
    7769    else
    7870    {
    79         diagQualify_(n) = BDIAGUNNECESSARY;
    80         diagPunct_       (n) = 0;
    81     }
    82     return(diagQualify_(n));
     71      if     (x>=testAccept_(n)) { diagQualify_(n) = BDIAGGOOD;  }
     72      else if(x>=testRefuse_(n)) { diagQualify_(n) = BDIAGACCEPT; }
     73      else                       { diagQualify_(n) = BDIAGREFUSE; }
     74    }
     75    diagPunct_(n) = x;
     76  }
     77  else
     78  {
     79    diagQualify_(n) = BDIAGUNNECESSARY;
     80    diagPunct_   (n) = 0;
     81  }
     82  return(diagQualify_(n));
    8383}
    8484
    … …  
    8787//--------------------------------------------------------------------
    8888{
    89     return(2*BNormalDist::Dist01(Abs(x)) - BDat(1));
    90 }
     89  return(2*BNormalDist::Dist01(Abs(x)) - BDat(1));
     90}
     91
    9192
    9293//--------------------------------------------------------------------
    … …  
    9798//--------------------------------------------------------------------
    9899{
    99     diagValue_(n) = resACor_.Get(0,1);
    100     BDat y = diagValue_(n)*Sqrt(N_);
    101     BDat x = AbsNormal(y);
    102 //Std(BText("\nFirst regular autocorrelation = ")+y.Name());
    103 //Std(BText("\nProbability                   = ")+x.Name());
    104     return(Qualify(n,x,BTRUE));
     100  Std(BText("\nTRACE FstRegAutCorTest "));
     101  BInt m = 1;
     102  BDat x,y;
     103  if(m<resACor_.Rows()) { y = resACor_(m-1,1)*Sqrt(N_); }
     104  diagValue_(n) = y;
     105  x = AbsNormal(y);
     106  return(Qualify(n,x,BTRUE));
    105107}
    106108
    … …  
    112114//--------------------------------------------------------------------
    113115{
    114     BDat x;
    115     BBool necessary = (periodicity_>1);
    116     if(necessary)
    117     {
    118         diagValue_(n) = resACor_.Get(periodicity_-1,1);
    119         BDat y = diagValue_(n)*Sqrt(N_);
    120         x = AbsNormal(y);
    121 //  Std(BText("\nFirst seasonal autocorrelation = ")+y.Name());
    122 //  Std(BText("\nProbability                    = ")+x.Name());
    123     }
    124     return(Qualify(n,x,necessary));
     116  Std(BText("\nTRACE FstSeaAutCorTest "));
     117  BInt m = periodicity_;
     118  BDat x,y;
     119  BBool necessary = (periodicity_>1);
     120  if(necessary)
     121  {
     122    if(m<resACor_.Rows()) { y = resACor_(m-1,1)*Sqrt(N_); }
     123    diagValue_(n) = y;
     124    x = AbsNormal(y);
     125  }
     126  return(Qualify(n,x,necessary));
    125127}
    126128
    … …  
    132134//--------------------------------------------------------------------
    133135{
    134     diagValue_(n) = resACor_.Get(1,1);
    135     BDat y = diagValue_(n)*Sqrt(N_);
    136     BDat x = AbsNormal(y);
    137 //Std(BText("\nSecond regular autocorrelation = ")+y.Name());
    138 //Std(BText("\nProbability                    = ")+x.Name());
    139     return(Qualify(n,x,BTRUE));
     136  Std(BText("\nTRACE SndRegAutCorTest "));
     137  BInt m = 2;
     138  BDat x,y;
     139  if(m<resACor_.Rows()) { y = resACor_(m-1,1)*Sqrt(N_); }
     140  diagValue_(n) = y;
     141  x = AbsNormal(y);
     142  return(Qualify(n,x,BTRUE));
    140143}
    141144
    … …  
    147150//--------------------------------------------------------------------
    148151{
    149     BDat x;
    150     BBool necessary = (periodicity_>1);
    151     if(necessary)
    152     {
    153         BInt m = Minimum(2*periodicity_-1,resACor_.Rows()-1);
    154         BDat y = resACor_(m,1)*Sqrt(N_);
    155         diagValue_(n) = y;
    156         x = AbsNormal(y);
    157 //  Std(BText("\nSecond seasonal autocorrelation = ")+y.Name());
    158 //  Std(BText("\nProbability                     = ")+x.Name());
    159     }
    160     return(Qualify(n,x,necessary));
     152  Std(BText("\nTRACE SndSeaAutCorTest "));
     153  BInt m = 2*periodicity_;
     154  BDat x,y;
     155  BBool necessary = (periodicity_>1);
     156  if(necessary)
     157  {
     158    if(m<resACor_.Rows()) { y = resACor_(m-1,1)*Sqrt(N_); }
     159    diagValue_(n) = y;
     160    x = AbsNormal(y);
     161  }
     162  return(Qualify(n,x,necessary));
    161163}
    162164
    … …  
    170172//--------------------------------------------------------------------
    171173{
    172     BInt        m = aCorNum_;
    173     BInt        p = arParam_;
    174     BInt        q = maParam_;
    175     RegularBoxPierceLjung_ = BoxPierceLjung(residuous_, m, p, q);
    176     BChiSquareDist Chi(m-p-q);
    177     diagValue_(n) = RegularBoxPierceLjung_;
    178     BDat x =  Chi.Dist(RegularBoxPierceLjung_);
     174  Std(BText("\nTRACE RegBoxPierceLjungTest "));
     175  BInt  m = aCorNum_;
     176  BInt  p = arParam_;
     177  BInt  q = maParam_;
     178  Std(BText("\nTRACE RegBoxPierceLjungTest N=")+A_.Data().Size());
     179  Std(BText("\nTRACE RegBoxPierceLjungTest m=")+m);
     180  Std(BText("\nTRACE RegBoxPierceLjungTest p=")+p);
     181  Std(BText("\nTRACE RegBoxPierceLjungTest q=")+q);
     182  RegularBoxPierceLjung_ = BoxPierceLjung(A_.Data(), m, p, q);
     183  Std(BText("\nTRACE RegBoxPierceLjungTest RegularBoxPierceLjung_=")+RegularBoxPierceLjung_);
     184  BChiSquareDist Chi(m-p-q);
     185  diagValue_(n) = RegularBoxPierceLjung_;
     186  BDat x =  Chi.Dist(RegularBoxPierceLjung_);
     187  Std(BText("\nTRACE RegBoxPierceLjungTest x=")+x);
    179188//Std(BText("\nRegular Box-Pierce-Ljung = ")+RegularBoxPierceLjung_.Name());
    180 //Std(BText("\nProbability              = ")+x.Name());
    181     return(Qualify(n,x,BTRUE));
     189//Std(BText("\nProbability    = ")+x.Name());
     190  return(Qualify(n,x,BTRUE));
    182191}
    183192
    … …  
    190199//--------------------------------------------------------------------
    191200{
    192     BBool necessary = (periodicity_ > 1) &&
    193         (residuous_.Size()/periodicity_ > 2);
    194     BDat x;
    195     if(necessary)
    196     {
    197         BInt  m = aCorNum_/periodicity_;
    198         BInt  p = arFactors_[1].Size()-1;
    199         BInt  q = maFactors_[1].Size()-1;
    200         BArray<BDat> seaRes(N_/periodicity_);
    201 /*
    202   Std(BText("\nm = ")+m);
    203   Std(BText("\np = ")+p);
    204   Std(BText("\nq = ")+q);
    205 */
    206         for(BInt i=0; i<seaRes.Size(); i++)
    207         {
    208             seaRes[i] = residuous_[i*periodicity_];
    209 //    Std(BText("\nseaRes[")+i+"]="+seaRes[i]);
    210         }
    211         SeasonalBoxPierceLjung_ = BoxPierceLjung(seaRes, m, p, q);
    212         diagValue_(n) = SeasonalBoxPierceLjung_;
    213         BChiSquareDist Chi(m-p-q);
    214         x =      Chi.Dist(SeasonalBoxPierceLjung_);
    215 /*
    216   Std(BText("\nSeasonal Box-Pierce-Ljung = ")+
    217   SeasonalBoxPierceLjung_.Name());
    218   Std(BText("\nProbability                 = ")+x.Name());
    219 */
    220     }
    221     return(Qualify(n,x,necessary));
     201  Std(BText("\nTRACE SeaBoxPierceLjungTest "));
     202  BBool necessary = (periodicity_ > 1) &&
     203  (A_.Data().Size()/periodicity_ > 2);
     204  BDat x;
     205  if(necessary)
     206  {
     207    int j = arFactors_.Size()-1;
     208    BInt  m = aCorNum_/periodicity_;
     209    BInt  p = arFactors_[j].Size()-1;
     210    BInt  q = maFactors_[j].Size()-1;
     211    BArray<BDat> seaRes(N_/periodicity_);
     212    for(int i=0; i<seaRes.Size(); i++)
     213    {
     214      seaRes[i] = A_.Data()[i*periodicity_];
     215  //  Std(BText("\nseaRes[")+i+"]="+seaRes[i]);
     216    }
     217    SeasonalBoxPierceLjung_ = BoxPierceLjung(seaRes, m, p, q);
     218    diagValue_(n) = SeasonalBoxPierceLjung_;
     219    BChiSquareDist Chi(m-p-q);
     220    x = Chi.Dist(SeasonalBoxPierceLjung_);
     221  /*
     222    Std(BText("\nSeasonal Box-Pierce-Ljung = ")+
     223    SeasonalBoxPierceLjung_.Name());
     224    Std(BText("\nProbability       = ")+x.Name());
     225  */
     226  }
     227  return(Qualify(n,x,necessary));
    222228}
    223229
    … …  
    229235//--------------------------------------------------------------------
    230236{
     237  Std(BText("\nTRACE KullbackLeiblerDistanceTest "));
    231238  static BReal paso = 0.25;
    232239  static BDat  KullbackLeiblerAvr  = -5.349629861;
    233240  static BDat  KullbackLeiblerStDs =  3*0.544121702;
    234241
    235   BArray<BDat> res01 = residuous_;
    236   BInt         i,j;
    237   BDat         max=0;
     242  BArray<BDat> res01 = A_.Data();
     243  BInt         i,j;
     244  BDat         max=0;
    238245  for(i=0; i<res01.Size(); i++)
    239246  {
    … …  
    255262  for(j=0; j<numInterv; j++)
    256263  {
    257     BReal from  = paso*j;
     264    BReal from  = paso*j;
    258265    BReal until = paso*(j+1);
    259266    BDat pt    = 2*(BNormalDist::Dist01(until)-BNormalDist::Dist01(from));
    … …  
    266273  x = AbsNormal(x);
    267274
    268 //Std(BText("\nNumber of intervals      = ")+numInterv);
     275//Std(BText("\nNumber of intervals  = ")+numInterv);
    269276//Std(BText("\nKullback-Leibler Distance = ")+epsilon.Name());
    270 //Std(BText("\nProbability              = ")+x.Name());
     277//Std(BText("\nProbability    = ")+x.Name());
    271278  return(Qualify(n,x,BTRUE));
    272279}
    … …  
    279286//--------------------------------------------------------------------
    280287{
    281     BInt numInterv = BInt(Sqrt(N_));
     288  Std(BText("\nTRACE FisherDistanceTest "));
     289  BInt numInterv = BInt(Sqrt(N_));
    282290//Std(BText("\nnumInterv=")+numInterv);
    283    
    284     BChiSquareDist Chi(numInterv);
    285     BArray<BDat>        res01  = residuous_;
    286     BMatrix<BDat> freq01;
    287     BInt               i,j;
    288     for(i=0; i<res01.Size(); i++) { res01[i] = res01[i]/standardError_; }
    289     Frequency(res01, freq01, numInterv);
    290     BDat epsilon = 0;
    291     BDat from    = 0;
    292     for(j=0; j<numInterv; j++)
    293     {
    294         BDat until = 1;
    295         if(j!=numInterv-1) { until = BNormalDist::Dist01(freq01(j,0)); }
    296         BDat ft    = res01.Size()*(until-from);
    297         BDat fe    = freq01(j,1);
    298         BDat d     = ((fe-ft)*(fe-ft))/ft;
    299         if(d.IsKnown()) { epsilon += d; }
     291 
     292  BChiSquareDist Chi(numInterv);
     293  BArray<BDat>  res01  = A_.Data();
     294  BMatrix<BDat> freq01;
     295  BInt         i,j;
     296  for(i=0; i<res01.Size(); i++) { res01[i] = res01[i]/standardError_; }
     297  Frequency(res01, freq01, numInterv);
     298  BDat epsilon = 0;
     299  BDat from    = 0;
     300  for(j=0; j<numInterv; j++)
     301  {
     302    BDat until = 1;
     303    if(j!=numInterv-1) { until = BNormalDist::Dist01(freq01(j,0)); }
     304    BDat ft    = res01.Size()*(until-from);
     305    BDat fe    = freq01(j,1);
     306    BDat d     = ((fe-ft)*(fe-ft))/ft;
     307    if(d.IsKnown()) { epsilon += d; }
    300308//  Std(BText("\nj=")+j+"\td="+d.Name());
    301         from  = until;
    302     }
    303     diagValue_(n) = epsilon;
     309    from  = until;
     310  }
     311  diagValue_(n) = epsilon;
    304312//Std(BText("\nepsilon")=epsilon.Name());
    305     BDat x;
    306     if(epsilon.IsKnown()) { x = Chi.Dist(epsilon); }
     313  BDat x;
     314  if(epsilon.IsKnown()) { x = Chi.Dist(epsilon); }
    307315//Std(BText("\nNumber of intervals = ")+numInterv);
    308 //Std(BText("\Fisher Distance      = ")+epsilon.Name());
    309 //Std(BText("\nProbability         = ")+x.Name());
    310     return(Qualify(n,x,BTRUE));
     316//Std(BText("\Fisher Distance     = ")+epsilon.Name());
     317//Std(BText("\nProbability     = ")+x.Name());
     318  return(Qualify(n,x,BTRUE));
    311319}
    312320
    … …  
    319327//--------------------------------------------------------------------
    320328{
    321     static BTStudentDist T(1);
    322     BInt        i;
    323     BDat        min=1000000;
    324     for(i=0; i<param_.Size(); i++)
    325     {
    326         BDat t = Abs(param_[i]/paramSD_[i]);
    327         if(t.IsKnown() && (min > t))
    328         {
    329             min = t;
    330             diagValue_(n) = param_[i]/paramSD_[i];
    331         }
    332     }
    333     BDat x = 2*(1-T.Dist(min));
     329  Std(BText("\nTRACE MinSignificationTest "));
     330  BTStudentDist T(N_-numParam_);
     331  BInt i;
     332  BDat min=BDat::PosInf();
     333  for(i=0; i<numParam_; i++)
     334  {
     335    BDat t = Abs(param_[i]/paramSD_[i]);
     336    if(!t.IsKnown())
     337    {
     338      diagValue_(n) = t;
     339      return(Qualify(n,1,BTRUE));
     340    }
     341    else if(min > t)
     342    {
     343      min = t;
     344      diagValue_(n) = min;
     345    }
     346  }
     347  BDat refuseProb = 2*(1-T.Dist(min));
    334348//Std(BText("\nMinimum signification = ")+min.Name());
    335 //Std(BText("\nProbability         = ")+x.Name());
    336     return(Qualify(n,x,BTRUE));
     349//Std(BText("\nProbability     = ")+x.Name());
     350  return(Qualify(n,refuseProb,BTRUE));
    337351}
    338352
    … …  
    344358//--------------------------------------------------------------------
    345359{
    346     BDat        max=0;
    347     BInt        N = paramCor_.Rows();
    348     for(BInt i=0; i<N; i++)
    349     {
    350         for(BInt j=0; j<i; j++)
    351         {
    352             BDat c = Abs(paramCor_(i,j));
    353             if(max < c)
    354             {
    355                 diagValue_(n) = paramCor_(i,j);
    356                 max = c;
    357             }
    358         }
    359     }
    360 //BDat t = Sqrt(N-2)*max/Sqrt(1-(max^2));
    361 //BTStudentDist T(N-2);
    362 //BDat x = 1-2*T.Dist(Abs(t));
    363 //Std(BText("\nMaximum correlation = ")+max.Name());
    364 //Std(BText("\nProbability         = ")+x.Name());
    365     return(Qualify(n,max,BTRUE));
     360  Std(BText("\nTRACE MaxCorrelationTest "));
     361  BInt  N = paramCor_.Rows();
     362  if(!N || !paramCor_.Columns())
     363  {
     364    diagValue_(n) = BDat::Unknown();
     365    return(Qualify(n,1,BTRUE));
     366  }
     367  BDat  max=0;
     368  for(BInt i=0; i<N; i++)
     369  {
     370    for(BInt j=0; j<i; j++)
     371    {
     372      BDat c = Abs(paramCor_(i,j));
     373      if(max < c)
     374      {
     375        diagValue_(n) = c;
     376        max = c;
     377      }
     378    }
     379  }
     380  return(Qualify(n,max,BTRUE));
    366381}
    367382
    … …  
    373388//--------------------------------------------------------------------
    374389{
    375     static BTStudentDist T(1);
    376     BBool necessary = (diagQualify_[7] != BDIAGREFUSE) &&
    377         (diagQualify_[8] != BDIAGREFUSE);
    378     BDat maxMixed = 0;
    379     if(necessary)
    380     {
    381         for(BInt i=0; i<paramCor_.Rows(); i++)
    382         {
    383             BDat        corr=0;
    384             for(BInt j=0; j<paramCor_.Rows(); j++)
    385             {
    386                 if(i!=j)
    387                 {
    388                     BDat c = Abs(paramCor_(i,j));
    389                     if(corr < c) { corr = c; }
    390                 }
    391             }
    392             BDat sign    = 2*(1-T.Dist(Abs(param_[i]/paramSD_[i])));
    393             BDat mixed = Sqrt(Abs(sign*corr));
    394             if(maxMixed<mixed) { maxMixed = mixed; }
    395         }
    396         diagValue_(n) = maxMixed;
     390  Std(BText("\nTRACE MixedSignCorrTest "));
     391  BTStudentDist T(N_-numParam_);
     392  BBool necessary = (diagQualify_[7] != BDIAGREFUSE) &&
     393                    (diagQualify_[8] != BDIAGREFUSE);
     394  BDat maxMixed = 0;
     395  if(necessary)
     396  {
     397    if(!paramCor_.Rows() || !paramCor_.Columns())
     398    {
     399      diagValue_(n) = BDat::Unknown();
     400      return(Qualify(n,1,BTRUE));
     401    }
     402    for(BInt i=0; i<paramCor_.Rows(); i++)
     403    {
     404      BDat  corr=0;
     405      for(BInt j=0; j<paramCor_.Rows(); j++)
     406      {
     407        if(i!=j)
     408        {
     409          BDat c = Abs(paramCor_(i,j));
     410          if(corr < c) { corr = c; }
     411        }
     412      }
     413      BDat sign   = 2*(1-T.Dist(Abs(param_[i]/paramSD_[i])));
     414      BDat mixed = Sqrt(Abs(sign*corr));
     415      if(maxMixed<mixed) { maxMixed = mixed; }
     416    }
     417    diagValue_(n) = maxMixed;
    397418//  Std(BText("\nMaximum mixed correlation-signification = ")+
    398 //      maxMixed.Name());
    399     }
    400     return(Qualify(n,maxMixed,necessary));
    401 }
     419//  maxMixed.Name());
     420  }
     421  return(Qualify(n,maxMixed,necessary));
     422}
     423
    402424
    403425
    … …  
    409431//--------------------------------------------------------------------
    410432{
    411 //Std("\nBModel::UnitRootsProbTest");
    412   BDat                  x = 0;
    413   BInt                  num = arParam_+maParam_;
    414   BBool necessary = 0; /*num!=0; */
     433  Std(BText("\nTRACE UnitRootsProbTest "));
     434  BDat x = 0;
     435  BInt num = arParam_+maParam_;
     436  BBool necessary = num!=0;
    415437  if(necessary)
    416438  {
     439    if(!paramCor_.Rows() || !paramCor_.Columns())
     440    {
     441      diagValue_(n) = BDat::Unknown();
     442      return(Qualify(n,1,BTRUE));
     443    }
    417444//  Std(BText("\nARMA Parameters Number = ")+num);
    418     BInt                   m   = paramCor_.Rows() - num;
     445    int facNum = arFactors_.Size();
     446    BInt m = paramCor_.Rows() - num;
    419447//  Std(BText("\nNON ARMA Parameters Number = ")+m);
    420     BInt                   i, j, k, iter;
    421     BSymMatrix<BDat>       cov (num);
    422     BMatrix<BDat>          param(num,1);
    423     BArray<BDat>           unarys(4);
    424     BArray<BPol>           factors(4);
    425 
    426     factors[0] = arFactors_[0];
    427     factors[1] = arFactors_[1];
    428     factors[2] = maFactors_[0];
    429     factors[3] = maFactors_[1];
    430     unarys [0] = 0;
    431     unarys [1] = 0;
    432     unarys [2] = 0;
    433     unarys [3] = 0;
     448    BInt i, j, k, iter;
     449    BSymMatrix<BDat> cov (num);
     450    BMatrix<BDat> param(num,1);
     451    BArray<BDat> unarys(2*facNum);
     452    BArray<BPol> factors(2*facNum);
     453
     454    for(j=0;j<facNum;j++)
     455    {
     456      factors[j]        = arFactors_[j];
     457      factors[j+facNum] = maFactors_[j];
     458      unarys [j]        = 0;
     459      unarys [j+facNum] = 0;
     460    };
    434461
    435462    i=0;
    436463//  Std(BText("\nFactors =\n"));
    437     for(i=j=0; j<4; j++)
     464    for(i=j=0; j<2*facNum; j++)
    438465    {
    439466//    Std(factors[j].Name());
    440467      for(k=1; k<factors[j].Size(); k++)
    441468      {
    442         param(i++,0)=-factors[j][k].Coef();
     469        param(i++,0)=-factors[j][k].Coef();
    443470      }
    444471    }
    445 
    446 
    447472//  Std(BText("\nParam =\n")+param.Name());
    448 
    449473    for(i=0; i<num; i++)
    450474    {
    … …  
    452476      for(j=0; j<i; j++)
    453477      {
    454         cov(i,j) = paramSD_(m+i)*paramSD_(m+j)*paramCor_.Get(m+i,m+j);
     478        cov(i,j) = paramSD_(m+i)*paramSD_(m+j)*paramCor_.Get(m+i,m+j);
    455479      }
    456480    }
    … …  
    459483    BMultivarDist dist(&n01,param,cov);
    460484
    461     BInt maxIter = 100;
     485    BInt maxIter = 1000;
    462486    for(iter=1; iter<=maxIter; iter++)
    463487    {
    464488      BMatrix<BDat> simula = dist.Random();
    465       for(i=j=0; j<4; j++)
     489      for(i=j=0; j<2*facNum; j++)
    466490      {
    467         for(k=1; k<factors[j].Size(); k++,i++)
    468         { factors[j][k].PutCoef(-simula(0,i)); }
    469 //        Std(BText("\nFactor[")+j+"]"+factors[j].Name() +
    470 //                " Stationary = " + factors[j].Stationary());
    471         if(!(IsStationary(factors[j]))) { unarys[j]+=1; }
     491        for(k=1; k<factors[j].Size(); k++,i++)
     492        { factors[j][k].PutCoef(-simula(0,i)); }
     493        if(!(IsStationary(factors[j]))) { unarys[j]+=1; }
    472494      }
    473495    }
    474496    x = 0;
    475     for(j=0; j<4; j++)
     497    for(j=0; j<2*facNum; j++)
    476498    {
    477499  //  Std(BText("unarys[")+j+"] = " + unarys[j]);
    … …  
    499521{
    500522//Std(Out()+"\n\nRunning test "+n+" " +testTitle_(n));
    501     switch(n)
    502     {
    503         case  0 : return(FstRegAutCorTest          (n));
    504         case  1 : return(FstSeaAutCorTest          (n));
    505         case  2 : return(SndRegAutCorTest          (n));
    506         case  3 : return(SndSeaAutCorTest          (n));
    507         case  4 : return(RegBoxPierceLjungTest     (n));
    508         case  5 : return(SeaBoxPierceLjungTest     (n));
    509         case  6 : return(FisherDistanceTest        (n));
    510         case  7 : return(MinSignificationTest      (n));
    511         case  8 : return(MaxCorrelationTest        (n));
    512         case  9 : return(MixedSignCorrTest         (n));
    513         case 10 : return(UnitRootsProbTest         (n));
    514         default : return(-1);
    515     }
     523  switch(n)
     524  {
     525    case  0 : return(FstRegAutCorTest     (n));
     526    case  1 : return(FstSeaAutCorTest     (n));
     527    case  2 : return(SndRegAutCorTest     (n));
     528    case  3 : return(SndSeaAutCorTest     (n));
     529    case  4 : return(RegBoxPierceLjungTest     (n));
     530    case  5 : return(SeaBoxPierceLjungTest     (n));
     531    case  6 : return(FisherDistanceTest     (n));
     532    case  7 : return(MinSignificationTest      (n));
     533    case  8 : return(MaxCorrelationTest     (n));
     534    case  9 : return(MixedSignCorrTest     (n));
     535    case 10 : return(UnitRootsProbTest     (n));
     536    default : return(-1);
     537  }
    516538}
    517539
    … …  
    533555  {
    534556    Error(I2("Cannot find Set DiagnosticsBounds.",
    535              "No se encuentra el conjunto DiagnosticsBounds."));
     557             "No se encuentra el conjunto DiagnosticsBounds."));
    536558    return;
    537559  }
    … …  
    540562  {
    541563    Error(I2("Wrong number of elements in Set DiagnosticsBounds.",
    542              "Error en el número de elementos del conjunto "
    543              "DiagnosticsBounds."));
     564             "Error en el número de elementos del conjunto "
     565             "DiagnosticsBounds."));
    544566    return;
    545567  }
    … …  
    549571  {
    550572    BSet& s = Set(set[n+1]);
    551     diagValue_  (n) = BDat::Unknown();
    552     diagPunct_  (n) = 0;
     573    diagValue_  (n) = BDat::Unknown();
     574    diagPunct_  (n) = 0;
    553575    testAccept_ (n) = Dat(s[1]);
    554576    testRefuse_ (n) = Dat(s[2]);
    555577    diagQualify_(n) = RunTest(n);
    556578    if(qualification_<diagQualify_(n)) { qualification_=diagQualify_(n); }
    557     BDat x0 = BDat(1)/3;
    558     BDat y0 = BDat(2)/3;
     579    BDat x0 = 1.0/3.0;
     580    BDat y0 = 2.0/3.0;
    559581    BDat x  = testAccept_ (n);
    560582    BDat y  = testRefuse_ (n);
    561     BDat a  = (y0-x0)/(y-x);
    562     BDat b  = x0 - a*x;
    563     BDat q  = a*diagPunct_(n)+b;
     583    BDat p  = diagPunct_(n);
     584    BDat q;
     585         if(p<=x) { q = p * x0/x; }
     586    else if(p<=y) { q = x0 + (p-x) * (y0-x0)/(y-x); }
     587    else          { q = y0 + (p-y) * ( 1-y0)/(1-y); }
    564588    arithmeticQualification_ += q;
    565589/*
    566     Std( BText("\nBounds[")     + testAccept_ (n) +
    567                ", "             + testRefuse_ (n) +
    568                "]; Value   = "  + diagValue_  (n) +
    569                "; Prob   = "    + diagPunct_  (n) +
    570                "; Qualify = "   + qualifTitle_[diagQualify_ (n)] +";"+
    571                "; Qualify(1/3,2/3) = " + q.Name() +";\n");
     590    Std( BText("\nBounds[")  + testAccept_ (n) +
     591         ", "    + testRefuse_ (n) +
     592         "]; Value   = "  + diagValue_  (n) +
     593         "; Prob   = "  + diagPunct_  (n) +
     594         "; Qualify = "  + qualifTitle_[diagQualify_ (n)] +";"+
     595         "; Qualify(1/3,2/3) = " + q.Name() +";\n");
    572596*/
    573597  }
    574598  if(periodicity_>1) { arithmeticQualification_ /= testTitle_.Size(); }
    575   else               { arithmeticQualification_ /= (testTitle_.Size()-3); }
     599  else               { arithmeticQualification_ /= (testTitle_.Size()-3); }
    576600/*
    577601  Std(BText("\nModel Average Qualify(1/3, 2/3) = ") +
    578             arithmeticQualification_.Name() +";\n");
     602      arithmeticQualification_.Name() +";\n");
    579603*/
    580604}

  • tol/btol/bmodel/model.cpp

    rc92ed76 r946360a  
    150150                                  "Distancia a la normal de Fisher");
    151151    BModel::testTitle_[ 7] = I2("Minimum Signification of Parameters",
    152                                 "Máxima Significación de los Parámetros");
     152                                "Mínima Significación de los Parámetros");
    153153    BModel::testTitle_[ 8] = I2("Maximum Correlation of Parameters",
    154154                                "Máxima Correlación de los Parámetros");
    … …  
    172172     nonLinFun_(NIL),
    173173     nonLinPar_(),
    174      DI_(NULL)
     174     DI_(NULL),
     175     aCorNum_(0)
    175176{
    176177    ModelInitialize();
    … …  
    188189    testAccept_ .ReallocBuffer(testTitle_.Size());
    189190    testRefuse_ .ReallocBuffer(testTitle_.Size());
     191   
    190192}
    191193
    … …  
    208210: aborted_(BFALSE),
    209211  constant_(BDat::Unknown()),
    210   DI_(NULL)
     212  DI_(NULL),
     213  aCorNum_(0)
    211214{
    212215  inData_ .AddList  (inputSeries);
    … …  
    245248BModel::BModel(BSet& set)
    246249//--------------------------------------------------------------------
    247     : aborted_(BFALSE), constant_(BDat::Unknown()), DI_(NULL)
     250    : aborted_(BFALSE), constant_(BDat::Unknown()), DI_(NULL),
     251     aCorNum_(0)
    248252{
    249253    ModelInitialize();
    … …  
    758762  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,DataNumber,N_));
    759763  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,VarNumber,param_.Size()));
    760   LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Average,Average(residuous_)));
     764  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Average,Average(A_.Data())));
    761765  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,RSS,rss));
    762766  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Variance,var));
    763767  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Sigma,standardError_));
    764   LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Asymmetry,AsymmetryCoefficient(residuous_)));
    765   LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Kurtosis,Kurtosis(residuous_)));
     768  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Asymmetry,AsymmetryCoefficient(A_.Data())));
     769  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Kurtosis,Kurtosis(A_.Data())));
    766770  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,Pearson,pearson));
    767771  LstFastAppend(lsta,auxa,tmp=NewFND(BDat,R2,r2));
Note: See TracChangeset for help on using the changeset viewer.