From 83fda9a587b382ba2252b13a31a07c830792ebaf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?L=C3=A9a=20Saviot?= <lea.saviot@numworks.com>
Date: Wed, 21 Aug 2019 14:45:27 +0200
Subject: [PATCH] [apps/proba] Fix Chi Square computations

Because the values can be very small or very big, computations should
not be made sequentially, to prevent rounding errors.
For instace, for degrees of freedom = 70, coefficient() would return 0
event though the cumulativeDistributiveFunctionAtAbscissa was not 0.
---
 .../distribution/chi_squared_distribution.cpp | 40 ++++++++-----------
 .../distribution/chi_squared_distribution.h   |  2 -
 .../probability/distribution/distribution.cpp |  4 +-
 3 files changed, 19 insertions(+), 27 deletions(-)

diff --git a/apps/probability/distribution/chi_squared_distribution.cpp b/apps/probability/distribution/chi_squared_distribution.cpp
index 464c222b0..ce7769f53 100644
--- a/apps/probability/distribution/chi_squared_distribution.cpp
+++ b/apps/probability/distribution/chi_squared_distribution.cpp
@@ -11,41 +11,40 @@ float ChiSquaredDistribution::xMin() const {
 }
 
 float ChiSquaredDistribution::xMax() const {
-  assert(m_parameter1 != 0);
+  assert(m_parameter1 != 0.0f);
   return (m_parameter1 + 5.0f * std::sqrt(m_parameter1)) * (1.0f + k_displayRightMarginRatio);
 }
 
 float ChiSquaredDistribution::yMax() const {
   float result;
-  if (m_parameter1/2.0 <= 1 + FLT_EPSILON) {
+  if (m_parameter1/2.0f <= 1.0f + FLT_EPSILON) {
     result = 0.5f;
   } else {
-    result = evaluateAtAbscissa(m_parameter1 - 1.0) * 1.2;
+    result = evaluateAtAbscissa(m_parameter1 - 1.0f) * 1.2f;
   }
   return result * (1.0f + k_displayTopMarginRatio);
 }
 
 float ChiSquaredDistribution::evaluateAtAbscissa(float x) const {
-  if (x < 0) {
+  if (x < 0.0f) {
     return NAN;
   }
-  const float halfk = m_parameter1/2;
-  const float halfX = x/2;
-  return coefficient() * std::pow(halfX, halfk-1) * std::exp(-halfX);
+  const float halfk = m_parameter1/2.0f;
+  const float halfX = x/2.0f;
+  return std::exp(-lgamma(halfk) - halfX + (halfk-1.0f) * std::log(halfX)) / 2.0f;
 }
 
 bool ChiSquaredDistribution::authorizedValueAtIndex(float x, int index) const {
   assert(index == 0);
-  return x > 0 && x == (int)x;
+  return x > 0.0f && x == (float)((int)x);
 }
 
 double ChiSquaredDistribution::cumulativeDistributiveFunctionAtAbscissa(double x) const {
-  if (x < 0) {
-    return 0;
+  if (x < 0.0) {
+    return 0.0;
   }
-  const float halfk = m_parameter1/2;
-  double result = 0;
-  if (regularizedGamma(halfk, x/2, k_regularizedGammaPrecision, k_maxRegularizedGammaIterations, &result)) {
+  double result = 0.0;
+  if (regularizedGamma(m_parameter1/2.0, x/2.0, k_regularizedGammaPrecision, k_maxRegularizedGammaIterations, &result)) {
     return result;
   }
   return NAN;
@@ -70,19 +69,14 @@ double ChiSquaredDistribution::cumulativeDistributiveInverseForProbability(doubl
    *         => b^(k/2) > k/2 * 2^(k/2) * proba * gamma(k/2)
    *         => exp(k/2 * ln(b)) > k/2 * 2^(k/2) * proba * gamma(k/2)
    *         => b > exp(2/k * ln(k/2 * 2^(k/2) * proba * gamma(k/2))) */
-
   const double k = m_parameter1;
-  const double ceilKOver2 = std::ceil(k/2);
-  const double kOver2Minus1 = k/2.0 - 1.0;
-  double xmax = m_parameter1 > 2.0 ?
-    2 * *probability * std::exp(std::lgamma(ceilKOver2)) / (exp(-kOver2Minus1) * std::pow(kOver2Minus1, kOver2Minus1)) :
+  const double ceilKOver2 = std::ceil(k/2.0f);
+  const double kOver2Minus1 = k/2.0f - 1.0f;
+  double xmax = m_parameter1 > 2.0f ?
+    2.0 * *probability * std::exp(std::lgamma(ceilKOver2)) / (exp(-kOver2Minus1) * std::pow(kOver2Minus1, kOver2Minus1)) :
     30.0; // Ad hoc value
+  xmax = std::isnan(xmax) ? 1000000000.0 : xmax;
   return cumulativeDistributiveInverseForProbabilityUsingIncreasingFunctionRoot(probability, DBL_EPSILON, maxDouble(xMax(), xmax));
 }
 
-float ChiSquaredDistribution::coefficient() const {
-  const float halfk = m_parameter1/2;
-  return 1 / (2 * std::exp(std::lgamma(halfk)));
-}
-
 }
diff --git a/apps/probability/distribution/chi_squared_distribution.h b/apps/probability/distribution/chi_squared_distribution.h
index ebd6c05e1..319f96a2d 100644
--- a/apps/probability/distribution/chi_squared_distribution.h
+++ b/apps/probability/distribution/chi_squared_distribution.h
@@ -30,8 +30,6 @@ public:
   bool authorizedValueAtIndex(float x, int index) const override;
   double cumulativeDistributiveFunctionAtAbscissa(double x) const override;
   double cumulativeDistributiveInverseForProbability(double * probability) override;
-private:
-  float coefficient() const;
 };
 
 }
diff --git a/apps/probability/distribution/distribution.cpp b/apps/probability/distribution/distribution.cpp
index 0a6aabdad..3949b974b 100644
--- a/apps/probability/distribution/distribution.cpp
+++ b/apps/probability/distribution/distribution.cpp
@@ -129,10 +129,10 @@ double Distribution::evaluateAtDiscreteAbscissa(int k) const {
 
 double Distribution::cumulativeDistributiveInverseForProbabilityUsingIncreasingFunctionRoot(double * probability, double ax, double bx) {
   assert(ax < bx);
-  if (*probability >= 1) {
+  if (*probability >= 1.0) {
     return INFINITY;
   }
-  if (*probability <= 0) {
+  if (*probability <= 0.0) {
     return -INFINITY;
   }
   Poincare::Coordinate2D result = Poincare::Solver::IncreasingFunctionRoot(