Skip to content
Snippets Groups Projects
Commit 602c0490 authored by Wouter Verkerke's avatar Wouter Verkerke
Browse files

  o RooPoisson

    - Add flag to disable rounding (request M.Baak)


  o RooJeffreysPrior, RooNonCentralChiSquare

    - Two new pdfs from Kyle to be used in RooStats


git-svn-id: http://root.cern.ch/svn/root/trunk@37207 27541ba8-7e3a-0410-8455-c3a389f83636
parent 329b455b
Branches
Tags
No related merge requests found
...@@ -50,6 +50,8 @@ ...@@ -50,6 +50,8 @@
#pragma link C++ class RooUniform+ ; #pragma link C++ class RooUniform+ ;
#pragma link C++ class RooSpHarmonic+ ; #pragma link C++ class RooSpHarmonic+ ;
#pragma link C++ class RooLegendre+ ; #pragma link C++ class RooLegendre+ ;
#pragma link C++ class RooNonCentralChiSquare+ ;
#pragma link C++ class RooJeffreysPrior+ ;
#pragma link C++ class RooFunctorBinding+ ; #pragma link C++ class RooFunctorBinding+ ;
#pragma link C++ class RooFunctor1DBinding+ ; #pragma link C++ class RooFunctor1DBinding+ ;
#pragma link C++ class RooFunctorPdfBinding+ ; #pragma link C++ class RooFunctorPdfBinding+ ;
......
/*****************************************************************************
* Project: RooStats
* Package: RooStats
* File: $Id$
* author: Kyle Cranmer
*****************************************************************************/
#ifndef JEFFREYSPRIOR
#define JEFFREYSPRIOR
#include "RooAbsPdf.h"
#include "RooRealProxy.h"
#include "RooListProxy.h"
class RooRealVar;
class RooArgList ;
class RooJeffreysPrior : public RooAbsPdf {
public:
RooJeffreysPrior() ;
RooJeffreysPrior(const char *name, const char *title, RooAbsPdf& nominal, const RooArgList& paramSet, const RooArgList& obsSet) ;
virtual ~RooJeffreysPrior() ;
RooJeffreysPrior(const RooJeffreysPrior& other, const char* name = 0);
virtual TObject* clone(const char* newname) const { return new RooJeffreysPrior(*this, newname); }
const RooArgList& lowList() const { return _obsSet ; }
const RooArgList& paramList() const { return _paramSet ; }
Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ;
Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ;
protected:
RooRealProxy _nominal; // The nominal value
//RooAbsPdf* _nominal; // The nominal value
RooArgList _ownedList ; // List of owned components
RooListProxy _obsSet ; // Low-side variation
RooListProxy _paramSet ; // interpolation parameters
mutable TIterator* _paramIter ; //! Iterator over paramSet
mutable TIterator* _obsIter ; //! Iterator over lowSet
Double_t evaluate() const;
ClassDef(RooJeffreysPrior,1) // Sum of RooAbsReal objects
};
#endif
/*****************************************************************************
* Project: RooFit *
* @(#)root/roofit:$Id$ *
* *
* RooFit NonCentralChisquare PDF *
* *
* Author: Kyle Cranmer *
* *
*****************************************************************************/
#ifndef ROO_NONCENTRALCHISQUARE
#define ROO_NONCENTRALCHISQUARE
#include "RooAbsPdf.h"
#include "RooRealProxy.h"
#include "RooCategoryProxy.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
class RooNonCentralChiSquare : public RooAbsPdf {
public:
RooNonCentralChiSquare() {} ;
RooNonCentralChiSquare(const char *name, const char *title,
RooAbsReal& _x,
RooAbsReal& _k,
RooAbsReal& _lambda);
RooNonCentralChiSquare(const RooNonCentralChiSquare& other, const char* name=0) ;
virtual TObject* clone(const char* newname) const { return new RooNonCentralChiSquare(*this,newname); }
inline virtual ~RooNonCentralChiSquare() { }
void SetErrorTolerance(Double_t t) {fErrorTol = t;}
void SetMaxIters(Int_t mi) {fMaxIters = mi;}
void SetForceSum(Bool_t flag);
Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ;
Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ;
protected:
RooRealProxy x ;
RooRealProxy k ;
RooRealProxy lambda ;
Double_t fErrorTol;
Int_t fMaxIters;
Bool_t fForceSum;
mutable Bool_t fHasIssuedConvWarning;
mutable Bool_t fHasIssuedSumWarning;
Double_t evaluate() const ;
private:
ClassDef(RooNonCentralChiSquare,1) // non-central chisquare pdf
};
#endif
...@@ -17,8 +17,8 @@ ...@@ -17,8 +17,8 @@
class RooPoisson : public RooAbsPdf { class RooPoisson : public RooAbsPdf {
public: public:
RooPoisson() {} ; RooPoisson() { _noRounding = kFALSE ; } ;
RooPoisson(const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _mean); RooPoisson(const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _mean, Bool_t noRounding=kFALSE);
RooPoisson(const RooPoisson& other, const char* name=0) ; RooPoisson(const RooPoisson& other, const char* name=0) ;
virtual TObject* clone(const char* newname) const { return new RooPoisson(*this,newname); } virtual TObject* clone(const char* newname) const { return new RooPoisson(*this,newname); }
inline virtual ~RooPoisson() { } inline virtual ~RooPoisson() { }
...@@ -33,14 +33,15 @@ protected: ...@@ -33,14 +33,15 @@ protected:
RooRealProxy x ; RooRealProxy x ;
RooRealProxy mean ; RooRealProxy mean ;
Bool_t _noRounding ;
Double_t evaluate() const ; Double_t evaluate() const ;
Double_t evaluate(Double_t k) const; Double_t evaluate(Double_t k) const;
private: private:
ClassDef(RooPoisson,1) // A Poisson PDF ClassDef(RooPoisson,2) // A Poisson PDF
}; };
#endif #endif
/*****************************************************************************
*****************************************************************************/
//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// RooJeffreysPrior
// END_HTML
//
#include "RooFit.h"
#include "Riostream.h"
#include "Riostream.h"
#include <math.h>
#include "RooJeffreysPrior.h"
#include "RooAbsReal.h"
#include "RooAbsPdf.h"
#include "RooErrorHandler.h"
#include "RooArgSet.h"
#include "RooMsgService.h"
#include "RooFitResult.h"
#include "TMatrixDSym.h"
#include "RooDataHist.h"
#include "RooFitResult.h"
#include "RooNumIntConfig.h"
#include "RooRealVar.h"
ClassImp(RooJeffreysPrior)
;
using namespace RooFit;
//_____________________________________________________________________________
RooJeffreysPrior::RooJeffreysPrior(const char* name, const char* title,
RooAbsPdf& nominal,
const RooArgList& paramSet,
const RooArgList& obsSet) :
RooAbsPdf(name, title),
_nominal("nominal","nominal",this,nominal,kFALSE,kFALSE),
_obsSet("!obsSet","obs-side variation",this,kFALSE,kFALSE),
_paramSet("!paramSet","high-side variation",this)
{
//_obsSet("!obsSet","obs-side variation",this),
_obsIter = _obsSet.createIterator() ;
_paramIter = _paramSet.createIterator() ;
TIterator* inputIter1 = obsSet.createIterator() ;
RooAbsArg* comp ;
while((comp = (RooAbsArg*)inputIter1->Next())) {
if (!dynamic_cast<RooAbsReal*>(comp)) {
coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
<< " in first list is not of type RooAbsReal" << endl ;
RooErrorHandler::softAbort() ;
}
_obsSet.add(*comp) ;
// if (takeOwnership) {
// _ownedList.addOwned(*comp) ;
// }
}
delete inputIter1 ;
TIterator* inputIter3 = paramSet.createIterator() ;
while((comp = (RooAbsArg*)inputIter3->Next())) {
if (!dynamic_cast<RooAbsReal*>(comp)) {
coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
<< " in first list is not of type RooAbsReal" << endl ;
RooErrorHandler::softAbort() ;
}
_paramSet.add(*comp) ;
// if (takeOwnership) {
// _ownedList.addOwned(*comp) ;
// }
}
delete inputIter3 ;
// use a different integrator by default.
if(paramSet.getSize()==1)
this->specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
}
//_____________________________________________________________________________
RooJeffreysPrior::RooJeffreysPrior(const RooJeffreysPrior& other, const char* name) :
RooAbsPdf(other, name),
_nominal("!nominal",this,other._nominal),
_obsSet("!obsSet",this,other._obsSet),
_paramSet("!paramSet",this,other._paramSet)
{
// Copy constructor
_obsIter = _obsSet.createIterator() ;
_paramIter = _paramSet.createIterator() ;
// Member _ownedList is intentionally not copy-constructed -- ownership is not transferred
}
//_____________________________________________________________________________
RooJeffreysPrior::RooJeffreysPrior()
{
// Default constructor
_obsIter = NULL;
_paramIter = NULL;
}
//_____________________________________________________________________________
RooJeffreysPrior::~RooJeffreysPrior()
{
// Destructor
if (_obsIter) delete _obsIter ;
if (_paramIter) delete _paramIter ;
}
//_____________________________________________________________________________
Double_t RooJeffreysPrior::evaluate() const
{
// Calculate and return current value of self
RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
// create Asimov dataset
// _paramSet.Print("v");
// return sqrt(_paramSet.getRealValue("mu"));
*(_nominal.arg().getVariables()) = _paramSet;
/*
cout << "_______________"<<endl;
_paramSet.Print("v");
_nominal->getVariables()->Print("v");
cout << "_______________"<<endl;
*/
RooDataHist* data = ((RooAbsPdf&)(_nominal.arg())).generateBinned(_obsSet,ExpectedData());
// data->Print("v");
//RooFitResult* res = _nominal->fitTo(*data, Save(),PrintLevel(-1),Minos(kFALSE),SumW2Error(kTRUE));
RooFitResult* res = ((RooAbsPdf&)(_nominal.arg())).fitTo(*data, Save(),PrintLevel(-1),Minos(kFALSE),SumW2Error(kFALSE));
TMatrixDSym cov = res->covarianceMatrix();
cov.Invert();
double ret = sqrt(cov.Determinant());
/*
// for 1 parameter can avoid making TMatrix etc.
// but number of params may be > 1 with others held constant
if(_paramSet.getSize()==1){
RooRealVar* var = (RooRealVar*) _paramSet.first();
// also, the _paramSet proxy one does not pick up a different value
cout << "eval at "<< ret << " " << 1/(var->getError()) << endl;
// need to get the actual variable instance out of the pdf like below
var = (RooRealVar*) _nominal->getVariables()->find(var->GetName());
cout << "eval at "<< ret << " " << 1/(var->getError()) << endl;
}
*/
// res->Print();
delete data;
delete res;
RooMsgService::instance().setGlobalKillBelow(msglevel);
// cout << "eval at "<< ret << endl;
// _paramSet.Print("v");
return ret;
}
//_____________________________________________________________________________
Int_t RooJeffreysPrior::getAnalyticalIntegral(RooArgSet& /*allVars*/, RooArgSet& /*analVars*/, const char* /*rangeName*/) const
{
// if (matchArgs(allVars,analVars,x)) return 1 ;
// if (matchArgs(allVars,analVars,mean)) return 2 ;
// return 1;
return 0 ;
}
//_____________________________________________________________________________
Double_t RooJeffreysPrior::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
{
assert(code==1 );
//cout << "evaluating analytic integral" << endl;
return 1.;
}
/*****************************************************************************
* Project: RooFit *
* Package: RooFitModels *
* @(#)root/roofit:$Id: RooNonCentralChiSquare *
* Authors: *
* Kyle Cranmer
* *
*****************************************************************************/
//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// The PDF of the Non-Central Chi Square distribution for n degrees of freedom.
// It is the asymptotic distribution of the profile likeihood ratio test q_mu
// when a different mu' is true. It is Wald's generalization of Wilks' Theorem.
//
// See:
// Asymptotic formulae for likelihood-based tests of new physics
// By Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells
// http://arXiv.org/abs/arXiv:1007.1727
//
// Wikipedia:
// http://en.wikipedia.org/wiki/Noncentral_chi-square_distribution#Approximation
//
// It requires MathMore to evaluate for non-integer degrees of freedom, k.
//
// When the Mathmore library is available we can use the MathMore libraries impelmented using GSL.
// It makes use of the modified Bessel function of the first kind (for k > 2). For k < 2 it uses
// the hypergeometric function 0F1.
// When is not available we use explicit summation of normal chi-squared distributions
// The usage of the sum can be forced by calling SetForceSum(true);
//
// This implementation could be improved. BOOST has a nice implementation:
// http://live.boost.org/doc/libs/1_42_0/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/nc_chi_squared_dist.html
// http://wesnoth.repositoryhosting.com/trac/wesnoth_wesnoth/browser/trunk/include/boost/math/distributions/non_central_chi_squared.hpp?rev=6
// END_HTML
//
#include "Riostream.h"
#include "RooNonCentralChiSquare.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include <math.h>
#include "TMath.h"
//#include "RooNumber.h"
#include "Math/DistFunc.h"
#include "RooMsgService.h"
ClassImp(RooNonCentralChiSquare)
RooNonCentralChiSquare::RooNonCentralChiSquare(const char *name, const char *title,
RooAbsReal& _x,
RooAbsReal& _k,
RooAbsReal& _lambda) :
RooAbsPdf(name,title),
x("x","x",this,_x),
k("k","k",this,_k),
lambda("lambda","lambda",this,_lambda),
fErrorTol(1E-3),
fMaxIters(10),
fHasIssuedConvWarning(false),
fHasIssuedSumWarning(false)
{
#ifdef R__HAS_MATHMORE
ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
"MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
fForceSum = false;
#else
fForceSum = true;
#endif
}
RooNonCentralChiSquare::RooNonCentralChiSquare(const RooNonCentralChiSquare& other, const char* name) :
RooAbsPdf(other,name),
x("x",this,other.x),
k("k",this,other.k),
lambda("lambda",this,other.lambda),
fErrorTol(other.fErrorTol),
fMaxIters(other.fMaxIters),
fHasIssuedConvWarning(false),
fHasIssuedSumWarning(false)
{
#ifdef R__HAS_MATHMORE
ccoutD(InputArguments) << "RooNonCentralChiSquare::ctor(" << GetName() <<
"MathMore Available, will use Bessel function expressions unless SetForceSum(true) "<< endl ;
fForceSum = other.fForceSum;
#else
fForceSum = true;
#endif
}
void RooNonCentralChiSquare::SetForceSum(Bool_t flag) {
fForceSum = flag;
#ifndef R__HAS_MATHMORE
if (!fForceSum) {
ccoutD(InputArguments) << "RooNonCentralChiSquare::SetForceSum" << GetName() <<
"MathMore is not available- ForceSum must be on "<< endl ;
fForceSum = true;
}
#endif
}
Double_t RooNonCentralChiSquare::evaluate() const
{
// ENTER EXPRESSION IN TERMS OF VARIABLE ARGUMENTS HERE
// chi^2(0,k) gives inf and causes various problems
// truncate
Double_t xmin = x.min();
Double_t xmax = x.max();
double _x = x;
if(_x<=0){
// options for dealing with this
// return 0; // gives a funny dip
// _x = 1./RooNumber::infinity(); // too tall
_x = xmin + 1e-3*(xmax-xmin); // very small fraction of range
}
// special case (form below doesn't work when lambda==0)
if(lambda==0){
return ROOT::Math::chisquared_pdf(_x,k);
}
// three forms
// FIRST FORM
// \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
// could truncate sum
if ( fForceSum ){
if(!fHasIssuedSumWarning){
coutI(InputArguments) << "RooNonCentralChiSquare sum being forced" << endl ;
fHasIssuedSumWarning=true;
}
double sum = 0;
double ithTerm = 0;
double errorTol = fErrorTol;
int MaxIters = fMaxIters;
int iDominant = (int) TMath::Floor(lambda/2);
// cout <<"iDominant: " << iDominant << endl;
// do 0th term last
// if(iDominant==0) iDominant = 1;
for(int i = iDominant; ; ++i){
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
sum+=ithTerm;
// cout <<"progress: " << i << " " << ithTerm/sum << endl;
if(ithTerm/sum < errorTol)
break;
if( i>iDominant+MaxIters){
if(!fHasIssuedConvWarning){
fHasIssuedConvWarning=true;
coutW(Eval) << "RooNonCentralChiSquare did not converge: for x=" << x <<" k="<<k
<< ", lambda="<<lambda << " fractional error = " << ithTerm/sum
<< "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
<< endl;
}
break;
}
}
for(int i = iDominant - 1; i >= 0; --i){
// cout <<"Progress: " << i << " " << ithTerm/sum << endl;
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)*ROOT::Math::chisquared_pdf(_x,k+2*i)/TMath::Gamma(i+1);
sum+=ithTerm;
}
return sum;
}
// SECOND FORM (use MathMore function based on Bessel function (if k>2) or
// or regularized confluent hypergeometric limit function.
#ifdef R__HAS_MATHMORE__DISABLEFORNOW
return ROOT::Math::noncentral_chisquared_pdf(_x,k,lambda);
#else
coutF(Eval) << "RooNonCentralChisquare: ForceSum must be set" << endl;
return 0;
#endif
}
//_____________________________________________________________________________
Int_t RooNonCentralChiSquare::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
{
if (matchArgs(allVars,analVars,x)) return 1 ;
return 0 ;
}
//_____________________________________________________________________________
Double_t RooNonCentralChiSquare::analyticalIntegral(Int_t code, const char* rangeName) const
{
assert(code==1 );
// cout << "evaluating analytic integral" << endl;
Double_t xmin = x.min(rangeName);
Double_t xmax = x.max(rangeName);
// if xmin~0 and xmax big, then can return 1. b/c evaluate is normalized.
// special case (form below doesn't work when lambda==0)
if(lambda==0){
return (ROOT::Math::chisquared_cdf(xmax,k) - ROOT::Math::chisquared_cdf(xmin,k));
}
// three forms
// FIRST FORM
// \sum_i=0^\infty exp(-lambda/2) (\lamda/2)^i chi2(x,k+2i) / i!
// could truncate sum
double sum = 0;
double ithTerm = 0;
double errorTol = fErrorTol; // for nomralization allow slightly larger error
int MaxIters = fMaxIters; // for normalization use more terms
int iDominant = (int) TMath::Floor(lambda/2);
// cout <<"iDominant: " << iDominant << endl;
// iDominant=0;
for(int i = iDominant; ; ++i){
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
*( ROOT::Math::chisquared_cdf(xmax,k+2*i)/TMath::Gamma(i+1)
- ROOT::Math::chisquared_cdf(xmin,k+2*i)/TMath::Gamma(i+1) );
sum+=ithTerm;
// cout <<"progress: " << i << " " << ithTerm << " " << sum << endl;
if(ithTerm/sum < errorTol)
break;
if( i>iDominant+MaxIters){
if(!fHasIssuedConvWarning){
fHasIssuedConvWarning=true;
coutW(Eval) << "RooNonCentralChiSquare Normalization did not converge: for k="<<k
<< ", lambda="<<lambda << " fractional error = " << ithTerm/sum
<< "\n either adjust tolerance with SetErrorTolerance(tol) or max_iter with SetMaxIter(max_it)"
<< endl;
}
break;
}
}
for(int i = iDominant - 1; i >= 0; --i){
ithTerm =exp(-lambda/2.)*pow(lambda/2.,i)
*( ROOT::Math::chisquared_cdf(xmax,k+2*i)/TMath::Gamma(i+1)
-ROOT::Math::chisquared_cdf(xmin,k+2*i)/TMath::Gamma(i+1) );
sum+=ithTerm;
}
return sum;
}
...@@ -31,10 +31,12 @@ ClassImp(RooPoisson) ...@@ -31,10 +31,12 @@ ClassImp(RooPoisson)
//_____________________________________________________________________________ //_____________________________________________________________________________
RooPoisson::RooPoisson(const char *name, const char *title, RooPoisson::RooPoisson(const char *name, const char *title,
RooAbsReal& _x, RooAbsReal& _x,
RooAbsReal& _mean) : RooAbsReal& _mean,
Bool_t noRounding) :
RooAbsPdf(name,title), RooAbsPdf(name,title),
x("x","x",this,_x), x("x","x",this,_x),
mean("mean","mean",this,_mean) mean("mean","mean",this,_mean),
_noRounding(noRounding)
{ {
// Constructor // Constructor
} }
...@@ -45,7 +47,8 @@ RooPoisson::RooPoisson(const char *name, const char *title, ...@@ -45,7 +47,8 @@ RooPoisson::RooPoisson(const char *name, const char *title,
RooPoisson::RooPoisson(const RooPoisson& other, const char* name) : RooPoisson::RooPoisson(const RooPoisson& other, const char* name) :
RooAbsPdf(other,name), RooAbsPdf(other,name),
x("x",this,other.x), x("x",this,other.x),
mean("mean",this,other.mean) mean("mean",this,other.mean),
_noRounding(other._noRounding)
{ {
// Copy constructor // Copy constructor
} }
...@@ -58,7 +61,7 @@ Double_t RooPoisson::evaluate() const ...@@ -58,7 +61,7 @@ Double_t RooPoisson::evaluate() const
{ {
// Implementation in terms of the TMath Poisson function // Implementation in terms of the TMath Poisson function
Double_t k = floor(x); Double_t k = _noRounding ? x : floor(x);
return TMath::Poisson(k,mean) ; return TMath::Poisson(k,mean) ;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment