Code Coverage
 
Lines
Functions and Methods
Classes and Traits
Total
89.80% covered (warning)
89.80%
44 / 49
80.00% covered (warning)
80.00%
4 / 5
 CRAP
0.00% covered (danger)
0.00%
0 / 1
Beta
89.80% covered (warning)
89.80%
44 / 49
80.00% covered (warning)
80.00%
4 / 5
 21.47
0.00% covered (danger)
0.00%
0 / 1
 __construct
n/a
0 / 0
n/a
0 / 0
 1
 incompleteBeta
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
 regularizedBeta
100.00% covered (success)
100.00%
10 / 10
100.00% covered (success)
100.00%
1 / 1
 7
 betaFraction
85.29% covered (warning)
85.29%
29 / 34
0.00% covered (danger)
0.00%
0 / 1
 7.16
 logBeta
100.00% covered (success)
100.00%
3 / 3
100.00% covered (success)
100.00%
1 / 1
 4
 beta
100.00% covered (success)
100.00%
1 / 1
100.00% covered (success)
100.00%
1 / 1
 1
1<?php
2/**
3 * Jingga
4 *
5 * PHP Version 8.1
6 *
7 * @package   phpOMS\Math\Functions
8 * @copyright Dennis Eichhorn
9 * @license   OMS License 2.0
10 * @version   1.0.0
11 * @link      https://jingga.app
12 */
13declare(strict_types=1);
14
15namespace phpOMS\Math\Functions;
16
17/**
18 * Beta function
19 *
20 * @package phpOMS\Math\Functions
21 * @license OMS License 2.0
22 * @link    https://jingga.app
23 * @since   1.0.0
24 */
25final class Beta
26{
27    /**
28     * Constructor.
29     *
30     * @since 1.0.0
31     * @codeCoverageIgnore
32     */
33    private function __construct()
34    {
35    }
36
37    /**
38     * Incomplete beta function
39     *
40     * @param float $x Value
41     * @param float $p p
42     * @param float $q q
43     *
44     * @return float
45     *
46     * @since 1.0.0
47     */
48    public static function incompleteBeta(float $x, float $p, float $q) : float
49    {
50        return self::regularizedBeta($x, $p, $q) * self::beta($p, $q);
51    }
52
53    /**
54     * Incomplete beta function
55     *
56     * @param float $x Value
57     * @param float $p p
58     * @param float $q q
59     *
60     * @return float
61     *
62     * @since 1.0.0
63     */
64    public static function regularizedBeta(float $x, float $p, float $q) : float
65    {
66        if ($x <= 0.0) {
67            return 0.0;
68        } elseif ($x >= 1.0) {
69            return 1.0;
70        } elseif ($p <= 0.0 || $q <= 0.0 || $p + $q > 10000000000.0) {
71            return 0.0;
72        }
73
74        $bGamma = \exp(-self::logBeta($p, $q) + $p * \log($x) + $q * \log(1.0 - $x));
75
76        // this uses the symmetry of the beta function
77        return ($x < ($p + 1.0) / ($p + $q + 2.0)
78            ? $bGamma * self::betaFraction($x, $p, $q) / $p
79            : 1.0 - $bGamma * self::betaFraction(1.0 - $x, $q, $p) / $q);
80    }
81
82    /**
83     * Fraction of the beta function
84     *
85     * @param float $x Value
86     * @param float $p p
87     * @param float $q q
88     *
89     * @see JSci
90     * @author Jaco van Kooten
91     * @license LGPL 2.1
92     * @since 1.0.0
93     */
94    private static function betaFraction(float $x, float $p, float $q) : float
95    {
96        $c      = 1.0;
97        $pqSum  = $p + $q;
98        $pPlus  = $p + 1.0;
99        $pMinus = $p - 1.0;
100        $h      = 1.0 - $pqSum * $x / $pPlus;
101
102        if (\abs($h) < 1.18e-37) {
103            $h = 1.18e-37;
104        }
105
106        $h     = 1.0 / $h;
107        $frac  = $h;
108        $m     = 1;
109        $delta = 0.0;
110
111        do {
112            $m2 = 2 * $m;
113            $d  = $m * ($q - $m) * $x / (($pMinus + $m2) * ($p + $m2));
114            $h  = 1.0 + $d * $h;
115            if (\abs($h) < 1.18e-37) {
116                $h = 1.18e-37;
117            }
118
119            $h = 1.0 / $h;
120            $c = 1.0 + $d / $c;
121            if (\abs($c) < 1.18e-37) {
122                $c = 1.18e-37;
123            }
124
125            $frac *= $h * $c;
126            $d     = -($p + $m) * ($pqSum + $m) * $x / (($p + $m2) * ($pPlus + $m2));
127            $h     = 1.0 + $d * $h;
128            if (\abs($h) < 1.18e-37) {
129                $h = 1.18e-37;
130            }
131
132            $h = 1.0 / $h;
133            $c = 1.0 + $d / $c;
134            if (\abs($c) < 1.18e-37) {
135                $c = 1.18e-37;
136            }
137
138            $delta = $h * $c;
139            $frac *= $delta;
140            ++$m;
141        } while ($m < 1000000 && \abs($delta - 1.0) > 8.88e-16);
142
143        return $frac;
144    }
145
146    /**
147     * Log of Beta
148     *
149     * @param float $p p
150     * @param float $q q
151     *
152     * @return float
153     *
154     * @since 1.0.0
155     */
156    public static function logBeta(float $p, float $q) : float
157    {
158        return $p <= 0.0 || $q <= 0.0 || $p + $q > 10000000000.0
159            ? 0.0
160            : Gamma::logGamma($p) + Gamma::logGamma($q) - Gamma::logGamma($p + $q);
161    }
162
163    /**
164     * Beta
165     *
166     * @param float $p p
167     * @param float $q q
168     *
169     * @return float
170     *
171     * @since 1.0.0
172     */
173    public static function beta(float $p, float $q) : float
174    {
175        return \exp(self::logBeta($p, $q));
176    }
177}