stressGoFTest.cxx

Go to the documentation of this file.

#include "TBenchmark.h"
#include "Math/GoFTest.h"
#include "Math/Functor.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TStopwatch.h"

#include "Math/GSLRndmEngines.h"
#if not defined(__CINT__)
#include "Math/Random.h"
#endif

#include "TRandom3.h"

#include <iostream>

#include <cassert>

/*N.B.: The tests' expected values (expectedDn and expectedA2) were computed on Pcphsft54.cern.ch i386 GNU/Linux computer (slc4_ia32_gcc34)
*/
struct GoFTStress {

   static enum EDebugLevelTypes {
      kNoDebug,       // no test results printing
      kBasicDebug,    // prints either "OK" or "FAIL" test results
      kStandardDebug, // prints summarized test results
   } fgDebugLevel;
   
   Int_t RunTests() {
      Int_t result = 0;
      result += UnitTest1();
      result += UnitTest2();
      result += UnitTest3();
      result += UnitTest4();
      result += UnitTest5();
      result += UnitTest6();
      result += UnitTest7();
      return result;
   }
      
   void PrintBenchmark() {
      if (fgDebugLevel == kNoDebug) {
         return;
      }
         
         //Print table with results
      Bool_t UNIX = strcmp(gSystem->GetName(), "Unix") == 0;
      printf("******************************************************************\n");
      if (UNIX) {
         FILE *fp = gSystem->OpenPipe("uname -a", "r");
         Char_t line[60];
         fgets(line,60,fp); line[59] = 0;
         printf("*  SYS: %s\n",line);
         gSystem->ClosePipe(fp);
      } else {
         const Char_t *os = gSystem->Getenv("OS");
         if (!os) {
         printf("*  SYS: Windows 95\n");
         } else {
         printf("*  SYS: %s %s \n",os,gSystem->Getenv("PROCESSOR_IDENTIFIER"));
         }
      }
         
      printf("*****************************************************************************************\n");
      gBenchmark->Print("GoFTestStress");
#ifdef __CINT__
      Double_t reftime = 0.02 ; // slc4rabacal interpreted (CPU time taken to run complete tests with ACLiC)
#else
      Double_t reftime = 0.04; // slc4rabacal compiled (CPU time taken to run complete tests)
#endif
      
      Double_t rootmarks = 800. * reftime / gBenchmark->GetCpuTime("GoFTestStress");
      
      printf("*****************************************************************************************\n");
      printf("*  ROOTMARKS = %6.1f   *  Root%-8s  %d/%d\n", rootmarks, gROOT->GetVersion(),
            gROOT->GetVersionDate(),gROOT->GetVersionTime());
      printf("*****************************************************************************************\n");
   }
   
   private:  
   
   Int_t UnitTest1() {
      std::cout << "UNIT TEST 1" << std::endl;
      
      UInt_t nsmps = 2;
      const UInt_t smpSize1 = 71;
      const UInt_t smpSize2 = 135;
   /*
      Data set adapted from the paper (1)
      "K-Sample Anderson-Darling Tests of Fit for continuous and discrete cases" by Scholz and Stephens
   */
      const Double_t smp1[smpSize1] = {194, 15, 41, 29, 33, 181, 413, 14, 58, 37, 100, 65, 9, 169, 447, 18, 4, 36, 201, 118, 34, 31, 18, 18, 67, 57, 62, 7, 22, 34, 90, 10, 60, 186, 61, 49, 14, 24, 56, 20, 79, 84, 44, 59, 29, 118, 25, 156, 310, 76, 26, 44, 23, 62, 130, 208, 70, 101, 208, 74, 57, 48, 29, 502, 12, 70, 21, 29, 386, 59, 27};
   
      const Double_t smp2[smpSize2] = {55, 320, 56, 104, 220, 239, 47, 246, 176, 182, 33, 23, 261, 87, 7, 120, 14, 62, 47, 225, 71, 246, 21, 42, 20, 5, 12, 120, 11, 3, 14, 71, 11, 14, 11, 16, 90, 1, 16, 52, 95, 97, 51, 11, 4, 141, 18, 142, 68, 77, 80, 1, 16, 106, 206, 82, 54, 31, 216, 46, 111, 39, 63, 18, 191, 18, 163, 24, 50, 44, 102, 72, 22, 39, 3, 15, 197, 188, 79, 88, 46, 5, 5, 36, 22, 139, 210, 97, 30, 23, 13, 14, 359, 9, 12, 270, 603, 3, 104, 2, 438, 50, 254, 5, 283, 35, 12, 487, 18, 100, 7, 98, 5, 85, 91, 43, 230, 3, 130, 102, 209, 14, 57, 54, 32, 67, 59, 134, 152, 27, 14, 230, 66, 61, 34};
   
      ROOT::Math::GoFTest* goft = new ROOT::Math::GoFTest( smpSize1, smp1, smpSize2, smp2 );
      
      Double_t A2 = goft->AndersonDarling2SamplesTest("t"); // standartized A2_akN
   //     Double_t pvalueAD = goft->AndersonDarling2SamplesTest();
      Double_t pvalueAD = (*goft)(ROOT::Math::GoFTest::kAD2s);
   
      Double_t expectedA2_akN = 1.58334 ; // A2_akN in (1)
   
      Double_t sigmaN = 0.754539;    // sigmaN in (1)
   
      Double_t zScore = 0.773108;    // zScore in (1)
   
      Int_t result = PrintResultAD2Samples(nsmps, A2, expectedA2_akN, sigmaN, zScore, pvalueAD);
   
   //     Double_t Dn = goft->KolmogorovSmirnov2SamplesTest("t"); // standartized A2_akN
      Double_t Dn = (*goft)(ROOT::Math::GoFTest::kKS2s, "t");
      Double_t pvalueKS = goft->KolmogorovSmirnov2SamplesTest();
   
      Double_t expectedDn = 0.139176;
   
      result += PrintResultKS(nsmps, Dn, expectedDn, pvalueKS);
   
      delete goft;
      return result;
   }
   
   Int_t UnitTest2() {
      std::cout << "UNIT TEST 2" << std::endl;
      
      const UInt_t nsmps = 2;
      const UInt_t smpSize = 16;
   /*
      Data sets adapted from the paper (1)
      "K-Sample Anderson-Darling Tests of Fit for continuous and discrete cases" by Scholz and Stephens
   */
      const Double_t smp1[smpSize] = {38.7, 41.5, 43.8, 44.5, 45.5, 46.0, 47.7, 58.0, 39.2, 39.3, 39.7, 41.4, 41.8, 42.9, 43.3, 45.8};
      const Double_t smp2[smpSize] = {34.0, 35.0, 39.0, 40.0, 43.0, 43.0, 44.0, 45.0, 34.0, 34.8, 34.8, 35.4, 37.2, 37.8, 41.2, 42.8};
   
      ROOT::Math::GoFTest* goft = new ROOT::Math::GoFTest( smpSize, smp1, smpSize, smp2 );
   
   //     Double_t A2 = goft->AndersonDarling2SamplesTest("t"); // standartized A2_akN
      Double_t A2 = (*goft)(ROOT::Math::GoFTest::kAD2s, "t");
      Double_t pvalueAD = goft->AndersonDarling2SamplesTest();
   
      Double_t expectedA2_akN = 4.5735; // unstandardized A2_akN in (1)
   
      Double_t sigmaN = 0.719388;   // sigmaN in (1)
   
      Double_t zScore = 4.96748;    // zScore in (1)
   
      Int_t result = PrintResultAD2Samples(nsmps, A2, expectedA2_akN, sigmaN, zScore, pvalueAD);
   
      Double_t Dn = goft->KolmogorovSmirnov2SamplesTest("t"); // standartized A2_akN
   //     Double_t pvalueKS = goft->KolmogorovSmirnov2SamplesTest();
      Double_t pvalueKS = (*goft)(ROOT::Math::GoFTest::kKS2s);
   
      Double_t expectedDn = 0.5;
   
      result += PrintResultKS(nsmps, Dn, expectedDn, pvalueKS);
   
      delete goft;
      return result;
   }
   
   Int_t UnitTest3() {
      std::cout << "UNIT TEST 3" << std::endl;
      
      UInt_t nEvents = 1000;
      UInt_t nsmps = 1;
         
      ROOT::Math::Random<ROOT::Math::GSLRngMT> r;
   
      Double_t* sample = new Double_t[nEvents];
   
      for (UInt_t i = 0; i < nEvents; ++i) { 
         Double_t data = r.LogNormal(5.0, 2.0);
         sample[i] = data;
         assert(sample[i] == data);
      }
      
      ROOT::Math::GoFTest* goft = new ROOT::Math::GoFTest(nEvents, sample, ROOT::Math::GoFTest::kLogNormal);
      
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug)
         std::cout << "LogNormal fitting" << std::endl;
      
   //     Double_t A2 = goft->AndersonDarlingTest("t");
      Double_t A2 = goft->operator()(ROOT::Math::GoFTest::kAD, "t");
      Double_t pvalueAD = goft->AndersonDarlingTest();
   
      Double_t expectedA2 = 0.422771;
   
      Int_t result = PrintResultAD1Sample(A2, expectedA2, pvalueAD);
   
      Double_t Dn = goft->KolmogorovSmirnovTest("t");
   //     Double_t pvalueKS = goft->KolmogorovSmirnovTest();
      Double_t pvalueKS = goft->operator()(ROOT::Math::GoFTest::kKS);
   
      Double_t expectedDn = 0.0204916;
   
      result += PrintResultKS(nsmps, Dn, expectedDn, pvalueKS);
   
      delete goft;
      return result;
   }
   
   Int_t UnitTest4() {
      std::cout << "UNIT TEST 4" << std::endl;
      
      UInt_t nEvents = 1000;
      UInt_t nsmps = 1;
      
      TRandom3 r;
   
      Double_t* sample = new Double_t[nEvents];
   
      for (UInt_t i = 0; i < nEvents; ++i) { 
         Double_t data = r.Exp(1.54);
         sample[i] = data;
         assert(sample[i] == data);
      }
      
      ROOT::Math::GoFTest* goft = new ROOT::Math::GoFTest(nEvents, sample, ROOT::Math::GoFTest::kExponential);
   
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug)
         std::cout << "**Exponential fitting**" << std::endl;
      
      Double_t A2 = goft->AndersonDarlingTest("t");
   //     Double_t pvalueAD = goft->AndersonDarlingTest();
      Double_t pvalueAD = goft->operator()();
   
      Double_t expectedA2 = 0.521153;
   
      Int_t result = PrintResultAD1Sample(A2, expectedA2, pvalueAD);
   
   //     Double_t Dn = goft->KolmogorovSmirnovTest("t");
      Double_t Dn = goft->operator()(ROOT::Math::GoFTest::kKS, "t");
      Double_t pvalueKS = goft->KolmogorovSmirnovTest();
   
      Double_t expectedDn = 0.0218148;
   
      result += PrintResultKS(nsmps, Dn, expectedDn, pvalueKS);
   
      delete goft;
      return result;
   }
   
   Int_t UnitTest5() {
      std::cout << "UNIT TEST 5" << std::endl;
      
      UInt_t nEvents = 1000;
      UInt_t nsmps = 1;
      
      TRandom3 r;
   
      Double_t* sample = new Double_t[nEvents];
   
      for (UInt_t i = 0; i < nEvents; ++i) { 
         Double_t data = r.Gaus(300, 50);
         sample[i] = data;
         assert(sample[i] == data);
      }
      
      ROOT::Math::GoFTest* goft = new ROOT::Math::GoFTest(nEvents, sample);
      goft->SetDistribution(ROOT::Math::GoFTest::kGaussian);
   
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug)
         std::cout << "**Gaussian fitting**" << std::endl;
      
      Double_t A2 =  goft->AndersonDarlingTest("t");
      Double_t pvalueAD = goft->AndersonDarlingTest();
   
      Double_t expectedA2 = 0.441755;
   
      Int_t result = PrintResultAD1Sample(A2, expectedA2, pvalueAD);
   
      Double_t Dn = goft->KolmogorovSmirnovTest("t");
      Double_t pvalueKS = goft->KolmogorovSmirnovTest();
   
      Double_t expectedDn = 0.0282508;
   
      result += PrintResultKS(nsmps, Dn, expectedDn, pvalueKS);
      
      delete goft;
      return result;
   }
   
   Int_t UnitTest6() {
      std::cout << "UNIT TEST 6" << std::endl;
      UInt_t nEvents = 1000;
      UInt_t nsmps = 1;
      
      TRandom3 r;
   
      Double_t* sample = new Double_t[nEvents];
   
      for (UInt_t i = 0; i < nEvents; ++i) { 
         Double_t data = r.Landau();
         sample[i] = data;
         assert(sample[i] == data);
      }
      
      ROOT::Math::Functor1D userCdf(&TMath::LandauI);
      ROOT::Math::GoFTest* goft = new ROOT::Math::GoFTest(nEvents, sample, userCdf, ROOT::Math::GoFTest::kCDF);
         
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug)
         std::cout << "**Landau fitting**" << std::endl;
      
      Double_t A2 =  goft->AndersonDarlingTest("t");
      Double_t pvalueAD = goft->AndersonDarlingTest();
   
      Double_t expectedA2 = 0.544658;
   
      Int_t result = PrintResultAD1Sample(A2, expectedA2, pvalueAD);
      
      Double_t Dn = goft->KolmogorovSmirnovTest("t");
      Double_t pvalueKS = goft->KolmogorovSmirnovTest();
   
      Double_t expectedDn = 0.0203432;
   
      result += PrintResultKS(nsmps, Dn, expectedDn, pvalueKS);
      
      delete goft;
      return result;
   }
   
   Int_t UnitTest7() {
      std::cout << "UNIT TEST 7" << std::endl;
      
      UInt_t nEvents = 1000;
      UInt_t nsmps = 1;
      
      TRandom3 r;
   
      Double_t* sample = new Double_t[nEvents];
   
      for (UInt_t i = 0; i < nEvents; ++i) { 
         Double_t data = r.Landau();
         sample[i] = data;
      }
      
      // need to specify min and max otherwise pdf does not converge
      ROOT::Math::GoFTest* goft = new ROOT::Math::GoFTest(nEvents, sample); 

      ROOT::Math::Functor1D userPdf(&TMath::Landau);
      // need to use a reasanble range for the Landau 
      // but must be bigger than xmin and xmax 
      double xmin = 3*TMath::MinElement(nEvents, sample);
      double xmax = 3*TMath::MaxElement(nEvents, sample);
         
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug)
         std::cout << "**Landau fitting**" << " in [ " << xmin << " , " << xmax << " ]" << std::endl;

      goft->SetUserDistribution(userPdf,ROOT::Math::GoFTest::kPDF, xmin, xmax);
      
      Double_t A2 =  goft->AndersonDarlingTest("t");
      Double_t pvalueAD = goft->AndersonDarlingTest();
   
      Double_t expectedA2 =  0.544658;
   
      // use larger tolerance due to truncation error of Landau
      Int_t result = PrintResultAD1Sample(A2, expectedA2, pvalueAD,0.002);
      
      Double_t Dn = goft->KolmogorovSmirnovTest("t");
      Double_t pvalueKS = goft->KolmogorovSmirnovTest();
   
      Double_t expectedDn = 0.0203432;
   
      result += PrintResultKS(nsmps, Dn, expectedDn, pvalueKS,0.001);
      
      delete goft;
      return result;
   }
   
   Int_t PrintResultAD2Samples(UInt_t nsmps, Double_t A2, Double_t expectedA2_akN, Double_t sigmaN, Double_t zScore, Double_t pvalue, Double_t limit = 0.0001) {
   
      Double_t A2_akN = A2  * sigmaN + (nsmps - 1);
   
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug) {
         std::cout << "Anderson-Darling 2-Samples Test" << std::endl;
   
         std::cout << "pvalue = " << pvalue << std::endl;
   
         std::cout << "zScore = " << zScore << std::endl;
      
         std::cout << "Expected A2 value = " << A2 << std::endl;
      
         std::cout << "A2_akN = " << A2_akN << std::endl;
      
         std::cout << "Expected A2_akN value = " << expectedA2_akN << std::endl;
      }
      if (TMath::Abs(A2_akN - expectedA2_akN) < limit * expectedA2_akN || TMath::Abs(A2 - zScore) < limit * zScore) {
         if (GoFTStress::fgDebugLevel >= GoFTStress::kBasicDebug) 
         std::cout << "+++++++ TEST OK +++++++ \n" << std::endl;
         return EXIT_SUCCESS;
      } else {
         if (GoFTStress::fgDebugLevel >= GoFTStress::kBasicDebug) 
         std::cout << "------- TEST FAIL ------- \n" << std::endl;
         return EXIT_FAILURE;
      }
   }
   
   Int_t PrintResultAD1Sample(Double_t A2, Double_t expectedA2, Double_t pvalue, Double_t limit = 0.0001) {
      
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug) {  
         std::cout << "Anderson-Darling 1-Sample Test" << std::endl;
         
         std::cout << "pvalue = " << pvalue << std::endl;
      
         std::cout << "A2 value = " << A2 << std::endl;
      
         std::cout << "Expected A2 value = " << expectedA2 << std::endl;
      }
      if (TMath::Abs(A2 - expectedA2) < limit * expectedA2) {
         if (GoFTStress::fgDebugLevel >= GoFTStress::kBasicDebug) 
         std::cout << "+++++++ TEST OK +++++++ \n" << std::endl;
         return EXIT_SUCCESS;
      } else {
         if (GoFTStress::fgDebugLevel >= GoFTStress::kBasicDebug) 
         std::cout << "------- TEST FAIL ------- \n" << std::endl;
         return EXIT_FAILURE;
      }
   }
   
   Int_t PrintResultKS(UInt_t nsmps, Double_t Dn, Double_t expectedDn, Double_t pvalue, Double_t limit = 0.0001) {
   
      if (GoFTStress::fgDebugLevel == GoFTStress::kStandardDebug) {
         std::cout << "Kolmogorov-Smirnov "<< nsmps << "-Samples Test" << std::endl;
         
         std::cout << "pvalue = " << pvalue << std::endl;
      
         std::cout << "Dn value = " << Dn << std::endl;
      
         std::cout << "Expected Dn value = " << expectedDn << std::endl;
      }
      if (TMath::Abs(Dn - expectedDn) < limit * expectedDn) {
         if (GoFTStress::fgDebugLevel >= GoFTStress::kBasicDebug)
         std::cout << "+++++++ TEST OK +++++++ \n" << std::endl;
         return EXIT_SUCCESS;
      } else {
         if (GoFTStress::fgDebugLevel >= GoFTStress::kBasicDebug) 
         std::cout << "------- TEST FAIL ------- \n" << std::endl;
         return EXIT_FAILURE;
      }
   }
};
   
#ifdef __MAKECINT__
#pragma link C++ class GoFTStress-;
#endif
   
GoFTStress::EDebugLevelTypes GoFTStress::fgDebugLevel;

Int_t RunTests(Int_t argc, Char_t* argv[]) {
   Int_t result = 0;
   Bool_t validOption(kFALSE);
   
   if (argc == 1) {
      GoFTStress::fgDebugLevel = GoFTStress::kBasicDebug;
      validOption = kTRUE;
   } else if (argc == 3 && strcmp(argv[1], "--o") == 0) {
      if (strcmp(argv[2], "0") == 0 || strcmp(argv[2], "no-debug") == 0) {
         GoFTStress::fgDebugLevel = GoFTStress::kNoDebug;
         validOption = kTRUE;
      } else if (strcmp(argv[2], "1") == 0 || strcmp(argv[2], "basic-debug") == 0) {
         GoFTStress::fgDebugLevel = GoFTStress::kBasicDebug;
         validOption = kTRUE;
         std::cout << "Debug level: basic." << std::endl;
      } else if (strcmp(argv[2], "2") == 0 || strcmp(argv[2], "standard-debug") == 0) {
         GoFTStress::fgDebugLevel = GoFTStress::kStandardDebug;
         validOption = kTRUE;
         std::cout << "Debug level: standard." << std::endl;
      } else if (strcmp(argv[2], "help") == 0) {
         std::cout << "Showing help: " << std::endl;
      }
   }
   
   if (!validOption) {
      std::cout << "Please type just one of the following debug levels in either formats preceded by \"--o\":" << std::endl;
      std::cout << "a) no-debug(0) [no test results printing]" << std::endl;
      std::cout << "b) basic-debug(1) [prints either \"OK\" or \"FAIL\" test results]" << std::endl;
      std::cout << "c) standard-debug(2) [prints summarized test results]" << std::endl;
      return 0;
   }
   
   if (GoFTStress::fgDebugLevel >= GoFTStress::kBasicDebug) { 
      std::cout << "*****************************************************************************************" << std::endl;
      std::cout << "*                           Goodness of Fit Test STRESS suite                           *" <<std::endl;
      std::cout << "*****************************************************************************************" << std::endl;
}

   gBenchmark = new TBenchmark();
   
   TStopwatch timer;
   timer.Start();
   
   gBenchmark->Start("GoFTestStress");
   
   GoFTStress* goftStressTest = new GoFTStress;
   
   result = goftStressTest->RunTests();
   
   gBenchmark->Stop("GoFTestStress");
   
   goftStressTest->PrintBenchmark();
   
   delete goftStressTest;
   
   delete gBenchmark;
   
   return result;
}

Int_t stressGoFTest(Int_t argc = 1 , Char_t* argv[] = 0) {
   return RunTests(argc, argv);
}
   
#if not defined(__CINT__) && not defined(__MAKECINT__)
Int_t main(Int_t argc, Char_t* argv[]) {
   return RunTests(argc, argv);
}
#endif

---