blob: 3c1e6feb52300854a6e5378d462f63585f7a8c57 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
|
import math
import random
import sys
if sys.version_info < (3, 2):
from repoze.lru import lru_cache
# from functools32 import lru_cache
else:
from functools import lru_cache
@lru_cache(1000)
def _get_poisson_var(parameter):
return Poisson(parameter)
class Poisson(object):
def __init__(self, parameter):
self.parameter = parameter
def c_p(k):
"""Compute the Poisson CDF for k iteratively."""
if k == 0:
return self._p(0)
else:
return self._compute_poisson_cdf(k - 1) + self._p(k)
self._compute_poisson_cdf = lru_cache(int(2.5*self.parameter) + 1)(c_p)
@staticmethod
def _get_random():
return random.random()
def _p(self, k):
# l^k * e^-l / k!
# Computed as exp(klog(l) - l - log(k!))
l = self.parameter
l_to_the_k = k * math.log(l)
k_fact = sum([math.log(i + 1) for i in range(k)])
return math.exp(l_to_the_k - l - k_fact)
def sample(self):
# The idea is:
# take a sample from U(0,1) and call it f.
# Let x be s.t. x = min_N F(x) > f
# where F is the cumulative distribution function of Poisson(parameter)
# return x
f = self._get_random()
# We compute x iteratively by computing
# \sum_k(P=k)
# where P ~ Poisson(parameter) and stopping as soon as
# it is greater than f.
# We use the cache to store results.
current_cdf = -1
current_x = -1
while current_cdf < f:
current_x += 1
current_cdf = self._compute_poisson_cdf(current_x)
return current_x
def poisson(parameter):
return _get_poisson_var(parameter).sample()
|
