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

[apps/proba] Fix Chi Square computations

Browse Source 
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.
This commit is contained in:
Léa Saviot
2019-08-21 14:45:27 +02:00
parent 04b0df9a72
commit 83fda9a587
3 changed files with 19 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
apps/probability/distribution/chi_squared_distribution.cpp
View File

@@ -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)));
}
}

2
apps/probability/distribution/chi_squared_distribution.h
View File

@@ -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;
};
}

4
apps/probability/distribution/distribution.cpp
View File

@@ -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(