The class for testing model hypotheses. More...

#include <BCGoFTest.h>

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 (const std::vector< double > &parameters)
double	LogAPrioriProbability (const std::vector< double > &)
void	MCMCUserIterationInterface ()
Private Attributes
std::vector< int >	fMapDataPoint
std::vector< int >	fMapDataValue
int	fPValueBelow
int	fPValueAbove
BCModel *	fTestModel
BCDataSet *	fTemporaryDataSet
double	fLogLikelihood
double	fLogLikelihoodMin
double	fLogLikelihoodMax
TH1D *	fHistogramLogProb

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
   MCMCSetNChains(5);
   MCMCSetNIterationsMax(100000);
   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 < GetNParameters(); ++i)
      fTestModel->SetDataBoundaries(
            fMapDataValue[i],
            GetParameter(i)->GetLowerLimit(),
            GetParameter(i)->GetUpperLimit());

   // delete temporary data set
   delete fTemporaryDataSet;
}

Member Function Documentation

double BCGoFTest::GetCalculatedPValue ( bool flag_histogram = false )

Calculated the p-value.

Parameters:

flag_histogram

A histogram is either filled or not.

Returns:: p-value

Definition at line 223 of file BCGoFTest.cxx.

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

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

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

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

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

   // run MCMC
   MarginalizeAll();

   // check for convergence
   if (MCMCGetNIterationsConvergenceGlobal() < 0.)
   {
      BCLog::OutDetail(" --> 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]

Returns:: distribution of log(likelihood)

Definition at line 61 of file BCGoFTest.h.

         { return fHistogramLogProb; };

BCModel* BCGoFTest::GetTestModel ( ) [inline]

Returns:: pointer to the tested model

Definition at line 66 of file BCGoFTest.h.

         { return fTestModel; };

double BCGoFTest::LogAPrioriProbability ( const 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 ( const std::vector< double > & parameters ) [virtual]

Implements 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]

Method executed for every iteration of the MCMC. User's code should be provided via overloading in the derived class

Reimplemented from BCIntegrate.

Definition at line 89 of file BCGoFTest.cxx.

{
   int nchains = MCMCGetNChains();

   for (int i = 0; i < nchains; ++i)
   {
      // get likelihood at the point of the original parameters
      double loglikelihood = 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]

Set the model to be tested.

Parameters:

testmodel

pointer to the model to be tested

Definition at line 76 of file BCGoFTest.h.

         { fTestModel = testmodel; };

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

Sets the set of parameters which the p-values is calculated for.

Parameters:

parameters

parameters

Returns:: error code

Definition at line 120 of file BCGoFTest.cxx.

{
   // check if the boundaries of the original data set exist.
   if (!fTestModel->GetFlagBoundaries())
   {
      BCLog::OutError("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
         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::OutError("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)
      SetDataBoundaries(
            i,
            fTestModel->GetParameter(i)->GetLowerLimit(),
            fTestModel->GetParameter(i)->GetUpperLimit());

   return 1;
}