p-value now 2-sided and added chi^2 CCDF

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

examples/exPValue.cpp

@@ -0,0 +1,23 @@
+#include <LatAnalyze/Core/Math.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char* argv[])
+{
+    double chi2, ndof;
+
+    if (argc != 3)
+    {
+        cerr << "usage: " << argv[0] << " <chi^2> <ndof>" << endl;
+        
+        return EXIT_FAILURE;
+    }
+    chi2 = strTo<double>(argv[1]);
+    ndof = strTo<double>(argv[2]);
+
+    cout << "Two-sided p-value: " << Math::chi2PValue(chi2, ndof) << endl;
+    cout << "chi^2 CCDF       : " << Math::chi2Ccdf(chi2, ndof) << endl;
+
+    return EXIT_SUCCESS;
+}

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);

lib/Core/Math.hpp

@@ -152,6 +152,7 @@ DECL_STD_FUNC(fabs)
 namespace MATH_NAMESPACE
 {
     extern DoubleFunction chi2PValue;
+    extern DoubleFunction chi2Ccdf;
 }
 
 END_LATAN_NAMESPACE

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(),
+    sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- chi^2 CCDF= %.2e -- p-value= %.2e",
-            static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
+        getChi2(), static_cast<int>(getNDof()), getChi2PerDof(), getCcdf(), 
+        getPValue());
     out << buf << endl;
     for (Index p = 0; p < pMax; ++p)
     {

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,

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(),
+    sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- chi^2 CCDF= %.2e -- p-value= %.2e",
-            static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
+        getChi2(), static_cast<int>(getNDof()), getChi2PerDof(), getCcdf(), 
+        getPValue());
     out << buf << endl;
     for (Index p = 0; p < pMax; ++p)
     {

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,