1
|
|
|
<?php
|
2
|
|
|
|
3
|
|
|
namespace Samsara\Fermat\Provider\Stats\Distribution;
|
4
|
|
|
|
5
|
|
|
use RandomLib\Factory;
|
6
|
|
|
use Samsara\Exceptions\UsageError\IntegrityConstraint;
|
7
|
|
|
use Samsara\Exceptions\UsageError\OptionalExit;
|
8
|
|
|
use Samsara\Fermat\Numbers;
|
9
|
|
|
use Samsara\Fermat\Provider\Stats\Distribution\Base\DistributionInterface;
|
10
|
|
|
use Samsara\Fermat\Types\Base\NumberInterface;
|
11
|
|
|
use Samsara\Fermat\Values\ImmutableNumber;
|
12
|
|
|
|
13
|
|
|
class Poisson
|
14
|
|
|
{
|
15
|
|
|
|
16
|
|
|
/**
|
17
|
|
|
* @var ImmutableNumber
|
18
|
|
|
*/
|
19
|
|
|
private $lambda;
|
20
|
|
|
|
21
|
|
View Code Duplication |
public function __construct($lambda)
|
|
|
|
|
22
|
|
|
{
|
23
|
|
|
$lambda = Numbers::makeOrDont(Numbers::IMMUTABLE, $lambda);
|
24
|
|
|
|
25
|
|
|
if (!$lambda->isPositive()) {
|
26
|
|
|
throw new IntegrityConstraint(
|
27
|
|
|
'Lambda must be positive',
|
28
|
|
|
'Provide a positive lambda',
|
29
|
|
|
'Poisson distributions work on time to occurrence; the mean time to occurrence (lambda) must be positive'
|
30
|
|
|
);
|
31
|
|
|
}
|
32
|
|
|
|
33
|
|
|
$this->lambda = $lambda;
|
|
|
|
|
34
|
|
|
}
|
35
|
|
|
|
36
|
|
|
public function probabilityOfKEvents($k)
|
37
|
|
|
{
|
38
|
|
|
|
39
|
|
|
return $this->pmf($k);
|
40
|
|
|
|
41
|
|
|
}
|
42
|
|
|
|
43
|
|
|
public function cdf($x)
|
44
|
|
|
{
|
45
|
|
|
|
46
|
|
|
$x = Numbers::makeOrDont(Numbers::IMMUTABLE, $x);
|
47
|
|
|
|
48
|
|
|
if (!$x->isNatural()) {
|
49
|
|
|
throw new IntegrityConstraint(
|
50
|
|
|
'Only integers are valid x values for Poisson distributions',
|
51
|
|
|
'Provide an integer value to calculate the CDF',
|
52
|
|
|
'Poisson distributions describe discrete occurrences; only integers are valid x values'
|
53
|
|
|
);
|
54
|
|
|
}
|
55
|
|
|
|
56
|
|
|
if (
|
57
|
|
|
function_exists('stats_cdf_poisson') &&
|
58
|
|
|
$this->lambda->isLessThanOrEqualTo(PHP_INT_MAX) &&
|
59
|
|
|
$x->isLessThanOrEqualTo(PHP_INT_MAX)
|
60
|
|
|
) {
|
61
|
|
|
return Numbers::makeOrDont(Numbers::IMMUTABLE, stats_cdf_poisson($x->getValue(), $this->lambda->getValue(), 1));
|
62
|
|
|
}
|
63
|
|
|
|
64
|
|
|
$cumulative = Numbers::makeZero();
|
65
|
|
|
|
66
|
|
|
for ($i = 0;$x->isGreaterThanOrEqualTo($i);$i++) {
|
67
|
|
|
$cumulative = $cumulative->add($this->pmf($i));
|
68
|
|
|
}
|
69
|
|
|
|
70
|
|
|
return $cumulative;
|
71
|
|
|
|
72
|
|
|
}
|
73
|
|
|
|
74
|
|
|
/**
|
75
|
|
|
* Not sure if the usage of stats_dens_pmf_poisson() is correct here because it is undocumented in
|
76
|
|
|
* the C code for the PHP extension and on PHP.net
|
77
|
|
|
*
|
78
|
|
|
* @param $x
|
79
|
|
|
*
|
80
|
|
|
* @return ImmutableNumber
|
81
|
|
|
* @throws IntegrityConstraint
|
82
|
|
|
*/
|
83
|
|
View Code Duplication |
public function pmf($x)
|
|
|
|
|
84
|
|
|
{
|
85
|
|
|
$x = Numbers::makeOrDont(Numbers::IMMUTABLE, $x);
|
86
|
|
|
|
87
|
|
|
if (!$x->isNatural()) {
|
88
|
|
|
throw new IntegrityConstraint(
|
89
|
|
|
'Only integers are valid x values for Poisson distributions',
|
90
|
|
|
'Provide an integer value to calculate the PMF',
|
91
|
|
|
'Poisson distributions describe discrete occurrences; only integers are valid x values'
|
92
|
|
|
);
|
93
|
|
|
}
|
94
|
|
|
|
95
|
|
|
if (
|
96
|
|
|
function_exists('stats_dens_pmf_poisson') &&
|
97
|
|
|
$x->isLessThanOrEqualTo(PHP_INT_MAX) &&
|
98
|
|
|
$this->lambda->isLessThanOrEqualTo(PHP_INT_MAX)
|
99
|
|
|
) {
|
100
|
|
|
return Numbers::make(Numbers::IMMUTABLE, stats_dens_pmf_poisson($x->getValue(), $this->lambda->getValue()));
|
101
|
|
|
}
|
102
|
|
|
|
103
|
|
|
$e = Numbers::makeE();
|
104
|
|
|
|
105
|
|
|
/** @var ImmutableNumber $pmf */
|
106
|
|
|
$pmf = $this->lambda->pow($x)->multiply($e->pow($this->lambda->multiply(-1)))->divide($x->factorial());
|
107
|
|
|
|
108
|
|
|
return $pmf;
|
109
|
|
|
}
|
110
|
|
|
|
111
|
|
|
public function rangePmf($x1, $x2)
|
112
|
|
|
{
|
113
|
|
|
$x1 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x1);
|
114
|
|
|
$x2 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x2);
|
115
|
|
|
|
116
|
|
|
if ($x1->equals($x2)) {
|
117
|
|
|
return Numbers::makeZero();
|
118
|
|
|
} elseif ($x1->isGreaterThan($x2)) {
|
|
|
|
|
119
|
|
|
$larger = $x1;
|
120
|
|
|
$smaller = $x2;
|
121
|
|
|
} else {
|
122
|
|
|
$larger = $x2;
|
123
|
|
|
$smaller = $x1;
|
124
|
|
|
}
|
125
|
|
|
|
126
|
|
|
if (!$larger->isNatural() || !$smaller->isNatural()) {
|
|
|
|
|
127
|
|
|
throw new IntegrityConstraint(
|
128
|
|
|
'Only integers are valid x values for Poisson distributions',
|
129
|
|
|
'Provide integer values to calculate the range PMF',
|
130
|
|
|
'Poisson distributions describe discrete occurrences; only integers are valid x values'
|
131
|
|
|
);
|
132
|
|
|
}
|
133
|
|
|
|
134
|
|
|
$cumulative = Numbers::makeZero();
|
135
|
|
|
|
136
|
|
|
for (;$larger->isGreaterThanOrEqualTo($smaller);$smaller->add(1)) {
|
|
|
|
|
137
|
|
|
$cumulative = $cumulative->add($this->pmf($smaller));
|
138
|
|
|
}
|
139
|
|
|
|
140
|
|
|
return $cumulative;
|
141
|
|
|
}
|
142
|
|
|
|
143
|
|
|
/**
|
144
|
|
|
* @return ImmutableNumber
|
145
|
|
|
*/
|
146
|
|
|
public function random()
|
147
|
|
|
{
|
148
|
|
|
|
149
|
|
|
if (
|
150
|
|
|
function_exists('stats_gen_rand_ipoisson') &&
|
151
|
|
|
$this->lambda->isLessThanOrEqualTo(PHP_INT_MAX) // WTF are you doing with a lambda larger than PHP_INT_MAX!?
|
152
|
|
|
) {
|
153
|
|
|
return Numbers::make(Numbers::IMMUTABLE, stats_gen_rand_ipoisson($this->lambda->getValue()));
|
154
|
|
|
} elseif ($this->lambda->isLessThanOrEqualTo(30)) {
|
155
|
|
|
return $this->knuthRandom();
|
156
|
|
|
} else {
|
157
|
|
|
return $this->methodPARandom();
|
158
|
|
|
}
|
159
|
|
|
|
160
|
|
|
}
|
161
|
|
|
|
162
|
|
|
/**
|
163
|
|
|
* WARNING: This function is of very limited use with Poisson distributions, and may represent a SIGNIFICANT
|
164
|
|
|
* performance hit for certain values of $min, $max, $lambda, and $maxIterations
|
165
|
|
|
*
|
166
|
|
|
* @param int|float|NumberInterface $min
|
167
|
|
|
* @param int|float|NumberInterface $max
|
168
|
|
|
* @param int $maxIterations
|
169
|
|
|
*
|
170
|
|
|
* @throws OptionalExit
|
171
|
|
|
* @return ImmutableNumber
|
172
|
|
|
*/
|
173
|
|
View Code Duplication |
public function rangeRandom($min = 0, $max = PHP_INT_MAX, $maxIterations = 20)
|
|
|
|
|
174
|
|
|
{
|
175
|
|
|
$i = 0;
|
176
|
|
|
|
177
|
|
|
do {
|
178
|
|
|
$randomNumber = $this->random();
|
179
|
|
|
$i++;
|
180
|
|
|
} while (($randomNumber->isGreaterThanOrEqualTo($max) || $randomNumber->isLessThanOrEqualTo($min)) && $i < $maxIterations);
|
181
|
|
|
|
182
|
|
|
if ($randomNumber->isGreaterThan($max) || $randomNumber->isLessThan($min)) {
|
183
|
|
|
throw new OptionalExit(
|
184
|
|
|
'All random numbers generated were outside of the requested range',
|
185
|
|
|
'A suitable random number, restricted by the $max ('.$max.') and $min ('.$min.'), could not be found within '.$maxIterations.' iterations'
|
186
|
|
|
);
|
187
|
|
|
} else {
|
188
|
|
|
return $randomNumber;
|
189
|
|
|
}
|
190
|
|
|
}
|
191
|
|
|
|
192
|
|
|
/**
|
193
|
|
|
* Method PA from The Computer Generation of Poisson Random Variables by A. C. Atkinson, 1979
|
194
|
|
|
* Journal of the Royal Statistical Society Series C, Vol. 28, No. 1, Pages 29-35
|
195
|
|
|
*
|
196
|
|
|
* As described by John D. Cook: http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/
|
197
|
|
|
*
|
198
|
|
|
* @return ImmutableNumber
|
199
|
|
|
*/
|
200
|
|
|
protected function methodPARandom()
|
201
|
|
|
{
|
202
|
|
|
$randFactory = new Factory();
|
203
|
|
|
$generator = $randFactory->getMediumStrengthGenerator();
|
204
|
|
|
/** @var ImmutableNumber $c */
|
205
|
|
|
$c = $this->lambda->pow(-1)->multiply(3.36)->multiply(-1)->add(0.767);
|
206
|
|
|
/** @var ImmutableNumber $beta */
|
207
|
|
|
$beta = Numbers::makePi()->divide($this->lambda->multiply(3)->sqrt());
|
208
|
|
|
/** @var ImmutableNumber $alpha */
|
209
|
|
|
$alpha = $this->lambda->multiply($beta);
|
210
|
|
|
/** @var ImmutableNumber $k */
|
211
|
|
|
$k = $c->ln(20)->subtract($this->lambda)->subtract($beta->ln(20));
|
212
|
|
|
/** @var ImmutableNumber $one */
|
213
|
|
|
$one = Numbers::makeOne();
|
214
|
|
|
/** @var ImmutableNumber $oneHalf */
|
215
|
|
|
$oneHalf = Numbers::make(Numbers::IMMUTABLE, '0.5');
|
216
|
|
|
/** @var ImmutableNumber $e */
|
217
|
|
|
$e = Numbers::makeE();
|
218
|
|
|
|
219
|
|
|
while (true) {
|
220
|
|
|
/** @var ImmutableNumber $u */
|
221
|
|
|
$u = $generator->generateInt() / PHP_INT_MAX;
|
222
|
|
|
/** @var ImmutableNumber $x */
|
223
|
|
|
$x = $alpha->subtract($one->subtract($u)->divide($u)->ln(20)->divide($beta));
|
|
|
|
|
224
|
|
|
/** @var ImmutableNumber $n */
|
225
|
|
|
$n = $x->add($oneHalf)->floor();
|
|
|
|
|
226
|
|
|
|
227
|
|
|
if ($n->isNegative()) {
|
228
|
|
|
continue;
|
229
|
|
|
}
|
230
|
|
|
|
231
|
|
|
/** @var ImmutableNumber $v */
|
232
|
|
|
$v = $generator->generateInt() / PHP_INT_MAX;
|
233
|
|
|
/** @var ImmutableNumber $y */
|
234
|
|
|
$y = $alpha->subtract($beta->multiply($x));
|
235
|
|
|
/** @var ImmutableNumber $lhs */
|
236
|
|
|
$lhs = $y->add($v->divide($e->pow($y)->add($one)->pow(2)));
|
237
|
|
|
/** @var ImmutableNumber $rhs */
|
238
|
|
|
$rhs = $k->add($n->multiply($this->lambda->ln(20)))->subtract($n->factorial()->ln(20));
|
|
|
|
|
239
|
|
|
|
240
|
|
|
if ($lhs->isLessThanOrEqualTo($rhs)) {
|
241
|
|
|
return $n;
|
242
|
|
|
}
|
243
|
|
|
|
244
|
|
|
/*
|
245
|
|
|
* At least attempt to free up some memory, since this particular method is extra hard on object instantiation
|
246
|
|
|
*/
|
247
|
|
|
unset($u);
|
248
|
|
|
unset($x);
|
249
|
|
|
unset($n);
|
250
|
|
|
unset($v);
|
251
|
|
|
unset($y);
|
252
|
|
|
unset($lhs);
|
253
|
|
|
unset($rhs);
|
254
|
|
|
}
|
255
|
|
|
}
|
256
|
|
|
|
257
|
|
|
/**
|
258
|
|
|
* @return ImmutableNumber
|
259
|
|
|
*/
|
260
|
|
|
protected function knuthRandom()
|
261
|
|
|
{
|
262
|
|
|
$randFactory = new Factory();
|
263
|
|
|
/** @var ImmutableNumber $L */
|
264
|
|
|
$L = Numbers::makeE()->pow($this->lambda->multiply(-1));
|
265
|
|
|
/** @var ImmutableNumber $k */
|
266
|
|
|
$k = Numbers::makeZero();
|
267
|
|
|
/** @var ImmutableNumber $p */
|
268
|
|
|
$p = Numbers::makeOne();
|
269
|
|
|
|
270
|
|
|
do {
|
271
|
|
|
$k = $k->add(1);
|
272
|
|
|
/** @var ImmutableNumber $u */
|
273
|
|
|
$u = $randFactory->getMediumStrengthGenerator()->generateInt() / PHP_INT_MAX;
|
274
|
|
|
$p = $p->multiply($u);
|
275
|
|
|
} while ($p->isGreaterThan($L));
|
276
|
|
|
|
277
|
|
|
return $k->subtract(1);
|
278
|
|
|
}
|
279
|
|
|
|
280
|
|
|
} |
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.
You can also find more detailed suggestions in the “Code” section of your repository.