diff --git a/examples/Makefile.am b/examples/Makefile.am index f7fac16..dcb5a7c 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -17,6 +17,7 @@ noinst_PROGRAMS = \ exMathInterpreter \ exMin \ exPlot \ + exPValue \ exRand \ exRootFinder @@ -60,6 +61,10 @@ exPlot_SOURCES = exPlot.cpp exPlot_CXXFLAGS = $(COM_CXXFLAGS) exPlot_LDFLAGS = -L../lib/.libs -lLatAnalyze +exPValue_SOURCES = exPValue.cpp +exPValue_CXXFLAGS = $(COM_CXXFLAGS) +exPValue_LDFLAGS = -L../lib/.libs -lLatAnalyze + exRand_SOURCES = exRand.cpp exRand_CXXFLAGS = $(COM_CXXFLAGS) exRand_LDFLAGS = -L../lib/.libs -lLatAnalyze diff --git a/examples/exPValue.cpp b/examples/exPValue.cpp new file mode 100644 index 0000000..7f6a745 --- /dev/null +++ b/examples/exPValue.cpp @@ -0,0 +1,23 @@ +#include + +using namespace std; +using namespace Latan; + +int main(int argc, char* argv[]) +{ + double chi2, ndof; + + if (argc != 3) + { + cerr << "usage: " << argv[0] << " " << endl; + + return EXIT_FAILURE; + } + chi2 = strTo(argv[1]); + ndof = strTo(argv[2]); + + cout << "Two-sided p-value: " << Math::chi2PValue(chi2, ndof) << endl; + cout << "chi^2 CCDF : " << Math::chi2Ccdf(chi2, ndof) << endl; + + return EXIT_SUCCESS; +} diff --git a/lib/Core/Math.cpp b/lib/Core/Math.cpp index cca44dd..25c8495 100644 --- a/lib/Core/Math.cpp +++ b/lib/Core/Math.cpp @@ -107,8 +107,14 @@ DEF_STD_FUNC_1ARG(fabs) // p-value auto chi2PValueVecFunc = [](const double arg[2]) +{ + return 2.*min(gsl_cdf_chisq_P(arg[0], arg[1]), gsl_cdf_chisq_Q(arg[0], arg[1])); +}; + +auto chi2CcdfVecFunc = [](const double arg[2]) { return gsl_cdf_chisq_Q(arg[0], arg[1]); }; DoubleFunction MATH_NAMESPACE::chi2PValue(chi2PValueVecFunc, 2); +DoubleFunction MATH_NAMESPACE::chi2Ccdf(chi2CcdfVecFunc, 2); diff --git a/lib/Core/Math.hpp b/lib/Core/Math.hpp index cac584a..9ea5008 100644 --- a/lib/Core/Math.hpp +++ b/lib/Core/Math.hpp @@ -152,6 +152,7 @@ DECL_STD_FUNC(fabs) namespace MATH_NAMESPACE { extern DoubleFunction chi2PValue; + extern DoubleFunction chi2Ccdf; } END_LATAN_NAMESPACE diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 8f04a7d..4c46f70 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -62,6 +62,11 @@ double SampleFitResult::getPValue(const Index s) const return Math::chi2PValue(getChi2(s), getNDof()); } +double SampleFitResult::getCcdf(const Index s) const +{ + return Math::chi2Ccdf(getChi2(s), getNDof()); +} + const DoubleFunction & SampleFitResult::getModel(const Index s, const Index j) const { @@ -98,8 +103,9 @@ void SampleFitResult::print(const bool printXsi, ostream &out) const Index pMax = printXsi ? size() : nPar_; DMat err = this->variance().cwiseSqrt(); - sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- p-value= %.2e", getChi2(), - static_cast(getNDof()), getChi2PerDof(), getPValue()); + sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- chi^2 CCDF= %.2e -- p-value= %.2e", + getChi2(), static_cast(getNDof()), getChi2PerDof(), getCcdf(), + getPValue()); out << buf << endl; for (Index p = 0; p < pMax; ++p) { diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index 5bf8cf2..d7a4c7b 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -49,6 +49,7 @@ public: double getNDof(void) const; Index getNPar(void) const; double getPValue(const Index s = central) const; + double getCcdf(const Index s = central) const; const DoubleFunction & getModel(const Index s = central, const Index j = 0) const; const DoubleFunctionSample & getModel(const PlaceHolder ph, diff --git a/lib/Statistics/XYStatData.cpp b/lib/Statistics/XYStatData.cpp index b3acca6..32e4595 100644 --- a/lib/Statistics/XYStatData.cpp +++ b/lib/Statistics/XYStatData.cpp @@ -55,6 +55,11 @@ double FitResult::getPValue(void) const return Math::chi2PValue(getChi2(), getNDof());; } +double FitResult::getCcdf(void) const +{ + return Math::chi2Ccdf(getChi2(), getNDof());; +} + const DoubleFunction & FitResult::getModel(const Index j) const { return model_[j]; @@ -66,8 +71,9 @@ void FitResult::print(const bool printXsi, ostream &out) const char buf[256]; Index pMax = printXsi ? size() : nPar_; - sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- p-value= %.2e", getChi2(), - static_cast(getNDof()), getChi2PerDof(), getPValue()); + sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- chi^2 CCDF= %.2e -- p-value= %.2e", + getChi2(), static_cast(getNDof()), getChi2PerDof(), getCcdf(), + getPValue()); out << buf << endl; for (Index p = 0; p < pMax; ++p) { diff --git a/lib/Statistics/XYStatData.hpp b/lib/Statistics/XYStatData.hpp index 1bdd277..e8275dc 100644 --- a/lib/Statistics/XYStatData.hpp +++ b/lib/Statistics/XYStatData.hpp @@ -47,6 +47,7 @@ public: double getNDof(void) const; Index getNPar(void) const; double getPValue(void) const; + double getCcdf(void) const; const DoubleFunction & getModel(const Index j = 0) const; // IO void print(const bool printXsi = false,