The class for testing model hypotheses. More...

#include <BCGoFTest.h>

Inheritance diagram for BCGoFTest:

Collaboration diagram for BCGoFTest:

Public Member Functions
Constructors and destructors

	BCGoFTest (const char *name)
	~BCGoFTest ()
Member functions (get)

double	GetCalculatedPValue (bool flag_histogram=false)
TH1D *	GetHistogramLogProb ()
BCModel *	GetTestModel ()
Member functions (set)

void	SetTestModel (BCModel *testmodel)
int	SetTestPoint (std::vector< double > parameters)
Member functions (miscellaneous methods)

double	LogLikelihood (std::vector< double > parameters)
double	LogAPrioriProbability (std::vector< double > parameters)
void	MCMCUserIterationInterface ()
Private Attributes
TH1D *	fHistogramLogProb
double	fLogLikelihood
double	fLogLikelihoodMax
double	fLogLikelihoodMin
std::vector< int >	fMapDataPoint
std::vector< int >	fMapDataValue
int	fPValueAbove
int	fPValueBelow
BCDataSet *	fTemporaryDataSet
BCModel *	fTestModel

Detailed Description

The class for testing model hypotheses.

Author:: Daniel Kollar; Kevin Kröninger

Version:: 1.0

Date:: 08.2008 This class is used for calculating the p-value of a model.

Definition at line 33 of file BCGoFTest.h.

Constructor & Destructor Documentation

BCGoFTest::BCGoFTest ( const char * name )

Default constructor.

Definition at line 20 of file BCGoFTest.cxx.

                                     : BCModel(name)
{
   // set original data set to zero
   fTemporaryDataSet = 0;

   // set test mode to zero
   fTestModel = 0;

   // reset pvalue and counter
   fPValue = 0;
   fPValueAbove = 0;
   fPValueBelow = 0;

   // reset loglikelihood and range
   fLogLikelihood = 0;
   fLogLikelihoodMin = 1e99;
   fLogLikelihoodMax = -1e99;

   // define new histogram
   fHistogramLogProb = 0;

   // set defaults for the MCMC
   this -> MCMCSetNChains(5);
   this -> MCMCSetNIterationsMax(100000);
   this -> MCMCSetNIterationsRun(2000);
}

BCGoFTest::~BCGoFTest ( )

Default destructor.

Definition at line 49 of file BCGoFTest.cxx.

{
   // restore original data set

   // get number of data points and values
   int ndatapoints = fTemporaryDataSet -> GetNDataPoints();
   int ndatavalues = fTemporaryDataSet -> GetDataPoint(0) -> GetNValues();

   for (int i = 0; i < ndatapoints; ++i)
      for (int j = 0; j < ndatavalues; ++j)
         fTestModel -> GetDataSet() -> GetDataPoint(i) -> SetValue(j, fTemporaryDataSet -> GetDataPoint(i) -> GetValue(j));

   // restore data point limits
   for (unsigned int i = 0; i < this -> GetNParameters(); ++i)
      fTestModel -> SetDataBoundaries(
            fMapDataValue[i],
            this -> GetParameter(i) -> GetLowerLimit(),
            this -> GetParameter(i) -> GetUpperLimit());

   // delete temporary data set
   delete fTemporaryDataSet;
}

Member Function Documentation

double BCGoFTest::GetCalculatedPValue ( bool flag_histogram = false )

Definition at line 223 of file BCGoFTest.cxx.

{
   // set histogram point to null
   fHistogramLogProb = 0;

   if (flag_histogram)
   {
      // modify MCMC for first run
//    this -> MCMCSetNIterationsMax(100000);
//    this -> MCMCSetNIterationsRun(10000);

      // perform first run to obtain limits for the log(likelihood)
      this -> MarginalizeAll();

      // modify MCMC for second run
//    this -> MCMCSetNIterationsMax(100000);
//    this -> MCMCSetNIterationsRun(10000);

      // create histogram
      double D = fLogLikelihoodMax - fLogLikelihoodMin;
      fHistogramLogProb = new TH1D(Form("hist_%s_logprob", this -> GetName().data()), ";ln(prob);N", 100, fLogLikelihoodMin - 0.1*D, fLogLikelihoodMax + 0.1*D);
      fHistogramLogProb -> SetStats(kFALSE);
   }
   else
   {
      // modify MCMC
//    this -> MCMCSetNIterationsMax(100000);
//    this -> MCMCSetNIterationsRun(10000);
   }

   // run MCMC
   this -> MarginalizeAll();

   // check for convergence
   if (this -> MCMCGetNIterationsConvergenceGlobal() < 0.)
   {
      BCLog::Out(BCLog::detail, BCLog::detail,
            " --> MCMC did not converge in evaluation of the p-value.");
      return -1;
   }

   // calculate p-value
   fPValue = double(fPValueBelow) / double(fPValueBelow + fPValueAbove);

   // return p-value
   return fPValue;
}

TH1D* BCGoFTest::GetHistogramLogProb ( ) [inline]

Definition at line 61 of file BCGoFTest.h.

         { return fHistogramLogProb; };

BCModel* BCGoFTest::GetTestModel ( ) [inline]

Definition at line 66 of file BCGoFTest.h.

         { return fTestModel; };

double BCGoFTest::LogAPrioriProbability ( std::vector< double > parameters ) [inline, virtual]

Returns natural logarithm of the prior probability. Method needs to be overloaded by the user.

Parameters:

parameters

A set of parameter values

Returns:: The prior probability p(parameters)

See also:: GetPrior(std::vector <double> parameters)

Reimplemented from BCModel.

Definition at line 91 of file BCGoFTest.h.

         { return 0; };

double BCGoFTest::LogLikelihood ( std::vector< double > parameter ) [virtual]

Calculates natural logarithm of the likelihood. Method needs to be overloaded by the user.

Parameters:

parameters

A set of parameter values

Returns:: Natural logarithm of the likelihood

Reimplemented from BCModel.

Definition at line 74 of file BCGoFTest.cxx.

{
   // set the original data set to the new parameters
   for (int i = 0; i < int(parameters.size()); ++i)
      fTestModel -> GetDataSet() -> GetDataPoint(fMapDataPoint[i]) -> SetValue(fMapDataValue[i], parameters.at(i));

   // calculate likelihood at the point of the original parameters
   double loglikelihood = fTestModel -> LogLikelihood(fDataSet -> GetDataPoint(0) -> GetValues());

   // return likelihood
   return loglikelihood;
}

void BCGoFTest::MCMCUserIterationInterface ( ) [virtual]

Reimplemented from BCIntegrate.

Definition at line 89 of file BCGoFTest.cxx.

{
   int nchains = this -> MCMCGetNChains();

   for (int i = 0; i < nchains; ++i)
   {
      // get likelihood at the point of the original parameters
      double loglikelihood = this -> MCMCGetLogProbx(i);

      // calculate pvalue
      if (loglikelihood < fLogLikelihood)
         fPValueBelow++;
      else
         fPValueAbove++;

      // if histogram exists already, then fill it ...
      if (fHistogramLogProb)
         fHistogramLogProb -> Fill(loglikelihood);
      // ...otherwise find range
      else
      {
         if (loglikelihood > fLogLikelihoodMax)
            fLogLikelihoodMax = loglikelihood;
         else if (loglikelihood < fLogLikelihoodMin)
            fLogLikelihoodMin = loglikelihood;
      }
   }
}

void BCGoFTest::SetTestModel ( BCModel * testmodel ) [inline]

Definition at line 76 of file BCGoFTest.h.

         { fTestModel = testmodel; };

int BCGoFTest::SetTestPoint ( std::vector< double > parameters )

Definition at line 120 of file BCGoFTest.cxx.

{
   // check if the boundaries of the original data set exist.
   if (!fTestModel -> GetFlagBoundaries())
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCGoFTest::SetTestDataPoint(). Boundaries of the original data set are not defined.");
      return 0;
   }

   // reset histogram
   if (fHistogramLogProb)
   {
      delete fHistogramLogProb;
      fHistogramLogProb = 0;
   }

   // reset variables
   fPValue = 0;
   fPValueAbove = 0;
   fPValueBelow = 0;

   // create temporary data set ...
   fTemporaryDataSet = new BCDataSet();

   // ... and fill with the original one

   // get number of data points and values
   int ndatapoints = fTestModel -> GetDataSet() -> GetNDataPoints();
   int ndatavalues = fTestModel -> GetDataSet() -> GetDataPoint(0) -> GetNValues();

   for (int i = 0; i < ndatapoints; ++i)
   {
      BCDataPoint * dp = new BCDataPoint(fTestModel -> GetDataSet() -> GetDataPoint(i) -> GetValues());
      fTemporaryDataSet -> AddDataPoint(dp);
   }

   // clear maps
   fMapDataPoint.clear();
   fMapDataValue.clear();

   int counter = 0;

   // remove parameters
   fParameterSet -> clear();
   delete fParameterSet;
   fParameterSet = new BCParameterSet;

   // loop through data points and values
   for (int i = 0; i < ndatapoints; ++i)
      for (int j = 0; j < ndatavalues; ++j)
      {
         if (fTestModel -> GetFixedDataAxis(j))
            continue;

         // add parameter to this model
         this -> AddParameter(
               Form("parameter_%i", counter),
               fTestModel -> GetDataPointLowerBoundary(j),
               fTestModel -> GetDataPointUpperBoundary(j));

         // add another element to the maps
         fMapDataPoint.push_back(i);
         fMapDataValue.push_back(j);

         // increase counter
         counter ++;
      }

   // check if there are any non-fixed data values left
   if (counter == 0)
   {
      BCLog::Out(BCLog::warning, BCLog::warning,"BCGoFTest::SetTestDataPoint(). No non-fixed data values left.");
      return 0;
   }

   // create a new data set containing the vector of parameters which
   // are to be tested
   BCDataPoint * datapoint = new BCDataPoint(parameters);
   BCDataSet * dataset = new BCDataSet();
   dataset -> AddDataPoint(datapoint);

   // calculate likelihood of the original data set
   fLogLikelihood = fTestModel -> LogLikelihood(parameters);

   // if data set has been set before, delete
   if (fDataSet)
      delete fDataSet;

   // set data set of this model
   fDataSet = dataset;

   // put proper range to new data set
   for (int i = 0; i < int(parameters.size()); ++i)
      this -> SetDataBoundaries(
            i,
            fTestModel -> GetParameter(i) -> GetLowerLimit(),
            fTestModel -> GetParameter(i) -> GetUpperLimit());

   return 1;
}