diff --git a/MixtComp/src/lib/Mixture/Simple/Gaussian/GaussianLikelihood.cpp b/MixtComp/src/lib/Mixture/Simple/Gaussian/GaussianLikelihood.cpp index c44c6387f..ccf3ce111 100644 --- a/MixtComp/src/lib/Mixture/Simple/Gaussian/GaussianLikelihood.cpp +++ b/MixtComp/src/lib/Mixture/Simple/Gaussian/GaussianLikelihood.cpp @@ -71,7 +71,12 @@ Real GaussianLikelihood::lnObservedProbability(int i, int k) const { Real supCdf = normal_.cdf(supBound, mean, sd); - logProba = std::log(supCdf - infCdf); + + Real proba = supCdf - infCdf; + if(proba < epsilonProba) + proba = epsilonProba; + + logProba = std::log(proba); } break; @@ -80,6 +85,10 @@ Real GaussianLikelihood::lnObservedProbability(int i, int k) const { Real supCdf = normal_.cdf(supBound, mean, sd); + + if(supCdf < epsilonProba) + supCdf = epsilonProba; + logProba = std::log(supCdf); } break; @@ -89,7 +98,11 @@ Real GaussianLikelihood::lnObservedProbability(int i, int k) const { Real infCdf = normal_.cdf(infBound, mean, sd); - logProba = std::log(1. - infCdf); + Real proba = 1. - infCdf; + if(proba < epsilonProba) + proba = epsilonProba; + + logProba = std::log(proba); } break; diff --git a/MixtComp/src/lib/Mixture/Simple/Weibull/WeibullLikelihood.cpp b/MixtComp/src/lib/Mixture/Simple/Weibull/WeibullLikelihood.cpp index ccb45d7af..e7bc7a3fa 100644 --- a/MixtComp/src/lib/Mixture/Simple/Weibull/WeibullLikelihood.cpp +++ b/MixtComp/src/lib/Mixture/Simple/Weibull/WeibullLikelihood.cpp @@ -58,7 +58,11 @@ Real WeibullLikelihood::lnObservedProbability(Index i, Index k) const { case missingRUIntervals_: { Real infBound = augData_.misData_(i).second[0]; Real infCdf = weibull_.cdf(infBound, kParam, lambda); - logProba = std::log(1.0 - infCdf); + Real proba = 1. - infCdf; + if(proba < epsilonProba) + proba = epsilonProba; + + logProba = std::log(proba); } break; @@ -67,7 +71,12 @@ Real WeibullLikelihood::lnObservedProbability(Index i, Index k) const { Real supBound = augData_.misData_(i).second[1]; Real infCdf = weibull_.cdf(infBound, kParam, lambda); Real supCdf = weibull_.cdf(supBound, kParam, lambda); - logProba = std::log(supCdf - infCdf); + + Real proba = supCdf - infCdf; + if(proba < epsilonProba) + proba = epsilonProba; + + logProba = std::log(proba); } break; diff --git a/MixtComp/src/lib/Various/Constants.cpp b/MixtComp/src/lib/Various/Constants.cpp index e29316860..1108dc859 100644 --- a/MixtComp/src/lib/Various/Constants.cpp +++ b/MixtComp/src/lib/Various/Constants.cpp @@ -39,10 +39,12 @@ const std::string eol = "\n"; const int nbSamplingAttempts = 10000; -const Real epsilon = 1.e-8; -const std::string epsilonStr = "1.e-8"; +const Real epsilon = 1.e-10; +const std::string epsilonStr = "1.e-10"; const Real logEpsilon = std::log(epsilon); +const Real epsilonProba = std::numeric_limits::epsilon(); + #if defined MINMODALITY const int minModality = MINMODALITY; #else diff --git a/MixtComp/src/lib/Various/Constants.h b/MixtComp/src/lib/Various/Constants.h index 57ec18934..a39df3847 100644 --- a/MixtComp/src/lib/Various/Constants.h +++ b/MixtComp/src/lib/Various/Constants.h @@ -40,6 +40,8 @@ extern const int nbSamplingAttempts; // number of sampling attempts, when not en extern const Real epsilon; // very small value of real to check for near zero values extern const std::string epsilonStr; // previous value, in scientific notation extern const Real logEpsilon; // log of very small value +extern const Real epsilonProba; // very small value of real to check for near zero values for probabilities + extern const int minModality; // minimal modality for categorical models (for example, 0-based or 1-based numbering) extern const int minIndex; diff --git a/MixtComp/src/utest/Functional/UTestFuncCSComputation.cpp b/MixtComp/src/utest/Functional/UTestFuncCSComputation.cpp index 004537b70..8bbdd1437 100644 --- a/MixtComp/src/utest/Functional/UTestFuncCSComputation.cpp +++ b/MixtComp/src/utest/Functional/UTestFuncCSComputation.cpp @@ -72,8 +72,8 @@ TEST(FuncCSComputation, regressionNoNoise) { regression(design, y, betaEstimated, sdEstimated); - ASSERT_EQ(true, betaEstimated.isApprox(beta, epsilon)); - ASSERT_NEAR(0., sdEstimated, epsilon); + EXPECT_EQ(true, betaEstimated.isApprox(beta, 1e-8)); + EXPECT_NEAR(0., sdEstimated, 1e-8); } TEST(FuncCSComputation, regressionNoise) {