From a280207449e9d0f066af3748cbc77c19eb292d92 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9a=20Saviot?= Date: Wed, 14 Aug 2019 09:53:05 +0200 Subject: [PATCH] [apps/proba] Student cumulativeDistributiveInverseForProbability --- apps/probability/Makefile | 4 + apps/probability/law/helper.cpp | 4 +- apps/probability/law/helper.h | 4 +- .../law/hypergeometric_function.cpp | 32 ++++++++ .../probability/law/hypergeometric_function.h | 20 +++++ .../law/incomplete_beta_function.cpp | 82 +++++++++++++++++++ .../law/incomplete_beta_function.h | 34 ++++++++ apps/probability/law/regularized_gamma.cpp | 2 +- apps/probability/law/student_law.cpp | 12 ++- apps/probability/law/student_law.h | 4 + 10 files changed, 190 insertions(+), 8 deletions(-) create mode 100644 apps/probability/law/hypergeometric_function.cpp create mode 100644 apps/probability/law/hypergeometric_function.h create mode 100644 apps/probability/law/incomplete_beta_function.cpp create mode 100644 apps/probability/law/incomplete_beta_function.h diff --git a/apps/probability/Makefile b/apps/probability/Makefile index 17573d725..887ba98fe 100644 --- a/apps/probability/Makefile +++ b/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\ diff --git a/apps/probability/law/helper.cpp b/apps/probability/law/helper.cpp index 8f6f565c1..9722ce9cf 100644 --- a/apps/probability/law/helper.cpp +++ b/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) { diff --git a/apps/probability/law/helper.h b/apps/probability/law/helper.h index 83f8924f8..c7f5b93fe 100644 --- a/apps/probability/law/helper.h +++ b/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); }; diff --git a/apps/probability/law/hypergeometric_function.cpp b/apps/probability/law/hypergeometric_function.cpp new file mode 100644 index 000000000..34a2ed59c --- /dev/null +++ b/apps/probability/law/hypergeometric_function.cpp @@ -0,0 +1,32 @@ +#include "hypergeometric_function.h" +#include "helper.h" +#include +#include +#include + + +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 +} diff --git a/apps/probability/law/hypergeometric_function.h b/apps/probability/law/hypergeometric_function.h new file mode 100644 index 000000000..5579387f8 --- /dev/null +++ b/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 + diff --git a/apps/probability/law/incomplete_beta_function.cpp b/apps/probability/law/incomplete_beta_function.cpp new file mode 100644 index 000000000..2911bc7a0 --- /dev/null +++ b/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 +#include + +#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.*/ +} diff --git a/apps/probability/law/incomplete_beta_function.h b/apps/probability/law/incomplete_beta_function.h new file mode 100644 index 000000000..f8a266862 --- /dev/null +++ b/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 diff --git a/apps/probability/law/regularized_gamma.cpp b/apps/probability/law/regularized_gamma.cpp index da210fd94..4ea23ef13 100644 --- a/apps/probability/law/regularized_gamma.cpp +++ b/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, diff --git a/apps/probability/law/student_law.cpp b/apps/probability/law/student_law.cpp index 3bbb14f14..b0300b756 100644 --- a/apps/probability/law/student_law.cpp +++ b/apps/probability/law/student_law.cpp @@ -1,6 +1,7 @@ #include "student_law.h" +#include "incomplete_beta_function.h" +#include "helper.h" #include -#include 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) { diff --git a/apps/probability/law/student_law.h b/apps/probability/law/student_law.h index 97dc2e519..c53fafe87 100644 --- a/apps/probability/law/student_law.h +++ b/apps/probability/law/student_law.h @@ -2,11 +2,15 @@ #define PROBABILITY_STUDENT_LAW_H #include "one_parameter_law.h" +#include 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; }