[apps/proba] Student cumulativeDistributiveInverseForProbability

apps/probability/Makefile

@@ -3,6 +3,7 @@ app_headers += apps/probability/app.h
 
 app_probability_test_src = $(addprefix apps/probability/,\
   law/erf_inv.cpp \
+  law/incomplete_beta_function.cpp \
   law/law.cpp \
   law/regularized_gamma.cpp \
   law/helper.cpp \
@@ -47,6 +48,9 @@ i18n_files += $(addprefix apps/probability/,\
   base.universal.i18n\
 )
 
+tests_src += $(addprefix apps/probability/law/,\
+  hypergeometric_function.cpp\
+)
 tests_src += $(addprefix apps/probability/test/,\
   erf_inv.cpp\
   regularized_gamma.cpp\

apps/probability/law/helper.cpp

@@ -59,7 +59,7 @@ bool Helper::ContinuedFractionEvaluation(ValueForIndex a, ValueForIndex b, int m
   return true;
 }
 
-bool Helper::InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater, double termLimit, int maxNumberOfIterations, double * result, double context1, double context2) {
+bool Helper::InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater, double termLimit, int maxNumberOfIterations, double * result, double context1, double context2, double context3, double context4) {
   double iterationCount = 0.0;
   double currentTerm = firstTerm;
   double sum = currentTerm;
@@ -68,7 +68,7 @@ bool Helper::InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater,
       && std::fabs(currentTerm/sum) > termLimit)
   {
     iterationCount+= 1.0;
-    currentTerm = termUpdater(currentTerm, iterationCount, context1, context2);
+    currentTerm = termUpdater(currentTerm, iterationCount, context1, context2, context3, context4);
     sum+= currentTerm;
   }
   if (iterationCount >= maxNumberOfIterations) {

apps/probability/law/helper.h

@@ -8,8 +8,8 @@ public:
   static bool ContinuedFractionEvaluation(ValueForIndex a, ValueForIndex b, int maxNumberOfIterations, double * result, double context1, double context2);
 
   // Infinite series
-  typedef double (*TermUpdater)(double previousTerm, double index, double s, double x);
-  static bool InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater, double termLimit, int maxNumberOfIterations, double * result, double context1, double context2);
+  typedef double (*TermUpdater)(double previousTerm, double index, double d1, double d2, double d3, double d4);
+  static bool InfiniteSeriesEvaluation(double firstTerm, TermUpdater termUpdater, double termLimit, int maxNumberOfIterations, double * result, double context1, double context2, double context3 = 0.0, double context4 = 0.0);
 
 };
 

apps/probability/law/hypergeometric_function.cpp

@@ -0,0 +1,32 @@
+#include "hypergeometric_function.h"
+#include "helper.h"
+#include <cmath>
+#include <float.h>
+#include <assert.h>
+
+
+bool hypergeometricFunction(double a, double b, double c, double z, double epsilon, int maxNumberOfIterations, double * result) {
+  // TODO Put interruption instead of maxNumberOfIterations
+  assert(!std::isnan(a) && !std::isnan(b) && !std::isnan(c) && !std::isnan(z));
+  if (z == 0.0) {
+    *result = 0.0;
+    return true;
+  }
+  /* For x < s + 1: Compute using the infinite series representation
+   * hypergeometricFunction(a,b,c,z) = sum((a)n * (b)n * z^n / ((c)n * n!...)
+   * With (a)n = a(a+1)..(a+n-1)*/
+  if (std::fabs(z) < 1.0) {
+    assert(c > 0.0);
+    return Helper::InfiniteSeriesEvaluation(
+        1.0,
+        [](double previousTerm, double index, double a, double b, double c, double z) { return previousTerm * (a + index - 1) * (b + index - 1) * z / ((c + index - 1) * index); },
+        epsilon,
+        maxNumberOfIterations,
+        result,
+        a,
+        b,
+        c,
+        z);
+  }
+  return NAN; //TODO
+}

apps/probability/law/hypergeometric_function.h

@@ -0,0 +1,20 @@
+#ifndef PROBABILITY_HYPERGEOMETRIC_FUNCTION_H
+#define PROBABILITY_HYPERGEOMETRIC_FUNCTION_H
+
+/* This code can be used to compute the Student law for |x| < root(k).
+ * We do not use it because we want to cover more x, but we keep in case we need
+ * it later.  */
+
+/* hypergeometricFunction(a, b, c, z) = 2F1(a, b, c, z)
+ *           an * bn    z^n
+ *    = sum( -------- * --- )
+ *             cn        n!
+ *
+ * with :
+ * an = 1 if n = 0
+ * an = (a*(a+1)*...*(a+n-1) otherwise */
+
+bool hypergeometricFunction(double a, double b, double c, double z, double epsilon, int maxNumberOfIterations, double * result);
+
+#endif
+

apps/probability/law/incomplete_beta_function.cpp

@@ -0,0 +1,82 @@
+/*
+ * zlib License
+ *
+ * Regularized Incomplete Beta Function
+ *
+ * Copyright (c) 2016, 2017 Lewis Van Winkle
+ * http://CodePlea.com
+ *
+ * This software is provided 'as-is', without any express or implied
+ * warranty. In no event will the authors be held liable for any damages
+ * arising from the use of this software.
+ *
+ * Permission is granted to anyone to use this software for any purpose,
+ * including commercial applications, and to alter it and redistribute it
+ * freely, subject to the following restrictions:
+ *
+ * 1. The origin of this software must not be misrepresented; you must not
+ *    claim that you wrote the original software. If you use this software
+ *    in a product, an acknowledgement in the product documentation would be
+ *    appreciated but is not required.
+ * 2. Altered source versions must be plainly marked as such, and must not be
+ *    misrepresented as being the original software.
+ * 3. This notice may not be removed or altered from any source distribution.
+ */
+
+// WARNING: this code has been modified
+
+#include "incomplete_beta_function.h"
+#include <math.h>
+#include <cmath>
+
+#define STOP 1.0e-8
+#define TINY 1.0e-30
+
+double IncompleteBetaFunction(double a, double b, double x) {
+    if (x < 0.0 || x > 1.0) return NAN;
+
+    /*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/
+    if (x > (a+1.0)/(a+b+2.0)) {
+        return (1.0-IncompleteBetaFunction(b,a,1.0-x)); /*Use the fact that beta is symmetrical.*/
+    }
+
+    /*Find the first part before the continued fraction.*/
+    const double lbeta_ab = std::lgamma(a)+std::lgamma(b)-std::lgamma(a+b);
+    const double front = std::exp(std::log(x)*a+std::log(1.0-x)*b-lbeta_ab) / a;
+
+    /*Use Lentz's algorithm to evaluate the continued fraction.*/
+    double f = 1.0, c = 1.0, d = 0.0;
+
+    //TODO LEA: Use Helper::ContinuedFractionEvaluation
+    int i, m;
+    for (i = 0; i <= 200; ++i) {
+        m = i/2;
+
+        double numerator;
+        if (i == 0) {
+            numerator = 1.0; /*First numerator is 1.0.*/
+        } else if (i % 2 == 0) {
+            numerator = (m*(b-m)*x)/((a+2.0*m-1.0)*(a+2.0*m)); /*Even term.*/
+        } else {
+            numerator = -((a+m)*(a+b+m)*x)/((a+2.0*m)*(a+2.0*m+1)); /*Odd term.*/
+        }
+
+        /*Do an iteration of Lentz's algorithm.*/
+        d = 1.0 + numerator * d;
+        if (std::fabs(d) < TINY) d = TINY;
+        d = 1.0 / d;
+
+        c = 1.0 + numerator / c;
+        if (std::fabs(c) < TINY) c = TINY;
+
+        const double cd = c*d;
+        f *= cd;
+
+        /*Check for stop.*/
+        if (std::fabs(1.0-cd) < STOP) {
+            return front * (f-1.0);
+        }
+    }
+
+    return NAN; /*Needed more loops, did not converge.*/
+}

apps/probability/law/incomplete_beta_function.h

@@ -0,0 +1,34 @@
+#ifndef PROBABILITY_INCOMPLETE_BETA_FUNCTION_H
+#define PROBABILITY_INCOMPLETE_BETA_FUNCTION_H
+
+/*
+ * zlib License
+ *
+ * Regularized Incomplete Beta Function
+ *
+ * Copyright (c) 2016, 2017 Lewis Van Winkle
+ * http://CodePlea.com
+ *
+ * This software is provided 'as-is', without any express or implied
+ * warranty. In no event will the authors be held liable for any damages
+ * arising from the use of this software.
+ *
+ * Permission is granted to anyone to use this software for any purpose,
+ * including commercial applications, and to alter it and redistribute it
+ * freely, subject to the following restrictions:
+ *
+ * 1. The origin of this software must not be misrepresented; you must not
+ *    claim that you wrote the original software. If you use this software
+ *    in a product, an acknowledgement in the product documentation would be
+ *    appreciated but is not required.
+ * 2. Altered source versions must be plainly marked as such, and must not be
+ *    misrepresented as being the original software.
+ * 3. This notice may not be removed or altered from any source distribution.
+ */
+
+
+// WARNING: this code has been modified
+
+double IncompleteBetaFunction(double a, double b, double x);
+
+#endif

apps/probability/law/regularized_gamma.cpp

@@ -67,7 +67,7 @@ bool regularizedGamma(double s, double x, double epsilon, int maxNumberOfIterati
   double infiniteSeriesValue = 0.0;
   if (!Helper::InfiniteSeriesEvaluation(
         1.0/s,
-        [](double previousTerm, double index, double s, double x) { return previousTerm * x / (s + index); },
+        [](double previousTerm, double index, double s, double x, double d1, double d2) { return previousTerm * x / (s + index); },
         epsilon,
         maxNumberOfIterations,
         &infiniteSeriesValue,

apps/probability/law/student_law.cpp

@@ -1,6 +1,7 @@
 #include "student_law.h"
+#include "incomplete_beta_function.h"
+#include "helper.h"
 #include <cmath>
-#include <float.h>
 
 namespace Probability {
 
@@ -26,8 +27,13 @@ bool StudentLaw::authorizedValueAtIndex(float x, int index) const {
 }
 
 double StudentLaw::cumulativeDistributiveFunctionAtAbscissa(double x) const {
-  return 0;
-  //TODO
+  if (x == 1) {
+    return 0.5;
+  }
+  const float k = m_parameter1;
+  const double sqrtXSquaredPlusK = std::sqrt(x*x + k);
+  double t = (x + sqrtXSquaredPlusK) / (2.0 * sqrtXSquaredPlusK);
+  return IncompleteBetaFunction(k/2.0, k/2.0, t);
 }
 
 double StudentLaw::cumulativeDistributiveInverseForProbability(double * probability) {

apps/probability/law/student_law.h

@@ -2,11 +2,15 @@
 #define PROBABILITY_STUDENT_LAW_H
 
 #include "one_parameter_law.h"
+#include <float.h>
 
 namespace Probability {
 
 class StudentLaw : public OneParameterLaw {
 public:
+  static constexpr int k_maxHypergeometricFunctionIterations = 1000; // TODO LEA factorize with Chi Squared
+  static constexpr double k_hypergeometricFunctionPrecision = DBL_EPSILON;
+
   StudentLaw() : OneParameterLaw(1.0f) {}
   I18n::Message title() override { return I18n::Message::StudentLaw; }
   Type type() const override { return Type::Student; }