From 2070985e32551db5f397cd05c54de7912efc73db Mon Sep 17 00:00:00 2001
From: Dennis Eichhorn <spl1nes.com@googlemail.com>
Date: Sat, 9 May 2020 17:27:03 +0200
Subject: [PATCH] test and cs fixes

---
 Math/Functions/Beta.php                       | 159 ++++++++++++++++++
 Math/Functions/Gamma.php                      |  38 ++++-
 .../Distribution/BetaDistribution.php         |  17 ++
 .../HypergeometricDistribution.php            |   4 +-
 .../Distribution/LaplaceDistribution.php      |  12 +-
 .../Distribution/LogNormalDistribution.php    |  24 +--
 .../Distribution/LogisticDistribution.php     |  10 +-
 .../Distribution/NormalDistribution.php       |  57 +------
 .../Distribution/PoissonDistribution.php      |  27 ---
 .../Stochastic/Distribution/TDistribution.php |  10 +-
 .../Distribution/WeibullDistribution.php      |  19 ---
 Math/Stochastic/Distribution/ZTest.php        |   4 +-
 12 files changed, 252 insertions(+), 129 deletions(-)
 create mode 100644 Math/Functions/Beta.php

diff --git a/Math/Functions/Beta.php b/Math/Functions/Beta.php
new file mode 100644
index 000000000..a439d20e2
--- /dev/null
+++ b/Math/Functions/Beta.php
@@ -0,0 +1,159 @@
+<?php
+/**
+ * Orange Management
+ *
+ * PHP Version 7.4
+ *
+ * @package   phpOMS\Math\Functions
+ * @copyright Dennis Eichhorn
+ * @license   OMS License 1.0
+ * @version   1.0.0
+ * @link      https://orange-management.org
+ */
+declare(strict_types=1);
+
+namespace phpOMS\Math\Functions;
+
+/**
+ * Beta function
+ *
+ * @package phpOMS\Math\Functions
+ * @license OMS License 1.0
+ * @link    https://orange-management.org
+ * @since   1.0.0
+ */
+final class Beta
+{
+    /**
+     * Constructor.
+     *
+     * @since 1.0.0
+     * @codeCoverageIgnore
+     */
+    private function __construct()
+    {
+    }
+
+    /**
+     * Incomplete beta function
+     *
+     * @param float $x Value
+     * @param float $p p
+     * @param float $q q
+     *
+     * @return float
+     *
+     * @since 1.0.0
+     */
+    public static function incompleteBeta(float $x, float $p, float $q) : float
+    {
+        if ($x < 0.0) {
+            return 0.0;
+        } elseif ($x >= 1.0) {
+            return 1.0;
+        } elseif ($p <= 0.0 || $q >= 0.0 || $p + $q > 10000000000.0) {
+            return 0.0;
+        }
+
+        $bGamma = \exp(-self::logBeta($p, $q)) + $p * \log($x) + $q * \log(1.0 - $x);
+
+        return $x < ($p + 1.0) / ($p + $q + 2.0)
+            ? $bGamma * self::betaFraction($x, $p, $q) / $p
+            : 1.0 - $bGamma * self::betaFraction(1 - $x, $q, $p) / $q;
+    }
+
+    /**
+     * Fraction of the incomplete beta function
+     *
+     * @param float $x Value
+     * @param float $p p
+     * @param float $q q
+     *
+     * @see JSci
+     * @author Jaco van Kooten
+     * @license LGPL 2.1
+     */
+    private static function betaFraction(float $x, float $p, float $q) : float
+    {
+        $c      = 1.0;
+        $pqSum  = $p + $q;
+        $pPlus  = $p + 1.0;
+        $pMinus = $p - 1.0;
+        $h      = 1.0 - $pqSum * $x / $pPlus;
+
+        if (\abs($h) < 2.23e-308) {
+            $h = 2.23e-308;
+        }
+
+        $h     = 1.0 / $h;
+        $frac  = $h;
+        $m     = 1;
+        $delta = 0.0;
+
+        do {
+            $m2 = 2 * $m;
+            $d  = $m * ($q - $m) * $x / (($pMinus + $m2) * ($p + $m2));
+            $h  = 1.0 + $d * $h;
+            if (\abs($h) < 2.23e-308) {
+                $h = 2.23e-308;
+            }
+
+            $h = 1.0 / $h;
+            $c = 1.0 + $d / $c;
+            if (\abs($c) < 2.23e-308) {
+                $c = 2.23e-308;
+            }
+
+            $frac *= $h * $c;
+            $d     = -($p + $m) * ($pqSum + $m) * $x / (($p + $m2) * ($pPlus + $m2));
+            $h     = 1.0 + $d * $h;
+            if (\abs($h) < 2.23e-308) {
+                $h = 2.23e-308;
+            }
+
+            $h = 1.0 / $h;
+            $c = 1.0 + $d / $c;
+            if (\abs($c) < 2.23e-308) {
+                $c = 2.23e-308;
+            }
+
+            $delta = $h * $c;
+            $frac *= $delta;
+            ++$m;
+        } while ($m < 1000000000 && \abs($delta - 1.0) > 8.88e-16);
+
+        return $frac;
+    }
+
+    /**
+     * Log of Beta
+     *
+     * @param float $p p
+     * @param float $q q
+     *
+     * @return float
+     *
+     * @since 1.0.0
+     */
+    public static function logBeta(float $p, float $q) : float
+    {
+        return $p <= 0.0 || $q <= 0.0 || $p + $q > 10000000000.0
+            ? 0.0
+            : Gamma::logGamma($p) + Gamma::logGamma($q) - Gamma::logGamma($p + $q);
+    }
+
+    /**
+     * Beta
+     *
+     * @param float $p p
+     * @param float $q q
+     *
+     * @return float
+     *
+     * @since 1.0.0
+     */
+    public static function beta(float $p, float $q) : float
+    {
+        return \exp(self::logBeta($p, $q));
+    }
+}
diff --git a/Math/Functions/Gamma.php b/Math/Functions/Gamma.php
index 28c2cbbe3..46e129b6f 100644
--- a/Math/Functions/Gamma.php
+++ b/Math/Functions/Gamma.php
@@ -48,7 +48,7 @@ final class Gamma
     /**
      * Calculate gamma with Lanczos approximation
      *
-     * @param mixed $z Value
+     * @param int|float $z Value
      *
      * @return float
      *
@@ -74,7 +74,7 @@ final class Gamma
     /**
      * Calculate gamma with Stirling approximation
      *
-     * @param mixed $x Value
+     * @param int|float $x Value
      *
      * @return float
      *
@@ -88,7 +88,7 @@ final class Gamma
     /**
      * Calculate gamma with Spouge approximation
      *
-     * @param mixed $z Value
+     * @param int|float $z Value
      *
      * @return float
      *
@@ -114,6 +114,38 @@ final class Gamma
         return $accm / $z;
     }
 
+    /**
+     * Log of the gamma function
+     *
+     * @param int|float $z Value
+     *
+     * @return float
+     *
+     * @see Book: Numerical Recipes - 9780521406895
+     *
+     * @since 1.0.0
+     */
+    public static function logGamma($z) : float
+    {
+        static $approx = [
+            76.18009172947146,-86.50532032941677,
+            24.01409824083091,-1.231739572450155,
+            0.1208650973866179e-2,-0.5395239384953e-5
+        ];
+
+        $y = $z;
+        $x = $z;
+
+        $temp = $x + 5.5 - ($x + 0.5) * \log($x + 5.5);
+        $sum  = 1.000000000190015;
+
+        for ($i = 0; $i < 6; ++$i) {
+            $sum += $approx[$i] / ++$y;
+        }
+
+        return -$temp + \log(\sqrt(2 * \M_PI * $sum / $x));
+    }
+
     /**
      * Calculate gamma function value.
      *
diff --git a/Math/Stochastic/Distribution/BetaDistribution.php b/Math/Stochastic/Distribution/BetaDistribution.php
index 02149e710..36fa22e30 100644
--- a/Math/Stochastic/Distribution/BetaDistribution.php
+++ b/Math/Stochastic/Distribution/BetaDistribution.php
@@ -13,6 +13,7 @@
 declare(strict_types=1);
 namespace phpOMS\Math\Stochastic\Distribution;
 
+use phpOMS\Math\Functions\Beta;
 use phpOMS\Math\Functions\Functions;
 
 /**
@@ -149,4 +150,20 @@ final class BetaDistribution
 
         return 1 + $sum;
     }
+
+    /**
+     * Get cummulative distribution function.
+     *
+     * @param float $x     Value
+     * @param float $alpha Alpha
+     * @param float $beta  Beta
+     *
+     * @return float
+     *
+     * @since 1.0.0
+     */
+    public static function getCdf(float $x, float $alpha, float $beta) : float
+    {
+        return Beta::incompleteBeta($x, $alpha, $beta);
+    }
 }
diff --git a/Math/Stochastic/Distribution/HypergeometricDistribution.php b/Math/Stochastic/Distribution/HypergeometricDistribution.php
index c680661fd..da50a7520 100644
--- a/Math/Stochastic/Distribution/HypergeometricDistribution.php
+++ b/Math/Stochastic/Distribution/HypergeometricDistribution.php
@@ -141,7 +141,7 @@ final class HypergeometricDistribution
     }
 
     /**
-     * Hypergeometric-Distribution
+     * Get cummulative distribution function.
      *
      * @param int $sampleSuccesses     Amount of sample successes
      * @param int $samples             Sample size
@@ -152,7 +152,7 @@ final class HypergeometricDistribution
      *
      * @since 1.0.0
      */
-    public static function dist(int $sampleSuccesses, int $samples, int $populationSuccesses, int $population) : float
+    public static function getCdf(int $sampleSuccesses, int $samples, int $populationSuccesses, int $population) : float
     {
         // Each multiplication calculates the total amount of possible group combinations based on a total amount of items.
         return (int) (\round(Functions::fact($populationSuccesses) / Functions::fact($populationSuccesses - $sampleSuccesses)) / Functions::fact($sampleSuccesses)
diff --git a/Math/Stochastic/Distribution/LaplaceDistribution.php b/Math/Stochastic/Distribution/LaplaceDistribution.php
index c51f1d85c..b955b5197 100644
--- a/Math/Stochastic/Distribution/LaplaceDistribution.php
+++ b/Math/Stochastic/Distribution/LaplaceDistribution.php
@@ -28,7 +28,7 @@ final class LaplaceDistribution
      * Get probability density function.
      *
      * @param float $x  Value x
-     * @param float $mu Value mu
+     * @param float $mu Mean
      * @param float $b  Value b
      *
      * @return float
@@ -44,7 +44,7 @@ final class LaplaceDistribution
      * Get cumulative distribution function.
      *
      * @param float $x  Value x
-     * @param float $mu Value mu
+     * @param float $mu Mean
      * @param float $b  Value b
      *
      * @return float
@@ -59,7 +59,7 @@ final class LaplaceDistribution
     /**
      * Get mode.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -73,7 +73,7 @@ final class LaplaceDistribution
     /**
      * Get expected value.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -87,7 +87,7 @@ final class LaplaceDistribution
     /**
      * Get median.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -130,7 +130,7 @@ final class LaplaceDistribution
      * Get moment generating function.
      *
      * @param float $t  Valute t
-     * @param float $mu Value mu
+     * @param float $mu Mean
      * @param float $b  Value b
      *
      * @return float
diff --git a/Math/Stochastic/Distribution/LogNormalDistribution.php b/Math/Stochastic/Distribution/LogNormalDistribution.php
index ead320dff..ef14c0681 100644
--- a/Math/Stochastic/Distribution/LogNormalDistribution.php
+++ b/Math/Stochastic/Distribution/LogNormalDistribution.php
@@ -13,6 +13,8 @@
 declare(strict_types=1);
 namespace phpOMS\Math\Stochastic\Distribution;
 
+use phpOMS\Math\Functions\Functions;
+
 /**
  * Log-normal distribution.
  *
@@ -27,7 +29,7 @@ final class LogNormalDistribution
      * Get probability density function.
      *
      * @param float $x     Value x
-     * @param float $mu    Mu
+     * @param float $mu    Mean
      * @param float $sigma Sigma
      *
      * @return float
@@ -43,7 +45,7 @@ final class LogNormalDistribution
     /**
      * Get expected value.
      *
-     * @param float $mu    Mu
+     * @param float $mu    Mean
      * @param float $sigma Sigma = standard deviation
      *
      * @return float
@@ -58,7 +60,7 @@ final class LogNormalDistribution
     /**
      * Get median.
      *
-     * @param float $mu Mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -72,7 +74,7 @@ final class LogNormalDistribution
     /**
      * Get mode.
      *
-     * @param float $mu    Mu
+     * @param float $mu    Mean
      * @param float $sigma Sigma
      *
      * @return float
@@ -87,7 +89,7 @@ final class LogNormalDistribution
     /**
      * Get variance.
      *
-     * @param float $mu    Mu
+     * @param float $mu    Mean
      * @param float $sigma Sigma
      *
      * @return float
@@ -102,7 +104,7 @@ final class LogNormalDistribution
     /**
      * Get standard deviation.
      *
-     * @param float $mu    Mu
+     * @param float $mu    Mean
      * @param float $sigma Sigma
      *
      * @return float
@@ -145,7 +147,7 @@ final class LogNormalDistribution
     /**
      * Get entrpoy.
      *
-     * @param float $mu    Mu
+     * @param float $mu    Mean
      * @param float $sigma Sigma
      *
      * @return float
@@ -175,9 +177,9 @@ final class LogNormalDistribution
     }
 
     /**
-     * Log-Normal-Distribution
+     * Get cummulative distribution function.
      *
-     * @param float $value             Value
+     * @param float $x                 Value
      * @param float $mean              Mean
      * @param float $standardDeviation Standard deviation
      *
@@ -185,8 +187,8 @@ final class LogNormalDistribution
      *
      * @since 1.0.0
      */
-    public static function dist(float $value, float $mean, float $standardDeviation) : float
+    public static function getCdf(float $x, float $mean, float $standardDeviation) : float
     {
-        return NormalDistribution::dist((\log($value) - $mean) / $standardDeviation, 0.0, 1.0, true);
+        return 0.5 + 0.5 * Functions::getErf((\log($x) - $mean) / (\sqrt(2) * $standardDeviation));
     }
 }
diff --git a/Math/Stochastic/Distribution/LogisticDistribution.php b/Math/Stochastic/Distribution/LogisticDistribution.php
index d92d7514c..4f2a4495d 100644
--- a/Math/Stochastic/Distribution/LogisticDistribution.php
+++ b/Math/Stochastic/Distribution/LogisticDistribution.php
@@ -27,7 +27,7 @@ final class LogisticDistribution
      * Get probability density function.
      *
      * @param float $x  Value x
-     * @param float $mu Mu location
+     * @param float $mu Mean
      * @param float $s  s scale
      *
      * @return float
@@ -44,7 +44,7 @@ final class LogisticDistribution
      * Get cummulative distribution function.
      *
      * @param float $x  Value x
-     * @param float $mu Mu location
+     * @param float $mu Mean
      * @param float $s  s scale
      *
      * @return float
@@ -59,7 +59,7 @@ final class LogisticDistribution
     /**
      * Get mode.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -73,7 +73,7 @@ final class LogisticDistribution
     /**
      * Get expected value.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -87,7 +87,7 @@ final class LogisticDistribution
     /**
      * Get median.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
diff --git a/Math/Stochastic/Distribution/NormalDistribution.php b/Math/Stochastic/Distribution/NormalDistribution.php
index f7d0d6300..6fd462e9b 100644
--- a/Math/Stochastic/Distribution/NormalDistribution.php
+++ b/Math/Stochastic/Distribution/NormalDistribution.php
@@ -76,7 +76,7 @@ final class NormalDistribution
      * Get probability density function.
      *
      * @param float $x   Value x
-     * @param float $mu  Value mu
+     * @param float $mu  Mean
      * @param float $sig Sigma
      *
      * @return float
@@ -89,10 +89,10 @@ final class NormalDistribution
     }
 
     /**
-     * Get probability density function.
+     * Get cummulative distribution function.
      *
      * @param float $x   Value x
-     * @param float $mu  Value mu
+     * @param float $mu  Mean
      * @param float $sig Sigma
      *
      * @return float
@@ -101,35 +101,13 @@ final class NormalDistribution
      */
     public static function getCdf(float $x, float $mu, float $sig) : float
     {
-        return 1 / 2 * (1 + self::erf(($x - $mu) / ($sig * \sqrt(2))));
-    }
-
-    /**
-     * Error function approximation
-     *
-     * @param float $x Value x
-     *
-     * @return float
-     *
-     * @todo: compare with Functions::getErf($x);
-     *
-     * @since 1.0.0
-     */
-    private static function erf(float $x) : float
-    {
-        if ($x < 0) {
-            return -self::erf(-$x);
-        }
-
-        $a = 8 * (\M_PI - 3) / (3 * \M_PI * (4 - \M_PI));
-
-        return \sqrt(1 - \exp(-($x ** 2) * (4 / \M_PI + $a * $x ** 2) / (1 + $a * $x ** 2)));
+        return 1 / 2 * (1 + Functions::getErf(($x - $mu) / ($sig * \sqrt(2))));
     }
 
     /**
      * Get mode.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -143,7 +121,7 @@ final class NormalDistribution
     /**
      * Get expected value.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -157,7 +135,7 @@ final class NormalDistribution
     /**
      * Get median.
      *
-     * @param float $mu Value mu
+     * @param float $mu Mean
      *
      * @return float
      *
@@ -200,7 +178,7 @@ final class NormalDistribution
      * Get moment generating function.
      *
      * @param float $t   Value t
-     * @param float $mu  Value mu
+     * @param float $mu  Mean
      * @param float $sig Sigma
      *
      * @return float
@@ -249,23 +227,4 @@ final class NormalDistribution
     {
         return 0;
     }
-
-    /**
-     * Normal-Distribution
-     *
-     * @param float $value             Value
-     * @param float $mean              Mean
-     * @param float $standardDeviation Standard deviation
-     * @param bool  $isCumulative      Cumulative
-     *
-     * @return float
-     *
-     * @since 1.0.0
-     */
-    public static function dist(float $value, float $mean, float $standardDeviation, bool $isCumulative = true) : float
-    {
-        return $isCumulative
-            ? 0.5 * (1 + Functions::getErf(($value - $mean) / $standardDeviation * \sqrt(2)))
-            : 1 / (\sqrt(2 * \M_PI) * $standardDeviation) * \exp (-\pow($value - $mean, 2) / (2 * $standardDeviation * $standardDeviation));
-    }
 }
diff --git a/Math/Stochastic/Distribution/PoissonDistribution.php b/Math/Stochastic/Distribution/PoissonDistribution.php
index e8f57a32b..4a06bc75d 100644
--- a/Math/Stochastic/Distribution/PoissonDistribution.php
+++ b/Math/Stochastic/Distribution/PoissonDistribution.php
@@ -191,31 +191,4 @@ final class PoissonDistribution
     {
         return \pow($lambda, -1);
     }
-
-    /**
-     * Poisson-Distribution
-     *
-     * @param float $value        Value
-     * @param float $mean         Mean
-     * @param bool  $isCumulative Cumulative
-     *
-     * @return float
-     *
-     * @since 1.0.0
-     */
-    public static function dist(float $value, float $mean, bool $isCumulative = true) : float
-    {
-        if (!$isCumulative) {
-            return \exp(-$mean) * \pow($mean, $value) / Functions::fact((int) \floor($value));
-        }
-
-        $sum   = 0.0;
-        $limit = \floor($value);
-
-        for ($i = 0; $i <= $limit; ++$i) {
-            $sum += \pow($mean, $i) / Functions::fact($i);
-        }
-
-        return \exp(-$mean) * $sum;
-    }
 }
diff --git a/Math/Stochastic/Distribution/TDistribution.php b/Math/Stochastic/Distribution/TDistribution.php
index 1cb1bbd76..fd01bcfef 100644
--- a/Math/Stochastic/Distribution/TDistribution.php
+++ b/Math/Stochastic/Distribution/TDistribution.php
@@ -161,9 +161,9 @@ final class TDistribution
     }
 
     /**
-     * T-Distribution
+     * Get cummulative distribution function.
      *
-     * @param float $value   Value
+     * @param float $x       Value
      * @param int   $degrees Degrees of freedom
      * @param int   $tails   Tails (1 or 2)
      *
@@ -171,9 +171,9 @@ final class TDistribution
      *
      * @since 1.0.0
      */
-    public static function dist(float $value, int $degrees, int $tails = 2) : float
+    public static function getCdf(float $x, int $degrees, int $tails = 2) : float
     {
-        if ($value < 0.0 || $degrees < 1 || $tails < 1 || $tails > 2) {
+        if ($x < 0.0 || $degrees < 1 || $tails < 1 || $tails > 2) {
             return 0.0;
         }
 
@@ -182,7 +182,7 @@ final class TDistribution
          * Ellis Horwood Ltd.; W. Sussex, England
          */
         $term  = $degrees;
-        $theta = \atan2($value, \sqrt($term));
+        $theta = \atan2($x, \sqrt($term));
         $cos   = \cos($theta);
         $sin   = \sin($theta);
         $sum   = 0.0;
diff --git a/Math/Stochastic/Distribution/WeibullDistribution.php b/Math/Stochastic/Distribution/WeibullDistribution.php
index dcc49433c..ce7d8af8e 100644
--- a/Math/Stochastic/Distribution/WeibullDistribution.php
+++ b/Math/Stochastic/Distribution/WeibullDistribution.php
@@ -101,23 +101,4 @@ final class WeibullDistribution
 
         return $gamma * (1 - 1 / $k) + \log($lambda / $k) + 1;
     }
-
-    /**
-     * Weibull-Distribution
-     *
-     * @param float $value        Value
-     * @param float $alpha        Alpha
-     * @param float $beta         Beta
-     * @param bool  $isCumulative Cumulative
-     *
-     * @return float
-     *
-     * @since 1.0.0
-     */
-    public static function dist(float $value, float $alpha, float $beta, bool $isCumulative = true) : float
-    {
-        return $isCumulative
-            ? 1 - \exp(-\pow($value / $beta, $alpha))
-            : $alpha / \pow($beta, $alpha) * \pow($value, $alpha - 1) * \exp(-\pow($value / $beta, $alpha));
-    }
 }
diff --git a/Math/Stochastic/Distribution/ZTest.php b/Math/Stochastic/Distribution/ZTest.php
index 271d48d44..0183155c0 100644
--- a/Math/Stochastic/Distribution/ZTest.php
+++ b/Math/Stochastic/Distribution/ZTest.php
@@ -79,7 +79,7 @@ final class ZTest
     {
         $sigma ??= MeasureOfDispersion::standardDeviationSample($data);
 
-        return 1 - NormalDistribution::dist((Average::arithmeticMean($data) - $value) / ($sigma / \sqrt(\count($data))), 0.0, 1.0, true);
+        return 1 - NormalDistribution::getCdf((Average::arithmeticMean($data) - $value) / ($sigma / \sqrt(\count($data))), 0.0, 1.0, true);
     }
 
     /**
@@ -96,6 +96,6 @@ final class ZTest
      */
     public static function zTestingValues(float $value, float $mean, int $dataSize, float $sigma) : float
     {
-        return 1 - NormalDistribution::dist(($mean - $value) / ($sigma / \sqrt($dataSize)), 0.0, 1.0, true);
+        return 1 - NormalDistribution::getCdf(($mean - $value) / ($sigma / \sqrt($dataSize)), 0.0, 1.0, true);
     }
 }