From 3dd5112a0b061e2790ab34f4b319ffc21c62711f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9a=20Saviot?= Date: Tue, 27 Aug 2019 10:52:47 +0200 Subject: [PATCH] [poincare] Use reg incomplete beta function in binomial distribution --- apps/probability/Makefile | 1 - .../distribution/binomial_distribution.cpp | 19 ++++------------- .../distribution/student_distribution.cpp | 4 ++-- poincare/Makefile | 1 + .../regularized_incomplete_beta_function.h | 9 +++++--- poincare/src/binomial_distribution.cpp | 21 +++---------------- .../regularized_incomplete_beta_function.cpp | 10 ++++++--- poincare/test/approximation.cpp | 2 +- 8 files changed, 24 insertions(+), 43 deletions(-) rename apps/probability/distribution/incomplete_beta_function.h => poincare/include/poincare/regularized_incomplete_beta_function.h (83%) rename apps/probability/distribution/incomplete_beta_function.cpp => poincare/src/regularized_incomplete_beta_function.cpp (90%) diff --git a/apps/probability/Makefile b/apps/probability/Makefile index 57f946a83..053ef6766 100644 --- a/apps/probability/Makefile +++ b/apps/probability/Makefile @@ -6,7 +6,6 @@ app_probability_test_src = $(addprefix apps/probability/,\ distribution/chi_squared_distribution.cpp \ distribution/geometric_distribution.cpp \ distribution/helper.cpp \ - distribution/incomplete_beta_function.cpp \ distribution/distribution.cpp \ distribution/regularized_gamma.cpp \ distribution/student_distribution.cpp \ diff --git a/apps/probability/distribution/binomial_distribution.cpp b/apps/probability/distribution/binomial_distribution.cpp index 8a9b386a3..ed1e8f213 100644 --- a/apps/probability/distribution/binomial_distribution.cpp +++ b/apps/probability/distribution/binomial_distribution.cpp @@ -54,22 +54,11 @@ float BinomialDistribution::yMax() const { bool BinomialDistribution::authorizedValueAtIndex(float x, int index) const { if (index == 0) { - /* As the cumulative probability are computed by looping over all discrete - * abscissa within the interesting range, the complexity of the cumulative - * probability is linear with the size of the range. Here we cap the maximal - * size of the range to 10000. If one day we want to increase or get rid of - * this cap, we should implement the explicit formula of the cumulative - * probability (which depends on an incomplete beta function) to make the - * comlexity O(1). */ - if (x != (int)x || x < 0.0f || x > 99999.0f) { - return false; - } - return true; + // n must be a positive integer + return (x == (int)x) && x >= 0.0f; } - if (x < 0.0f || x > 1.0f) { - return false; - } - return true; + // p must be between 0 and 1 + return (x >= 0.0f) && (x <= 1.0f); } double BinomialDistribution::cumulativeDistributiveInverseForProbability(double * probability) { diff --git a/apps/probability/distribution/student_distribution.cpp b/apps/probability/distribution/student_distribution.cpp index fd7c420d1..39c510552 100644 --- a/apps/probability/distribution/student_distribution.cpp +++ b/apps/probability/distribution/student_distribution.cpp @@ -1,5 +1,5 @@ #include "student_distribution.h" -#include "incomplete_beta_function.h" +#include #include "helper.h" #include @@ -36,7 +36,7 @@ double StudentDistribution::cumulativeDistributiveFunctionAtAbscissa(double x) c 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); + return Poincare::RegularizedIncompleteBetaFunction(k/2.0, k/2.0, t); } double StudentDistribution::cumulativeDistributiveInverseForProbability(double * probability) { diff --git a/poincare/Makefile b/poincare/Makefile index b71437ebe..aea534ca3 100644 --- a/poincare/Makefile +++ b/poincare/Makefile @@ -35,6 +35,7 @@ poincare_src += $(addprefix poincare/src/,\ exception_checkpoint.cpp \ helpers.cpp \ normal_distribution.cpp \ + regularized_incomplete_beta_function.cpp \ ) poincare_src += $(addprefix poincare/src/,\ diff --git a/apps/probability/distribution/incomplete_beta_function.h b/poincare/include/poincare/regularized_incomplete_beta_function.h similarity index 83% rename from apps/probability/distribution/incomplete_beta_function.h rename to poincare/include/poincare/regularized_incomplete_beta_function.h index f8a266862..3936c1892 100644 --- a/apps/probability/distribution/incomplete_beta_function.h +++ b/poincare/include/poincare/regularized_incomplete_beta_function.h @@ -1,5 +1,5 @@ -#ifndef PROBABILITY_INCOMPLETE_BETA_FUNCTION_H -#define PROBABILITY_INCOMPLETE_BETA_FUNCTION_H +#ifndef POINCARE_REGULARIZED_INCOMPLETE_BETA_FUNCTION_H +#define POINCARE_REGULARIZED_INCOMPLETE_BETA_FUNCTION_H /* * zlib License @@ -28,7 +28,10 @@ // WARNING: this code has been modified +namespace Poincare { -double IncompleteBetaFunction(double a, double b, double x); +double RegularizedIncompleteBetaFunction(double a, double b, double x); + +} #endif diff --git a/poincare/src/binomial_distribution.cpp b/poincare/src/binomial_distribution.cpp index 9b6585484..6c2e1f44d 100644 --- a/poincare/src/binomial_distribution.cpp +++ b/poincare/src/binomial_distribution.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -49,16 +50,8 @@ T BinomialDistribution::CumulativeDistributiveFunctionAtAbscissa(T x, T n, T p) if (x >= n) { return (T)1.0; } - bool isDouble = sizeof(T) == sizeof(double); - return Solver::CumulativeDistributiveFunctionForNDefinedFunction( - x, - [](double abscissa, Context * context, Poincare::Preferences::ComplexFormat complexFormat, Poincare::Preferences::AngleUnit angleUnit, const void * n, const void * p, const void * isDouble) { - if (*(bool *)isDouble) { - return (double)BinomialDistribution::EvaluateAtAbscissa(abscissa, *(reinterpret_cast(n)), *(reinterpret_cast(p))); - } - return (double)BinomialDistribution::EvaluateAtAbscissa(abscissa, *(reinterpret_cast(n)), *(reinterpret_cast(p))); - }, (Context *)nullptr, Preferences::ComplexFormat::Real, Preferences::AngleUnit::Degree, &n, &p, &isDouble); - // Context, complex format and angle unit are dummy values + T floorX = std::floor(x); + return RegularizedIncompleteBetaFunction(n-floorX, floorX+1.0, 1.0-p); } template @@ -100,14 +93,6 @@ T BinomialDistribution::CumulativeDistributiveInverseForProbability(T probabilit template bool BinomialDistribution::ParametersAreOK(T n, T p) { - /* TODO As the cumulative probability are computed by looping over all discrete - * abscissa within the interesting range, the complexity of the cumulative - * probability is linear with the size of the range. Here we cap the maximal - * size of the range to 10000. If one day we want to increase or get rid of - * this cap, we should implement the explicit formula of the cumulative - * probability (which depends on an incomplete beta function) to make the - * comlexity O(1). */ - //TODO LEA n doit être <= 99999.0f) { return !std::isnan(n) && !std::isnan(p) && !std::isinf(n) && !std::isinf(p) && (n == (int)n) && (n >= (T)0) diff --git a/apps/probability/distribution/incomplete_beta_function.cpp b/poincare/src/regularized_incomplete_beta_function.cpp similarity index 90% rename from apps/probability/distribution/incomplete_beta_function.cpp rename to poincare/src/regularized_incomplete_beta_function.cpp index 83b11928a..b7cc61179 100644 --- a/apps/probability/distribution/incomplete_beta_function.cpp +++ b/poincare/src/regularized_incomplete_beta_function.cpp @@ -25,19 +25,21 @@ // WARNING: this code has been modified -#include "incomplete_beta_function.h" +#include #include #include +namespace Poincare { + #define STOP 1.0e-8 #define TINY 1.0e-30 -double IncompleteBetaFunction(double a, double b, double x) { +double RegularizedIncompleteBetaFunction(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.*/ + return (1.0-RegularizedIncompleteBetaFunction(b,a,1.0-x)); /*Use the fact that beta is symmetrical.*/ } /*Find the first part before the continued fraction.*/ @@ -80,3 +82,5 @@ double IncompleteBetaFunction(double a, double b, double x) { return NAN; /*Needed more loops, did not converge.*/ } + +} diff --git a/poincare/test/approximation.cpp b/poincare/test/approximation.cpp index 83cc1c409..8d86303bb 100644 --- a/poincare/test/approximation.cpp +++ b/poincare/test/approximation.cpp @@ -235,7 +235,7 @@ QUIZ_CASE(poincare_approximation_function) { assert_expression_approximates_to("abs([[3+2𝐢,3+4𝐢][5+2𝐢,3+2𝐢]])", "[[3.605551275464,5][5.3851648071345,3.605551275464]]"); assert_expression_approximates_to("binomcdf(5.3, 9, 0.7)", "0.270341", Degree, Cartesian, 6); // FIXME: precision problem - assert_expression_approximates_to("binomcdf(5.3, 9, 0.7)", "0.270340902"); + assert_expression_approximates_to("binomcdf(5.3, 9, 0.7)", "0.270340902", Degree, Cartesian, 10); //FIXME precision problem assert_expression_approximates_to("binomial(10, 4)", "210"); assert_expression_approximates_to("binomial(10, 4)", "210");