UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-01-19 00:37:25 +01:00
Code Issues Packages Projects Releases Wiki Activity

[apps/proba] Fix Student cumulativeDistributiveInverseForProbability

Browse Source
This commit is contained in:
Léa Saviot
2019-08-19 15:15:19 +02:00
parent b26cd6a4fd
commit 56543e8886
2 changed files with 26 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

23
apps/probability/law/chi_squared_law.cpp
View File

@@ -60,12 +60,23 @@ double ChiSquaredLaw::cumulativeDistributiveInverseForProbability(double * proba
* => regularizedGamma(halfk, b/2, k_regularizedGammaPrecision) >proba
* => int(exp(-t)*t^(k/2 - 1), t, 0, b/2)/gamma(k/2) > proba
* => int(exp(-t)*t^(k/2 - 1), t, 0, b/2) > proba * gamma(k/2)
* true if: (b/2 - 0) * exp(-(k/2 - 1))*(k/2 - 1)^(k/2 - 1) > proba * gamma(k/2)
* => b/2 * exp(-(k/2 - 1) (1 + ln(k/2 - 1))) > proba * gamma(k/2)
* => b > 2 * proba * gamma(k/2) / exp(-(k/2 - 1) (1 + ln(k/2 - 1))) */
const double ceilKOver2 = std::ceil(m_parameter1/2);
const double kOver2Minus1 = std::ceil(m_parameter1/2);
double xmax = 2 * *probability * std::exp(std::lgamma(ceilKOver2)) / exp(-(kOver2Minus1)*(1 + std::log(kOver2Minus1)));
*
* true if: for (k/2 - 1) > 0 which is k > 2
* (b/2 - 0) * exp(-(k/2 - 1))*(k/2 - 1)^(k/2 - 1) > proba * gamma(k/2)
* => b/2 * exp(-(k/2 - 1) (1 - ln(k/2 - 1))) > proba * gamma(k/2)
* => b > 2 * proba * gamma(k/2) / exp(-(k/2 - 1) (1 + ln(k/2 - 1)))
* else
* 2/k * (b/2)^(k/2) > proba * gamma(k/2)
* => 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)) :
std::exp(2/k * std::log(k * std::pow(2.0, kOver2Minus1) * *probability * std::exp(std::lgamma(k/2))));
return cumulativeDistributiveInverseForProbabilityUsingBrentRoots(probability, DBL_EPSILON, maxDouble(xMax(), xmax));
}

11
apps/probability/law/student_law.cpp
View File

@@ -27,7 +27,7 @@ bool StudentLaw::authorizedValueAtIndex(float x, int index) const {
}
double StudentLaw::cumulativeDistributiveFunctionAtAbscissa(double x) const {
if (x == 1) {
if (x == 0) {
return 0.5;
}
const float k = m_parameter1;
@@ -37,7 +37,14 @@ double StudentLaw::cumulativeDistributiveFunctionAtAbscissa(double x) const {
}
double StudentLaw::cumulativeDistributiveInverseForProbability(double * probability) {
return cumulativeDistributiveInverseForProbabilityUsingBrentRoots(probability, xMin(), xMax());
if (*probability == 0.5) {
return 0.0;
}
const double small = DBL_EPSILON;
const double big = m_parameter1 >= 1 ? 50.0 : 50 / m_parameter1;
double xmin = *probability < 0.5 ? -big : small;
double xmax = *probability < 0.5 ? -small : big;
return cumulativeDistributiveInverseForProbabilityUsingBrentRoots(probability, xmin, xmax);
}
float StudentLaw::coefficient() const {