|
1
|
|
|
<?php
|
|
2
|
|
|
|
|
3
|
|
|
namespace Samsara\Fermat\Provider\Distribution;
|
|
4
|
|
|
|
|
5
|
|
|
use RandomLib\Factory;
|
|
6
|
|
|
use Samsara\Exceptions\SystemError\LogicalError\IncompatibleObjectState;
|
|
7
|
|
|
use Samsara\Exceptions\UsageError\IntegrityConstraint;
|
|
8
|
|
|
use Samsara\Exceptions\UsageError\OptionalExit;
|
|
9
|
|
|
use Samsara\Fermat\Numbers;
|
|
10
|
|
|
use Samsara\Fermat\Provider\Distribution\Base\Distribution;
|
|
11
|
|
|
use Samsara\Fermat\Types\Base\DecimalInterface;
|
|
12
|
|
|
use Samsara\Fermat\Types\Base\NumberInterface;
|
|
13
|
|
|
use Samsara\Fermat\Values\ImmutableNumber;
|
|
14
|
|
|
|
|
15
|
|
|
class Poisson extends Distribution
|
|
16
|
|
|
{
|
|
17
|
|
|
|
|
18
|
|
|
/**
|
|
19
|
|
|
* @var ImmutableNumber
|
|
20
|
|
|
*/
|
|
21
|
|
|
private $lambda;
|
|
22
|
|
|
|
|
23
|
|
|
/**
|
|
24
|
|
|
* Poisson constructor.
|
|
25
|
|
|
*
|
|
26
|
|
|
* @param int|float|DecimalInterface $lambda
|
|
27
|
|
|
*
|
|
28
|
|
|
* @throws IntegrityConstraint
|
|
29
|
|
|
*/
|
|
30
|
|
|
public function __construct($lambda)
|
|
31
|
|
|
{
|
|
32
|
|
|
$lambda = Numbers::makeOrDont(Numbers::IMMUTABLE, $lambda);
|
|
33
|
|
|
|
|
34
|
|
|
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
|
|
|
$this->lambda = $lambda;
|
|
|
|
|
|
|
43
|
|
|
}
|
|
44
|
|
|
|
|
45
|
|
|
/**
|
|
46
|
|
|
* @param int|float|DecimalInterface $k
|
|
47
|
|
|
*
|
|
48
|
|
|
* @return ImmutableNumber
|
|
49
|
|
|
* @throws IntegrityConstraint
|
|
50
|
|
|
* @throws IncompatibleObjectState
|
|
51
|
|
|
*/
|
|
52
|
|
|
public function probabilityOfKEvents($k): ImmutableNumber
|
|
53
|
|
|
{
|
|
54
|
|
|
|
|
55
|
|
|
return $this->pmf($k);
|
|
56
|
|
|
|
|
57
|
|
|
}
|
|
58
|
|
|
|
|
59
|
|
|
/**
|
|
60
|
|
|
* @param int|float|DecimalInterface $x
|
|
61
|
|
|
*
|
|
62
|
|
|
* @return ImmutableNumber
|
|
63
|
|
|
* @throws IntegrityConstraint
|
|
64
|
|
|
* @throws IncompatibleObjectState
|
|
65
|
|
|
*/
|
|
66
|
|
|
public function cdf($x): ImmutableNumber
|
|
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
|
|
|
$cumulative = Numbers::makeZero();
|
|
80
|
|
|
|
|
81
|
|
|
for ($i = 0;$x->isGreaterThanOrEqualTo($i);$i++) {
|
|
82
|
|
|
$cumulative = $cumulative->add($this->pmf($i));
|
|
83
|
|
|
}
|
|
84
|
|
|
|
|
85
|
|
|
return $cumulative;
|
|
86
|
|
|
|
|
87
|
|
|
}
|
|
88
|
|
|
|
|
89
|
|
|
/**
|
|
90
|
|
|
* @param int|float|DecimalInterface $x
|
|
91
|
|
|
*
|
|
92
|
|
|
* @return ImmutableNumber
|
|
93
|
|
|
* @throws IntegrityConstraint
|
|
94
|
|
|
* @throws IncompatibleObjectState
|
|
95
|
|
|
*/
|
|
96
|
|
|
public function pmf($x): ImmutableNumber
|
|
97
|
|
|
{
|
|
98
|
|
|
$x = Numbers::makeOrDont(Numbers::IMMUTABLE, $x);
|
|
99
|
|
|
|
|
100
|
|
|
if (!$x->isNatural()) {
|
|
101
|
|
|
throw new IntegrityConstraint(
|
|
102
|
|
|
'Only integers are valid x values for Poisson distributions',
|
|
103
|
|
|
'Provide an integer value to calculate the PMF',
|
|
104
|
|
|
'Poisson distributions describe discrete occurrences; only integers are valid x values'
|
|
105
|
|
|
);
|
|
106
|
|
|
}
|
|
107
|
|
|
|
|
108
|
|
|
$e = Numbers::makeE();
|
|
109
|
|
|
|
|
110
|
|
|
/** @var ImmutableNumber $pmf */
|
|
111
|
|
|
$pmf = $this->lambda->pow($x)->multiply($e->pow($this->lambda->multiply(-1)))->divide($x->factorial());
|
|
|
|
|
|
|
112
|
|
|
|
|
113
|
|
|
return $pmf;
|
|
114
|
|
|
}
|
|
115
|
|
|
|
|
116
|
|
|
/**
|
|
117
|
|
|
* @param int|float|DecimalInterface $x1
|
|
118
|
|
|
* @param int|float|DecimalInterface $x2
|
|
119
|
|
|
*
|
|
120
|
|
|
* @return ImmutableNumber
|
|
121
|
|
|
* @throws IntegrityConstraint
|
|
122
|
|
|
*/
|
|
123
|
|
|
public function rangePmf($x1, $x2): ImmutableNumber
|
|
124
|
|
|
{
|
|
125
|
|
|
$x1 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x1);
|
|
126
|
|
|
$x2 = Numbers::makeOrDont(Numbers::IMMUTABLE, $x2);
|
|
127
|
|
|
|
|
128
|
|
|
if ($x1->equals($x2)) {
|
|
|
|
|
|
|
129
|
|
|
return Numbers::makeZero();
|
|
130
|
|
|
} elseif ($x1->isGreaterThan($x2)) {
|
|
131
|
|
|
$larger = $x1;
|
|
132
|
|
|
$smaller = $x2;
|
|
133
|
|
|
} else {
|
|
134
|
|
|
$larger = $x2;
|
|
135
|
|
|
$smaller = $x1;
|
|
136
|
|
|
}
|
|
137
|
|
|
|
|
138
|
|
|
if (!$larger->isNatural() || !$smaller->isNatural()) {
|
|
139
|
|
|
throw new IntegrityConstraint(
|
|
140
|
|
|
'Only integers are valid x values for Poisson distributions',
|
|
141
|
|
|
'Provide integer values to calculate the range PMF',
|
|
142
|
|
|
'Poisson distributions describe discrete occurrences; only integers are valid x values'
|
|
143
|
|
|
);
|
|
144
|
|
|
}
|
|
145
|
|
|
|
|
146
|
|
|
$cumulative = Numbers::makeZero();
|
|
147
|
|
|
|
|
148
|
|
|
for (;$larger->isGreaterThanOrEqualTo($smaller);$smaller->add(1)) {
|
|
149
|
|
|
$cumulative = $cumulative->add($this->pmf($smaller));
|
|
|
|
|
|
|
150
|
|
|
}
|
|
151
|
|
|
|
|
152
|
|
|
return $cumulative;
|
|
153
|
|
|
}
|
|
154
|
|
|
|
|
155
|
|
|
/**
|
|
156
|
|
|
* @return ImmutableNumber
|
|
157
|
|
|
* @throws IntegrityConstraint
|
|
158
|
|
|
* @throws IncompatibleObjectState
|
|
159
|
|
|
*/
|
|
160
|
|
|
public function random(): ImmutableNumber
|
|
161
|
|
|
{
|
|
162
|
|
|
if ($this->lambda->isLessThanOrEqualTo(30)) {
|
|
163
|
|
|
return $this->knuthRandom();
|
|
164
|
|
|
} else {
|
|
165
|
|
|
return $this->methodPARandom();
|
|
166
|
|
|
}
|
|
167
|
|
|
}
|
|
168
|
|
|
|
|
169
|
|
|
/**
|
|
170
|
|
|
* WARNING: This function is of very limited use with Poisson distributions, and may represent a SIGNIFICANT
|
|
171
|
|
|
* performance hit for certain values of $min, $max, $lambda, and $maxIterations
|
|
172
|
|
|
*
|
|
173
|
|
|
* @param int|float|NumberInterface $min
|
|
174
|
|
|
* @param int|float|NumberInterface $max
|
|
175
|
|
|
* @param int $maxIterations
|
|
176
|
|
|
*
|
|
177
|
|
|
* @throws OptionalExit
|
|
178
|
|
|
* @return ImmutableNumber
|
|
179
|
|
|
*/
|
|
180
|
|
|
public function rangeRandom($min = 0, $max = PHP_INT_MAX, int $maxIterations = 20): ImmutableNumber
|
|
181
|
|
|
{
|
|
182
|
|
|
$i = 0;
|
|
183
|
|
|
|
|
184
|
|
|
do {
|
|
185
|
|
|
$randomNumber = $this->random();
|
|
186
|
|
|
$i++;
|
|
187
|
|
|
} while (($randomNumber->isGreaterThan($max) || $randomNumber->isLessThan($min)) && $i < $maxIterations);
|
|
188
|
|
|
|
|
189
|
|
|
if ($randomNumber->isGreaterThan($max) || $randomNumber->isLessThan($min)) {
|
|
190
|
|
|
throw new OptionalExit(
|
|
191
|
|
|
'All random numbers generated were outside of the requested range',
|
|
192
|
|
|
'A suitable random number, restricted by the $max ('.$max.') and $min ('.$min.'), could not be found within '.$maxIterations.' iterations'
|
|
|
|
|
|
|
193
|
|
|
);
|
|
194
|
|
|
} else {
|
|
195
|
|
|
return $randomNumber;
|
|
196
|
|
|
}
|
|
197
|
|
|
}
|
|
198
|
|
|
|
|
199
|
|
|
/**
|
|
200
|
|
|
* Method PA from The Computer Generation of Poisson Random Variables by A. C. Atkinson, 1979
|
|
201
|
|
|
* Journal of the Royal Statistical Society Series C, Vol. 28, No. 1, Pages 29-35
|
|
202
|
|
|
*
|
|
203
|
|
|
* As described by John D. Cook: http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/
|
|
204
|
|
|
*
|
|
205
|
|
|
* @return ImmutableNumber
|
|
206
|
|
|
* @throws IntegrityConstraint
|
|
207
|
|
|
* @throws IncompatibleObjectState
|
|
208
|
|
|
*/
|
|
209
|
|
|
protected function methodPARandom(): ImmutableNumber
|
|
210
|
|
|
{
|
|
211
|
|
|
$randFactory = new Factory();
|
|
212
|
|
|
$generator = $randFactory->getMediumStrengthGenerator();
|
|
213
|
|
|
/** @var ImmutableNumber $c */
|
|
214
|
|
|
$c = $this->lambda->pow(-1)->multiply(3.36)->multiply(-1)->add(0.767);
|
|
|
|
|
|
|
215
|
|
|
/** @var ImmutableNumber $beta */
|
|
216
|
|
|
$beta = Numbers::makePi()->divide($this->lambda->multiply(3)->sqrt());
|
|
|
|
|
|
|
217
|
|
|
/** @var ImmutableNumber $alpha */
|
|
218
|
|
|
$alpha = $this->lambda->multiply($beta);
|
|
219
|
|
|
/** @var ImmutableNumber $k */
|
|
220
|
|
|
$k = $c->ln(20)->subtract($this->lambda)->subtract($beta->ln(20));
|
|
|
|
|
|
|
221
|
|
|
/** @var ImmutableNumber $one */
|
|
222
|
|
|
$one = Numbers::makeOne();
|
|
223
|
|
|
/** @var ImmutableNumber $oneHalf */
|
|
224
|
|
|
$oneHalf = Numbers::make(Numbers::IMMUTABLE, '0.5');
|
|
225
|
|
|
/** @var ImmutableNumber $e */
|
|
226
|
|
|
$e = Numbers::makeE();
|
|
227
|
|
|
|
|
228
|
|
|
while (true) {
|
|
229
|
|
|
/** @var ImmutableNumber $u */
|
|
230
|
|
|
$u = $generator->generateInt() / PHP_INT_MAX;
|
|
231
|
|
|
/** @var ImmutableNumber $x */
|
|
232
|
|
|
$x = $alpha->subtract($one->subtract($u)->divide($u)->ln(20)->divide($beta));
|
|
|
|
|
|
|
233
|
|
|
/** @var ImmutableNumber $n */
|
|
234
|
|
|
$n = $x->add($oneHalf)->floor();
|
|
235
|
|
|
|
|
236
|
|
|
if ($n->isNegative()) {
|
|
237
|
|
|
continue;
|
|
238
|
|
|
}
|
|
239
|
|
|
|
|
240
|
|
|
/** @var ImmutableNumber $v */
|
|
241
|
|
|
$v = $generator->generateInt() / PHP_INT_MAX;
|
|
242
|
|
|
/** @var ImmutableNumber $y */
|
|
243
|
|
|
$y = $alpha->subtract($beta->multiply($x));
|
|
244
|
|
|
/** @var ImmutableNumber $lhs */
|
|
245
|
|
|
$lhs = $y->add($v->divide($e->pow($y)->add($one)->pow(2)));
|
|
|
|
|
|
|
246
|
|
|
/** @var ImmutableNumber $rhs */
|
|
247
|
|
|
$rhs = $k->add($n->multiply($this->lambda->ln(20)))->subtract($n->factorial()->ln(20));
|
|
248
|
|
|
|
|
249
|
|
|
if ($lhs->isLessThanOrEqualTo($rhs)) {
|
|
250
|
|
|
return $n;
|
|
251
|
|
|
}
|
|
252
|
|
|
|
|
253
|
|
|
/*
|
|
254
|
|
|
* At least attempt to free up some memory, since this particular method is extra hard on object instantiation
|
|
255
|
|
|
*/
|
|
256
|
|
|
unset($u);
|
|
257
|
|
|
unset($x);
|
|
258
|
|
|
unset($n);
|
|
259
|
|
|
unset($v);
|
|
260
|
|
|
unset($y);
|
|
261
|
|
|
unset($lhs);
|
|
262
|
|
|
unset($rhs);
|
|
263
|
|
|
}
|
|
|
|
|
|
|
264
|
|
|
}
|
|
265
|
|
|
|
|
266
|
|
|
/**
|
|
267
|
|
|
* @return ImmutableNumber
|
|
268
|
|
|
* @throws IntegrityConstraint
|
|
269
|
|
|
*/
|
|
270
|
|
|
protected function knuthRandom(): ImmutableNumber
|
|
271
|
|
|
{
|
|
272
|
|
|
$randFactory = new Factory();
|
|
273
|
|
|
/** @var ImmutableNumber $L */
|
|
274
|
|
|
$L = Numbers::makeE()->pow($this->lambda->multiply(-1));
|
|
275
|
|
|
/** @var ImmutableNumber $k */
|
|
276
|
|
|
$k = Numbers::makeZero();
|
|
277
|
|
|
/** @var ImmutableNumber $p */
|
|
278
|
|
|
$p = Numbers::makeOne();
|
|
279
|
|
|
|
|
280
|
|
|
do {
|
|
281
|
|
|
$k = $k->add(1);
|
|
282
|
|
|
/** @var ImmutableNumber $u */
|
|
283
|
|
|
$u = $randFactory->getMediumStrengthGenerator()->generateInt() / PHP_INT_MAX;
|
|
284
|
|
|
$p = $p->multiply($u);
|
|
285
|
|
|
} while ($p->isGreaterThan($L));
|
|
286
|
|
|
|
|
287
|
|
|
return $k->subtract(1);
|
|
288
|
|
|
}
|
|
289
|
|
|
|
|
290
|
|
|
} |
Our type inference engine has found a suspicous assignment of a value to a property. This check raises an issue when a value that can be of a mixed type is assigned to a property that is type hinted more strictly.
For example, imagine you have a variable
$accountIdthat can either hold an Id object or false (if there is no account id yet). Your code now assigns that value to theidproperty of an instance of theAccountclass. This class holds a proper account, so the id value must no longer be false.Either this assignment is in error or a type check should be added for that assignment.