Fix misinterpreted EF in log-N deviate formulas

The two-sided formulas
were a wrong approach for computing the distribution with error factors.
The formulas are corrected
to work with one-sided, lower-tailed bounds,
as provided in the MEF standard.

As a result of the correction,
the sample standard deviations of the Monte-Carlo simulations increased,
and the simulated probabilities became closer to the mean values.

In addition, the maximum value evaluation now computes 99.9% percentile
instead of 99% (just for the sake of simplicity of 3 sigmas).

Fixes #147

src/expression/random_deviate.cc

@@ -108,14 +108,12 @@ double LogNormalDeviate::GetSample() noexcept {
 double LogNormalDeviate::Max() noexcept {
   double sigma = ComputeScale(level_.Mean(), ef_.Mean());
   double mu = ComputeLocation(mean_.Max(), sigma);
-  return std::exp(
-      std::sqrt(2) * std::pow(boost::math::erfc(1 / 50), -1) * sigma + mu);
+  return std::exp(3 * sigma + mu);
 }
 
 double LogNormalDeviate::ComputeScale(double level, double ef) noexcept {
-  double p = level + (1 - level) / 2;
-  double z = std::sqrt(2) * boost::math::erfc_inv(2 * p);
-  return std::log(ef) / std::abs(z);
+  double z = -std::sqrt(2) * boost::math::erfc_inv(2 * level);
+  return std::log(ef) / z;
 }
 
 double LogNormalDeviate::ComputeLocation(double mean, double sigma) noexcept {

src/expression/random_deviate.h

@@ -93,7 +93,8 @@ class NormalDeviate : public RandomDeviate {
   Expression& sigma_;  ///< Standard deviation of normal distribution.
 };
 
-/// Log-normal distribution.
+/// Log-normal distribution defined by
+/// its expected value and error factor of certain confidence level.
 class LogNormalDeviate : public RandomDeviate {
  public:
   /// Setup for log-normal distribution.
@@ -105,7 +106,7 @@ class LogNormalDeviate : public RandomDeviate {
   ///                  sigma is the scale factor.
   ///                  E(x) = exp(mu + sigma^2 / 2)
   /// @param[in] ef  The error factor of the log-normal distribution.
-  ///                EF = exp(z * sigma)
+  ///                EF = exp(z_alpha * sigma)
   /// @param[in] level  The confidence level.
   LogNormalDeviate(const ExpressionPtr& mean, const ExpressionPtr& ef,
                    const ExpressionPtr& level);
@@ -115,7 +116,7 @@ class LogNormalDeviate : public RandomDeviate {
 
   double Mean() noexcept override { return mean_.Mean(); }
 
-  /// 99 percentile estimate.
+  /// 99.9 percentile estimate.
   double Max() noexcept override;
 
   double Min() noexcept override { return 0; }

tests/bench_bscu_tests.cc

@@ -48,12 +48,12 @@ TEST_P(RiskAnalysisTest, BSCU) {
 
   if (settings.approximation() == "rare-event") {
     EXPECT_NEAR(0.135372, p_total(), 1e-4);
-    EXPECT_NEAR(0.1448961, mean(), 5e-3);
-    EXPECT_NEAR(0.203192, sigma(), 5e-3);
+    EXPECT_NEAR(0.137, mean(), 5e-3);
+    EXPECT_NEAR(0.217, sigma(), 5e-3);
   } else {
     EXPECT_NEAR(0.1124087, p_total(), 1e-4);
-    EXPECT_NEAR(0.1212541, mean(), 5e-3);
-    EXPECT_NEAR(0.1646726, sigma(), 5e-3);
+    EXPECT_NEAR(0.117, mean(), 5e-3);
+    EXPECT_NEAR(0.183, sigma(), 5e-3);
   }
 }
 

tests/bench_small_tree_tests.cc

@@ -34,8 +34,15 @@ TEST_P(RiskAnalysisTest, SmallTree) {
   std::set<std::set<std::string>> mcs = {{"e1", "e2"}, {"e3", "e4"}};
   EXPECT_EQ(2, products().size());
   EXPECT_EQ(mcs, products());
-  EXPECT_NEAR(0.02569, mean(), 1e-3);
-  EXPECT_NEAR(0.018065, sigma(), 2e-3);
+  if (settings.approximation() == "rare-event") {
+    EXPECT_NEAR(0.02696, p_total(), 1e-5);
+    EXPECT_NEAR(0.0255, mean(), 1e-3);
+    EXPECT_NEAR(0.0225, sigma(), 2e-3);
+  } else {
+    EXPECT_NEAR(0.02678, p_total(), 1e-5);
+    EXPECT_NEAR(0.0253, mean(), 1e-3);
+    EXPECT_NEAR(0.022, sigma(), 2e-3);
+  }
 }
 
 }  // namespace test