Poisson::knuthRandom() - Code Metrics - JordanRL/Fermat - Measure and Improve Code Quality continuously with Scrutinizer

Poisson::knuthRandom() A
last analyzed 2022-10-22 04:55 UTC

↳ Parent: Poisson

Complexity

Conditions	2
Paths	1

Size

Total Lines	17
Code Lines	9

Duplication

Lines	0
Ratio	0 %

Code Coverage

Tests	0
CRAP Score	6

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	2
dl	0
loc	17
ccs	0
cts	0
cp	0
crap	6
rs	9.9666
c	1
b	0
f	0
eloc	9
nc	1
nop	0

<?php

namespace Samsara\Fermat\Stats\Values\Distribution;

use Samsara\Exceptions\SystemError\LogicalError\IncompatibleObjectState;
use Samsara\Exceptions\UsageError\IntegrityConstraint;
use Samsara\Exceptions\UsageError\OptionalExit;
use Samsara\Fermat\Core\Numbers;
use Samsara\Fermat\Stats\Types\Distribution;
use Samsara\Fermat\Core\Provider\RandomProvider;
use Samsara\Fermat\Core\Types\Base\Interfaces\Numbers\DecimalInterface;
use Samsara\Fermat\Core\Types\Base\Interfaces\Numbers\NumberInterface;
filter:
    dependency_paths: ["lib/*"]
use Samsara\Fermat\Core\Values\ImmutableDecimal;

class Poisson extends Distribution
{

    /**
     * @var ImmutableDecimal
     */
    private $lambda;

    /**
     * Poisson constructor.
     *
     * @param int|float|DecimalInterface $lambda
     *
     * @throws IntegrityConstraint
     */
    public function __construct($lambda)
    {
        $lambda = Numbers::makeOrDont(Numbers::IMMUTABLE, $lambda);

        if (!$lambda->isPositive()) {
            throw new IntegrityConstraint(
                'Lambda must be positive',
                'Provide a positive lambda',
                'Poisson distributions work on time to occurrence; the mean time to occurrence (lambda) must be positive'
            );
        }

        $this->lambda = $lambda;
    }

    /**
     * @param int|float|DecimalInterface $k
     *
     * @return ImmutableDecimal
     * @throws IntegrityConstraint
     * @throws IncompatibleObjectState
     */
    public function probabilityOfKEvents($k, int $scale = 10): ImmutableDecimal
    {

        return $this->pmf($k, $scale);
        
    }

    /**
     * @param int|float|DecimalInterface $x
     *
     * @return ImmutableDecimal
     * @throws IntegrityConstraint
     * @throws IncompatibleObjectState
     */
    public function cdf(int|float|DecimalInterface $x, int $scale = 10): ImmutableDecimal
    {

        $x = Numbers::makeOrDont(Numbers::IMMUTABLE, $x);

        if (!$x->isNatural()) {
            throw new IntegrityConstraint(
                'Only integers are valid x values for Poisson distributions',
                'Provide an integer value to calculate the CDF',
                'Poisson distributions describe discrete occurrences; only integers are valid x values'
            );
        }

        $internalScale = $scale + 2;

        $cumulative = Numbers::makeZero();

        for ($i = 0;$x->isGreaterThanOrEqualTo($i);$i++) {
            $cumulative =
                $cumulative->add(
                    $this->pmf(
                        $i,
                        $internalScale
                    )
                )->truncateToScale($scale);
        }

        return $cumulative;

    }

    /**
     * @param float|int|DecimalInterface $x
     *
     * @return ImmutableDecimal
     * @throws IntegrityConstraint
     * @throws IncompatibleObjectState
     */
    public function pmf(float|int|DecimalInterface $x, int $scale = 10): ImmutableDecimal
    {
        $x = Numbers::makeOrDont(Numbers::IMMUTABLE, $x);

        if (!$x->isNatural()) {
            throw new IntegrityConstraint(
                'Only integers are valid x values for Poisson distributions',
                'Provide an integer value to calculate the PMF',
                'Poisson distributions describe discrete occurrences; only integers are valid x values'
            );
        }

        $internalScale = $scale + 2;

        $e = Numbers::makeE($internalScale);

        /** @var ImmutableDecimal $pmf */
        $pmf =
            $this->lambda
            ->pow($x)
            ->multiply(
                $e->pow(
                    $this
                        ->lambda
                        ->multiply(-1)
                )
            )->divide(
                $x->factorial(),
                $internalScale
            )->truncateToScale($scale);

        return $pmf;
    }

    /**
     * @param int|float|DecimalInterface $x1
     * @param int|float|DecimalInterface $x2
     *
     * @return ImmutableDecimal
     * @throws IntegrityConstraint
     */
    public function rangePmf($x1, $x2): ImmutableDecimal
    {
        $x1 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x1);
        $x2 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x2);

        if ($x1->equals($x2)) {
            return Numbers::makeZero();
        } elseif ($x1->isGreaterThan($x2)) {
            $larger = $x1;
            $smaller = $x2;
        } else {
            $larger = $x2;
            $smaller = $x1;
        }

        if (!$larger->isNatural() || !$smaller->isNatural()) {
            throw new IntegrityConstraint(
                'Only integers are valid x values for Poisson distributions',
                'Provide integer values to calculate the range PMF',
                'Poisson distributions describe discrete occurrences; only integers are valid x values'
            );
        }

        $cumulative = Numbers::makeZero();

        for (;$larger->isGreaterThanOrEqualTo($smaller);$smaller = $smaller->add(1)) {
            $cumulative =
                $cumulative->add(
                    $this->pmf($smaller)
                );
        }

        return $cumulative;
    }

    /**
     * @return ImmutableDecimal
     * @throws IntegrityConstraint
     * @throws IncompatibleObjectState
     *
     * @codeCoverageIgnore
     */
    public function random(): ImmutableDecimal
    {
        if ($this->lambda->isLessThanOrEqualTo(30)) {
            return $this->knuthRandom();
        } else {
            return $this->methodPARandom();
        }
    }

    /**
     * WARNING: This function is of very limited use with Poisson distributions, and may represent a SIGNIFICANT
     * performance hit for certain values of $min, $max, $lambda, and $maxIterations
     *
     * @param int|float|NumberInterface $min
     * @param int|float|NumberInterface $max
     * @param int $maxIterations
     *
     * @return ImmutableDecimal
     * @throws OptionalExit
     * @throws IntegrityConstraint
     * @throws IncompatibleObjectState
     *
     * @codeCoverageIgnore
     */
    public function rangeRandom($min = 0, $max = PHP_INT_MAX, int $maxIterations = 20): ImmutableDecimal
    {
        $i = 0;

        do {
            $randomNumber = $this->random();
            $i++;
        } while (($randomNumber->isGreaterThan($max) || $randomNumber->isLessThan($min)) && $i < $maxIterations);

        if ($randomNumber->isGreaterThan($max) || $randomNumber->isLessThan($min)) {
            throw new OptionalExit(
                'All random numbers generated were outside of the requested range',
                'A suitable random number, restricted by the $max ('.$max.') and $min ('.$min.'), could not be found within '.$maxIterations.' iterations'
            );
        } else {
            return $randomNumber;
        }
    }

    /**
     * Method PA from The Computer Generation of Poisson Random Variables by A. C. Atkinson, 1979
     * Journal of the Royal Statistical Society Series C, Vol. 28, No. 1, Pages 29-35
     *
     * As described by John D. Cook: http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/
     *
     * @return ImmutableDecimal
     * @throws IntegrityConstraint
     * @throws IncompatibleObjectState
     *
     * @codeCoverageIgnore
     */
    protected function methodPARandom(): ImmutableDecimal
    {
        /** @var ImmutableDecimal $c */
        $c = $this->lambda->pow(-1)->multiply(3.36)->multiply(-1)->add(0.767);
        /** @var ImmutableDecimal $beta */
        $beta = Numbers::makePi()->divide($this->lambda->multiply(3)->sqrt());
        /** @var ImmutableDecimal $alpha */
        $alpha = $this->lambda->multiply($beta);
        /** @var ImmutableDecimal $k */
        $k = $c->ln(20)->subtract($this->lambda)->subtract($beta->ln(20));
        
        $one = Numbers::makeOne();
        /** @var ImmutableDecimal $oneHalf */
        $oneHalf = Numbers::make(Numbers::IMMUTABLE, '0.5');
        
        $e = Numbers::makeE();

        while (true) {
            /** @var ImmutableDecimal $u */
            $u = PolyfillProvider::randomInt(0, PHP_INT_MAX) / PHP_INT_MAX;
filter:
    dependency_paths: ["lib/*"]
            /** @var ImmutableDecimal $x */
            $x = $alpha->subtract($one->subtract($u)->divide($u)->ln(20)->divide($beta));
            /** @var ImmutableDecimal $n */
            $n = $x->add($oneHalf)->floor();

            if ($n->isNegative()) {
                continue;
            }

            /** @var ImmutableDecimal $v */
            $v = Numbers::make(Numbers::IMMUTABLE, PolyfillProvider::randomInt(0, PHP_INT_MAX))->divide(PHP_INT_MAX);
            /** @var ImmutableDecimal $y */
            $y = $alpha->subtract($beta->multiply($x));
            /** @var ImmutableDecimal $lhs */
            $lhs = $y->add($v->divide($e->pow($y)->add($one)->pow(2)));
            /** @var ImmutableDecimal $rhs */
            $rhs = $k->add($n->multiply($this->lambda->ln(20)))->subtract($n->factorial()->ln(20));

            if ($lhs->isLessThanOrEqualTo($rhs)) {
                return $n;
            }

            /*
             * At least attempt to free up some memory, since this particular method is extra hard on object instantiation
             */
            unset($u);
            unset($x);
            unset($n);
            unset($v);
            unset($y);
            unset($lhs);
            unset($rhs);
        }
    }

    /**
     * @return ImmutableDecimal
     * @throws IntegrityConstraint
     *
     * @codeCoverageIgnore
     */
    protected function knuthRandom(): ImmutableDecimal
    {
        /** @var ImmutableDecimal $L */
        $L = Numbers::makeE()->pow($this->lambda->multiply(-1));
        
        $k = Numbers::makeZero();
        
        $p = Numbers::makeOne();

        do {
            $k = $k->add(1);
            
            $u = RandomProvider::randomDecimal(10);
            $p = $p->multiply($u);
        } while ($p->isGreaterThan($L));

        return $k->subtract(1);
    }

}

1		<?php
2
3		namespace Samsara\Fermat\Stats\Values\Distribution;
4
5		use Samsara\Exceptions\SystemError\LogicalError\IncompatibleObjectState;
6		use Samsara\Exceptions\UsageError\IntegrityConstraint;
7		use Samsara\Exceptions\UsageError\OptionalExit;
8		use Samsara\Fermat\Core\Numbers;
9		use Samsara\Fermat\Stats\Types\Distribution;
10		use Samsara\Fermat\Core\Provider\RandomProvider;
11		use Samsara\Fermat\Core\Types\Base\Interfaces\Numbers\DecimalInterface;
12		use Samsara\Fermat\Core\Types\Base\Interfaces\Numbers\NumberInterface;
		0 ignored issues – show Bug introduced 2022-10-22 02:34 UTC by Report Bug Copy Issue Report The type `Samsara\Fermat\Core\Type...Numbers\NumberInterface` was not found. Maybe you did not declare it correctly or list all dependencies? The issue could also be caused by a filter entry in the build configuration. If the path has been excluded in your configuration, e.g. `excluded_paths: ["lib/"]`, you can move it to the dependency path list as follows: filter: dependency_paths: ["lib/"] For further information see https://scrutinizer-ci.com/docs/tools/php/php-scrutinizer/#list-dependency-paths Loading history...
13		use Samsara\Fermat\Core\Values\ImmutableDecimal;
14
15		class Poisson extends Distribution
16		{
17
18		/**
19		* @var ImmutableDecimal
20		*/
21		private $lambda;
22
23		/**
24		* Poisson constructor.
25		*
26		* @param int\|float\|DecimalInterface $lambda
27		*
28		* @throws IntegrityConstraint
29		*/
30	2	public function __construct($lambda)
31		{
32	2	$lambda = Numbers::makeOrDont(Numbers::IMMUTABLE, $lambda);
33
34	2	if (!$lambda->isPositive()) {
35		throw new IntegrityConstraint(
36		'Lambda must be positive',
37		'Provide a positive lambda',
38		'Poisson distributions work on time to occurrence; the mean time to occurrence (lambda) must be positive'
39		);
40		}
41
42	2	$this->lambda = $lambda;
43		}
44
45		/**
46		* @param int\|float\|DecimalInterface $k
47		*
48		* @return ImmutableDecimal
49		* @throws IntegrityConstraint
50		* @throws IncompatibleObjectState
51		*/
52		public function probabilityOfKEvents($k, int $scale = 10): ImmutableDecimal
53		{
54
55		return $this->pmf($k, $scale);
56
57		}
58
59		/**
60		* @param int\|float\|DecimalInterface $x
61		*
62		* @return ImmutableDecimal
63		* @throws IntegrityConstraint
64		* @throws IncompatibleObjectState
65		*/
66		public function cdf(int\|float\|DecimalInterface $x, int $scale = 10): ImmutableDecimal
67		{
68
69		$x = Numbers::makeOrDont(Numbers::IMMUTABLE, $x);
70
71		if (!$x->isNatural()) {
72		throw new IntegrityConstraint(
73		'Only integers are valid x values for Poisson distributions',
74		'Provide an integer value to calculate the CDF',
75		'Poisson distributions describe discrete occurrences; only integers are valid x values'
76		);
77		}
78
79		$internalScale = $scale + 2;
80
81		$cumulative = Numbers::makeZero();
82
83		for ($i = 0;$x->isGreaterThanOrEqualTo($i);$i++) {
84		$cumulative =
85		$cumulative->add(
86		$this->pmf(
87		$i,
88		$internalScale
89		)
90		)->truncateToScale($scale);
91		}
92
93		return $cumulative;
94
95		}
96
97		/**
98		* @param float\|int\|DecimalInterface $x
99		*
100		* @return ImmutableDecimal
101		* @throws IntegrityConstraint
102		* @throws IncompatibleObjectState
103		*/
104		public function pmf(float\|int\|DecimalInterface $x, int $scale = 10): ImmutableDecimal
105		{
106		$x = Numbers::makeOrDont(Numbers::IMMUTABLE, $x);
107
108		if (!$x->isNatural()) {
109		throw new IntegrityConstraint(
110		'Only integers are valid x values for Poisson distributions',
111		'Provide an integer value to calculate the PMF',
112		'Poisson distributions describe discrete occurrences; only integers are valid x values'
113		);
114		}
115
116		$internalScale = $scale + 2;
117
118		$e = Numbers::makeE($internalScale);
119
120		/** @var ImmutableDecimal $pmf */
121		$pmf =
122		$this->lambda
123		->pow($x)
124		->multiply(
125		$e->pow(
126		$this
127		->lambda
128		->multiply(-1)
129		)
130		)->divide(
131		$x->factorial(),
132		$internalScale
133		)->truncateToScale($scale);
134
135		return $pmf;
136		}
137
138		/**
139		* @param int\|float\|DecimalInterface $x1
140		* @param int\|float\|DecimalInterface $x2
141		*
142		* @return ImmutableDecimal
143		* @throws IntegrityConstraint
144		*/
145		public function rangePmf($x1, $x2): ImmutableDecimal
146		{
147		$x1 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x1);
148		$x2 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x2);
149
150		if ($x1->equals($x2)) {
151		return Numbers::makeZero();
152		} elseif ($x1->isGreaterThan($x2)) {
153		$larger = $x1;
154		$smaller = $x2;
155		} else {
156		$larger = $x2;
157		$smaller = $x1;
158		}
159
160		if (!$larger->isNatural() \|\| !$smaller->isNatural()) {
161		throw new IntegrityConstraint(
162		'Only integers are valid x values for Poisson distributions',
163		'Provide integer values to calculate the range PMF',
164		'Poisson distributions describe discrete occurrences; only integers are valid x values'
165		);
166		}
167
168		$cumulative = Numbers::makeZero();
169
170		for (;$larger->isGreaterThanOrEqualTo($smaller);$smaller = $smaller->add(1)) {
171		$cumulative =
172		$cumulative->add(
173		$this->pmf($smaller)
174		);
175		}
176
177		return $cumulative;
178		}
179
180		/**
181		* @return ImmutableDecimal
182		* @throws IntegrityConstraint
183		* @throws IncompatibleObjectState
184		*
185		* @codeCoverageIgnore
186		*/
187		public function random(): ImmutableDecimal
188		{
189		if ($this->lambda->isLessThanOrEqualTo(30)) {
190		return $this->knuthRandom();
191		} else {
192		return $this->methodPARandom();
193		}
194		}
195
196		/**
197		* WARNING: This function is of very limited use with Poisson distributions, and may represent a SIGNIFICANT
198		* performance hit for certain values of $min, $max, $lambda, and $maxIterations
199		*
200		* @param int\|float\|NumberInterface $min
201		* @param int\|float\|NumberInterface $max
202		* @param int $maxIterations
203		*
204		* @return ImmutableDecimal
205		* @throws OptionalExit
206		* @throws IntegrityConstraint
207		* @throws IncompatibleObjectState
208		*
209		* @codeCoverageIgnore
210		*/
211		public function rangeRandom($min = 0, $max = PHP_INT_MAX, int $maxIterations = 20): ImmutableDecimal
212		{
213		$i = 0;
214
215		do {
216		$randomNumber = $this->random();
217		$i++;
218		} while (($randomNumber->isGreaterThan($max) \|\| $randomNumber->isLessThan($min)) && $i < $maxIterations);
219
220		if ($randomNumber->isGreaterThan($max) \|\| $randomNumber->isLessThan($min)) {
221		throw new OptionalExit(
222		'All random numbers generated were outside of the requested range',
223		'A suitable random number, restricted by the $max ('.$max.') and $min ('.$min.'), could not be found within '.$maxIterations.' iterations'
224		);
225		} else {
226		return $randomNumber;
227		}
228		}
229
230		/**
231		* Method PA from The Computer Generation of Poisson Random Variables by A. C. Atkinson, 1979
232		* Journal of the Royal Statistical Society Series C, Vol. 28, No. 1, Pages 29-35
233		*
234		* As described by John D. Cook: http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/
235		*
236		* @return ImmutableDecimal
237		* @throws IntegrityConstraint
238		* @throws IncompatibleObjectState
239		*
240		* @codeCoverageIgnore
241		*/
242		protected function methodPARandom(): ImmutableDecimal
243		{
244		/** @var ImmutableDecimal $c */
245		$c = $this->lambda->pow(-1)->multiply(3.36)->multiply(-1)->add(0.767);
246		/** @var ImmutableDecimal $beta */
247		$beta = Numbers::makePi()->divide($this->lambda->multiply(3)->sqrt());
248		/** @var ImmutableDecimal $alpha */
249		$alpha = $this->lambda->multiply($beta);
250		/** @var ImmutableDecimal $k */
251		$k = $c->ln(20)->subtract($this->lambda)->subtract($beta->ln(20));
252
253		$one = Numbers::makeOne();
254		/** @var ImmutableDecimal $oneHalf */
255		$oneHalf = Numbers::make(Numbers::IMMUTABLE, '0.5');
256
257		$e = Numbers::makeE();
258
259		while (true) {
260		/** @var ImmutableDecimal $u */
261		$u = PolyfillProvider::randomInt(0, PHP_INT_MAX) / PHP_INT_MAX;
		0 ignored issues – show Bug introduced 2022-10-19 08:05 UTC by Report Bug Copy Issue Report The type `Samsara\Fermat\Stats\Val...bution\PolyfillProvider` was not found. Maybe you did not declare it correctly or list all dependencies? The issue could also be caused by a filter entry in the build configuration. If the path has been excluded in your configuration, e.g. `excluded_paths: ["lib/"]`, you can move it to the dependency path list as follows: filter: dependency_paths: ["lib/"] For further information see https://scrutinizer-ci.com/docs/tools/php/php-scrutinizer/#list-dependency-paths Loading history...
262		/** @var ImmutableDecimal $x */
263		$x = $alpha->subtract($one->subtract($u)->divide($u)->ln(20)->divide($beta));
264		/** @var ImmutableDecimal $n */
265		$n = $x->add($oneHalf)->floor();
266
267		if ($n->isNegative()) {
268		continue;
269		}
270
271		/** @var ImmutableDecimal $v */
272		$v = Numbers::make(Numbers::IMMUTABLE, PolyfillProvider::randomInt(0, PHP_INT_MAX))->divide(PHP_INT_MAX);
273		/** @var ImmutableDecimal $y */
274		$y = $alpha->subtract($beta->multiply($x));
275		/** @var ImmutableDecimal $lhs */
276		$lhs = $y->add($v->divide($e->pow($y)->add($one)->pow(2)));
277		/** @var ImmutableDecimal $rhs */
278		$rhs = $k->add($n->multiply($this->lambda->ln(20)))->subtract($n->factorial()->ln(20));
279
280		if ($lhs->isLessThanOrEqualTo($rhs)) {
281		return $n;
282		}
283
284		/*
285		* At least attempt to free up some memory, since this particular method is extra hard on object instantiation
286		*/
287		unset($u);
288		unset($x);
289		unset($n);
290		unset($v);
291		unset($y);
292		unset($lhs);
293		unset($rhs);
294		}
295		}
296
297		/**
298		* @return ImmutableDecimal
299		* @throws IntegrityConstraint
300		*
301		* @codeCoverageIgnore
302		*/
303		protected function knuthRandom(): ImmutableDecimal
304		{
305		/** @var ImmutableDecimal $L */
306		$L = Numbers::makeE()->pow($this->lambda->multiply(-1));
307
308		$k = Numbers::makeZero();
309
310		$p = Numbers::makeOne();
311
312		do {
313		$k = $k->add(1);
314
315		$u = RandomProvider::randomDecimal(10);
316		$p = $p->multiply($u);
317		} while ($p->isGreaterThan($L));
318
319		return $k->subtract(1);
320		}
321
322		}

JordanRL / Fermat

Poisson::knuthRandom() A last analyzed 2022-10-22 04:55 UTC

Complexity

Size

Duplication

Code Coverage

Importance

Duplication Side-by-Side

Filter issues like

Poisson::knuthRandom() A
last analyzed 2022-10-22 04:55 UTC