Beta::betaFraction() - Code Metrics - Inspection of "getColumnIterator() and getRowIterator() synonyms..." - PHPOffice/PhpSpreadsheet - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — master (#3302)

by Mark

created 2023-01-19 11:59 UTC

Beta::betaFraction() B

↳ Parent: Beta

Complexity

Conditions	8
Paths	34

Size

Total Lines	45
Code Lines	34

Duplication

Lines	0
Ratio	0 %

Code Coverage

Tests	30
CRAP Score	8.1867

Importance

Changes

Metric	Value
cc	8
eloc	34
nc	34
nop	3
dl	0
loc	45
ccs	30
cts	35
cp	0.8571
crap	8.1867
rs	8.1315
c	0
b	0
f	0

<?php

namespace PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;

use PhpOffice\PhpSpreadsheet\Calculation\ArrayEnabled;
use PhpOffice\PhpSpreadsheet\Calculation\Exception;
use PhpOffice\PhpSpreadsheet\Calculation\Functions;
use PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;

class Beta
{
    use ArrayEnabled;

    private const MAX_ITERATIONS = 256;

    private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;

    private const XMININ = 2.23e-308;

    /**
     * BETADIST.
     *
     * Returns the beta distribution.
     *
     * @param mixed $value Float value at which you want to evaluate the distribution
     *                      Or can be an array of values
     * @param mixed $alpha Parameter to the distribution as a float
     *                      Or can be an array of values
     * @param mixed $beta Parameter to the distribution as a float
     *                      Or can be an array of values
     * @param mixed $rMin as an float
     *                      Or can be an array of values
     * @param mixed $rMax as an float
     *                      Or can be an array of values
     *
     * @return array|float|string
     *         If an array of numbers is passed as an argument, then the returned result will also be an array
     *            with the same dimensions
     */
    public static function distribution($value, $alpha, $beta, $rMin = 0.0, $rMax = 1.0)
    {
        if (is_array($value) || is_array($alpha) || is_array($beta) || is_array($rMin) || is_array($rMax)) {
            return self::evaluateArrayArguments([self::class, __FUNCTION__], $value, $alpha, $beta, $rMin, $rMax);
        }

        $rMin = $rMin ?? 0.0;
        $rMax = $rMax ?? 1.0;

        try {
            $value = DistributionValidations::validateFloat($value);
            $alpha = DistributionValidations::validateFloat($alpha);
            $beta = DistributionValidations::validateFloat($beta);
            $rMax = DistributionValidations::validateFloat($rMax);
            $rMin = DistributionValidations::validateFloat($rMin);
        } catch (Exception $e) {
            return $e->getMessage();
        }

        if ($rMin > $rMax) {
            $tmp = $rMin;
            $rMin = $rMax;
            $rMax = $tmp;
        }
        if (($value < $rMin) || ($value > $rMax) || ($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax)) {
            return ExcelError::NAN();
        }

        $value -= $rMin;
        $value /= ($rMax - $rMin);

        return self::incompleteBeta($value, $alpha, $beta);
    }

    /**
     * BETAINV.
     *
     * Returns the inverse of the Beta distribution.
     *
     * @param mixed $probability Float probability at which you want to evaluate the distribution
     *                      Or can be an array of values
     * @param mixed $alpha Parameter to the distribution as a float
     *                      Or can be an array of values
     * @param mixed $beta Parameter to the distribution as a float
     *                      Or can be an array of values
     * @param mixed $rMin Minimum value as a float
     *                      Or can be an array of values
     * @param mixed $rMax Maximum value as a float
     *                      Or can be an array of values
     *
     * @return array|float|string
     *         If an array of numbers is passed as an argument, then the returned result will also be an array
     *            with the same dimensions
     */
    public static function inverse($probability, $alpha, $beta, $rMin = 0.0, $rMax = 1.0)
    {
        if (is_array($probability) || is_array($alpha) || is_array($beta) || is_array($rMin) || is_array($rMax)) {
            return self::evaluateArrayArguments([self::class, __FUNCTION__], $probability, $alpha, $beta, $rMin, $rMax);
        }

        $rMin = $rMin ?? 0.0;
        $rMax = $rMax ?? 1.0;

        try {
            $probability = DistributionValidations::validateProbability($probability);
            $alpha = DistributionValidations::validateFloat($alpha);
            $beta = DistributionValidations::validateFloat($beta);
            $rMax = DistributionValidations::validateFloat($rMax);
            $rMin = DistributionValidations::validateFloat($rMin);
        } catch (Exception $e) {
            return $e->getMessage();
        }

        if ($rMin > $rMax) {
            $tmp = $rMin;
            $rMin = $rMax;
            $rMax = $tmp;
        }
        if (($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax) || ($probability <= 0.0)) {
            return ExcelError::NAN();
        }

        return self::calculateInverse($probability, $alpha, $beta, $rMin, $rMax);
    }

    /**
     * @return float|string
     */
    private static function calculateInverse(float $probability, float $alpha, float $beta, float $rMin, float $rMax)
    {
        $a = 0;
        $b = 2;
        $guess = ($a + $b) / 2;

        $i = 0;
        while ((($b - $a) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
            $guess = ($a + $b) / 2;
            $result = self::distribution($guess, $alpha, $beta);
            if (($result === $probability) || ($result === 0.0)) {
                $b = $a;
            } elseif ($result > $probability) {
                $b = $guess;
            } else {
                $a = $guess;
            }
        }

        if ($i === self::MAX_ITERATIONS) {
            return ExcelError::NA();
        }

        return round($rMin + $guess * ($rMax - $rMin), 12);
    }

    /**
     * Incomplete beta function.
     *
     * @author Jaco van Kooten
     * @author Paul Meagher
     *
     * The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
     *
     * @param float $x require 0<=x<=1
     * @param float $p require p>0
     * @param float $q require q>0
     *
     * @return float 0 if x<0, p<=0, q<=0 or p+q>2.55E305 and 1 if x>1 to avoid errors and over/underflow
     */
    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) > self::LOG_GAMMA_X_MAX_VALUE)) {
            return 0.0;
        }

        $beta_gam = exp((0 - self::logBeta($p, $q)) + $p * log($x) + $q * log(1.0 - $x));
        if ($x < ($p + 1.0) / ($p + $q + 2.0)) {
            return $beta_gam * self::betaFraction($x, $p, $q) / $p;
        }

        return 1.0 - ($beta_gam * self::betaFraction(1 - $x, $q, $p) / $q);
    }

    // Function cache for logBeta function
    /** @var float */
    private static $logBetaCacheP = 0.0;

    /** @var float */
    private static $logBetaCacheQ = 0.0;

    /** @var float */
    private static $logBetaCacheResult = 0.0;

    /**
     * The natural logarithm of the beta function.
     *
     * @param float $p require p>0
     * @param float $q require q>0
     *
     * @return float 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
     *
     * @author Jaco van Kooten
     */
    private static function logBeta(float $p, float $q): float
    {
        if ($p != self::$logBetaCacheP || $q != self::$logBetaCacheQ) {
            self::$logBetaCacheP = $p;
            self::$logBetaCacheQ = $q;
            if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) {
                self::$logBetaCacheResult = 0.0;
            } else {
                self::$logBetaCacheResult = Gamma::logGamma($p) + Gamma::logGamma($q) - Gamma::logGamma($p + $q);
            }
        }

        return self::$logBetaCacheResult;
    }

    /**
     * Evaluates of continued fraction part of incomplete beta function.
     * Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
     *
     * @author Jaco van Kooten
     */
    private static function betaFraction(float $x, float $p, float $q): float
    {
        $c = 1.0;
        $sum_pq = $p + $q;
        $p_plus = $p + 1.0;
        $p_minus = $p - 1.0;
        $h = 1.0 - $sum_pq * $x / $p_plus;
        if (abs($h) < self::XMININ) {
            $h = self::XMININ;
        }
        $h = 1.0 / $h;
        $frac = $h;
        $m = 1;
        $delta = 0.0;
        while ($m <= self::MAX_ITERATIONS && abs($delta - 1.0) > Functions::PRECISION) {
            $m2 = 2 * $m;
            // even index for d
            $d = $m * ($q - $m) * $x / (($p_minus + $m2) * ($p + $m2));
            $h = 1.0 + $d * $h;
            if (abs($h) < self::XMININ) {
                $h = self::XMININ;
            }
            $h = 1.0 / $h;
            $c = 1.0 + $d / $c;
            if (abs($c) < self::XMININ) {
                $c = self::XMININ;
            }
            $frac *= $h * $c;
            // odd index for d
            $d = -($p + $m) * ($sum_pq + $m) * $x / (($p + $m2) * ($p_plus + $m2));
            $h = 1.0 + $d * $h;
            if (abs($h) < self::XMININ) {
                $h = self::XMININ;
            }
            $h = 1.0 / $h;
            $c = 1.0 + $d / $c;
            if (abs($c) < self::XMININ) {
                $c = self::XMININ;
            }
            $delta = $h * $c;
            $frac *= $delta;
            ++$m;
        }

        return $frac;
    }

    /*
    private static function betaValue(float $a, float $b): float
    {
        return (Gamma::gammaValue($a) * Gamma::gammaValue($b)) /
            Gamma::gammaValue($a + $b);
    }

    private static function regularizedIncompleteBeta(float $value, float $a, float $b): float
    {
        return self::incompleteBeta($value, $a, $b) / self::betaValue($a, $b);
    }
    */
}


1		<?php
2
3		namespace PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions;
4
5		use PhpOffice\PhpSpreadsheet\Calculation\ArrayEnabled;
6		use PhpOffice\PhpSpreadsheet\Calculation\Exception;
7		use PhpOffice\PhpSpreadsheet\Calculation\Functions;
8		use PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError;
9
10		class Beta
11		{
12		use ArrayEnabled;
13
14		private const MAX_ITERATIONS = 256;
15
16		private const LOG_GAMMA_X_MAX_VALUE = 2.55e305;
17
18		private const XMININ = 2.23e-308;
19
20		/**
21		* BETADIST.
22		*
23		* Returns the beta distribution.
24		*
25		* @param mixed $value Float value at which you want to evaluate the distribution
26		* Or can be an array of values
27		* @param mixed $alpha Parameter to the distribution as a float
28		* Or can be an array of values
29		* @param mixed $beta Parameter to the distribution as a float
30		* Or can be an array of values
31		* @param mixed $rMin as an float
32		* Or can be an array of values
33		* @param mixed $rMax as an float
34		* Or can be an array of values
35		*
36		* @return array\|float\|string
37		* If an array of numbers is passed as an argument, then the returned result will also be an array
38		* with the same dimensions
39		*/
40	30	public static function distribution($value, $alpha, $beta, $rMin = 0.0, $rMax = 1.0)
41		{
42	30	if (is_array($value) \|\| is_array($alpha) \|\| is_array($beta) \|\| is_array($rMin) \|\| is_array($rMax)) {
43	21	return self::evaluateArrayArguments([self::class, __FUNCTION__], $value, $alpha, $beta, $rMin, $rMax);
44		}
45
46	30	$rMin = $rMin ?? 0.0;
47	30	$rMax = $rMax ?? 1.0;
48
49		try {
50	30	$value = DistributionValidations::validateFloat($value);
51	29	$alpha = DistributionValidations::validateFloat($alpha);
52	28	$beta = DistributionValidations::validateFloat($beta);
53	27	$rMax = DistributionValidations::validateFloat($rMax);
54	26	$rMin = DistributionValidations::validateFloat($rMin);
55	5	} catch (Exception $e) {
56	5	return $e->getMessage();
57		}
58
59	25	if ($rMin > $rMax) {
60	1	$tmp = $rMin;
61	1	$rMin = $rMax;
62	1	$rMax = $tmp;
63		}
64	25	if (($value < $rMin) \|\| ($value > $rMax) \|\| ($alpha <= 0) \|\| ($beta <= 0) \|\| ($rMin == $rMax)) {
65	7	return ExcelError::NAN();
66		}
67
68	18	$value -= $rMin;
69	18	$value /= ($rMax - $rMin);
70
71	18	return self::incompleteBeta($value, $alpha, $beta);
72		}
73
74		/**
75		* BETAINV.
76		*
77		* Returns the inverse of the Beta distribution.
78		*
79		* @param mixed $probability Float probability at which you want to evaluate the distribution
80		* Or can be an array of values
81		* @param mixed $alpha Parameter to the distribution as a float
82		* Or can be an array of values
83		* @param mixed $beta Parameter to the distribution as a float
84		* Or can be an array of values
85		* @param mixed $rMin Minimum value as a float
86		* Or can be an array of values
87		* @param mixed $rMax Maximum value as a float
88		* Or can be an array of values
89		*
90		* @return array\|float\|string
91		* If an array of numbers is passed as an argument, then the returned result will also be an array
92		* with the same dimensions
93		*/
94	22	public static function inverse($probability, $alpha, $beta, $rMin = 0.0, $rMax = 1.0)
95		{
96	22	if (is_array($probability) \|\| is_array($alpha) \|\| is_array($beta) \|\| is_array($rMin) \|\| is_array($rMax)) {
97	21	return self::evaluateArrayArguments([self::class, __FUNCTION__], $probability, $alpha, $beta, $rMin, $rMax);
98		}
99
100	22	$rMin = $rMin ?? 0.0;
101	22	$rMax = $rMax ?? 1.0;
102
103		try {
104	22	$probability = DistributionValidations::validateProbability($probability);
105	19	$alpha = DistributionValidations::validateFloat($alpha);
106	18	$beta = DistributionValidations::validateFloat($beta);
107	17	$rMax = DistributionValidations::validateFloat($rMax);
108	16	$rMin = DistributionValidations::validateFloat($rMin);
109	7	} catch (Exception $e) {
110	7	return $e->getMessage();
111		}
112
113	15	if ($rMin > $rMax) {
114	1	$tmp = $rMin;
115	1	$rMin = $rMax;
116	1	$rMax = $tmp;
117		}
118	15	if (($alpha <= 0) \|\| ($beta <= 0) \|\| ($rMin == $rMax) \|\| ($probability <= 0.0)) {
119	6	return ExcelError::NAN();
120		}
121
122	9	return self::calculateInverse($probability, $alpha, $beta, $rMin, $rMax);
123		}
124
125		/**
126		* @return float\|string
127		*/
128	9	private static function calculateInverse(float $probability, float $alpha, float $beta, float $rMin, float $rMax)
129		{
130	9	$a = 0;
131	9	$b = 2;
132	9	$guess = ($a + $b) / 2;
133
134	9	$i = 0;
135	9	while ((($b - $a) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
136	9	$guess = ($a + $b) / 2;
137	9	$result = self::distribution($guess, $alpha, $beta);
138	9	if (($result === $probability) \|\| ($result === 0.0)) {
139	1	$b = $a;
140	9	} elseif ($result > $probability) {
141	9	$b = $guess;
142		} else {
143	9	$a = $guess;
144		}
145		}
146
147	9	if ($i === self::MAX_ITERATIONS) {
148		return ExcelError::NA();
149		}
150
151	9	return round($rMin + $guess * ($rMax - $rMin), 12);
152		}
153
154		/**
155		* Incomplete beta function.
156		*
157		* @author Jaco van Kooten
158		* @author Paul Meagher
159		*
160		* The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
161		*
162		* @param float $x require 0<=x<=1
163		* @param float $p require p>0
164		* @param float $q require q>0
165		*
166		* @return float 0 if x<0, p<=0, q<=0 or p+q>2.55E305 and 1 if x>1 to avoid errors and over/underflow
167		*/
168	24	public static function incompleteBeta(float $x, float $p, float $q): float
169		{
170	24	if ($x <= 0.0) {
171		return 0.0;
172	24	} elseif ($x >= 1.0) {
173	9	return 1.0;
174	24	} elseif (($p <= 0.0) \|\| ($q <= 0.0) \|\| (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) {
175		return 0.0;
176		}
177
178	24	$beta_gam = exp((0 - self::logBeta($p, $q)) + $p * log($x) + $q * log(1.0 - $x));
179	24	if ($x < ($p + 1.0) / ($p + $q + 2.0)) {
180	12	return $beta_gam * self::betaFraction($x, $p, $q) / $p;
181		}
182
183	19	return 1.0 - ($beta_gam * self::betaFraction(1 - $x, $q, $p) / $q);
184		}
185
186		// Function cache for logBeta function
187		/** @var float */
188		private static $logBetaCacheP = 0.0;
189
190		/** @var float */
191		private static $logBetaCacheQ = 0.0;
192
193		/** @var float */
194		private static $logBetaCacheResult = 0.0;
195
196		/**
197		* The natural logarithm of the beta function.
198		*
199		* @param float $p require p>0
200		* @param float $q require q>0
201		*
202		* @return float 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
203		*
204		* @author Jaco van Kooten
205		*/
206	24	private static function logBeta(float $p, float $q): float
207		{
208	24	if ($p != self::$logBetaCacheP \|\| $q != self::$logBetaCacheQ) {
209	20	self::$logBetaCacheP = $p;
210	20	self::$logBetaCacheQ = $q;
211	20	if (($p <= 0.0) \|\| ($q <= 0.0) \|\| (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) {
212		self::$logBetaCacheResult = 0.0;
213		} else {
214	20	self::$logBetaCacheResult = Gamma::logGamma($p) + Gamma::logGamma($q) - Gamma::logGamma($p + $q);
215		}
216		}
217
218	24	return self::$logBetaCacheResult;
219		}
220
221		/**
222		* Evaluates of continued fraction part of incomplete beta function.
223		* Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
224		*
225		* @author Jaco van Kooten
226		*/
227	24	private static function betaFraction(float $x, float $p, float $q): float
228		{
229	24	$c = 1.0;
230	24	$sum_pq = $p + $q;
231	24	$p_plus = $p + 1.0;
232	24	$p_minus = $p - 1.0;
233	24	$h = 1.0 - $sum_pq * $x / $p_plus;
234	24	if (abs($h) < self::XMININ) {
235		$h = self::XMININ;
236		}
237	24	$h = 1.0 / $h;
238	24	$frac = $h;
239	24	$m = 1;
240	24	$delta = 0.0;
241	24	while ($m <= self::MAX_ITERATIONS && abs($delta - 1.0) > Functions::PRECISION) {
242	24	$m2 = 2 * $m;
243		// even index for d
244	24	$d = $m * ($q - $m) * $x / (($p_minus + $m2) * ($p + $m2));
245	24	$h = 1.0 + $d * $h;
246	24	if (abs($h) < self::XMININ) {
247		$h = self::XMININ;
248		}
249	24	$h = 1.0 / $h;
250	24	$c = 1.0 + $d / $c;
251	24	if (abs($c) < self::XMININ) {
252		$c = self::XMININ;
253		}
254	24	$frac = $h $c;
255		// odd index for d
256	24	$d = -($p + $m) * ($sum_pq + $m) * $x / (($p + $m2) * ($p_plus + $m2));
257	24	$h = 1.0 + $d * $h;
258	24	if (abs($h) < self::XMININ) {
259		$h = self::XMININ;
260		}
261	24	$h = 1.0 / $h;
262	24	$c = 1.0 + $d / $c;
263	24	if (abs($c) < self::XMININ) {
264		$c = self::XMININ;
265		}
266	24	$delta = $h * $c;
267	24	$frac *= $delta;
268	24	++$m;
269		}
270
271	24	return $frac;
272		}
273
274		/*
275		private static function betaValue(float $a, float $b): float
276		{
277		return (Gamma::gammaValue($a) * Gamma::gammaValue($b)) /
278		Gamma::gammaValue($a + $b);
279		}
280
281		private static function regularizedIncompleteBeta(float $value, float $a, float $b): float
282		{
283		return self::incompleteBeta($value, $a, $b) / self::betaValue($a, $b);
284		}
285		*/
286		}
287

PHPOffice / PhpSpreadsheet

Pull Request — master (#3302)

Beta::betaFraction() B

Complexity

Size

Duplication

Code Coverage

Importance

Duplication Side-by-Side

Filter issues like