Question

I have a simple script to set up a Poisson distribution by constructing an array of "events" of probability = 0.1, and then counting the number of successes in each group of 10. It almost works, but the distribution is not quite right (P(0) should equal P(1), but is instead about 90% of P(1)). It's like there's an off-by-one kind of error, but I can't figure out what it is. The script uses the Counter class from here (because I have Python 2.6 and not 2.7) and the grouping uses itertools as discussed here. It's not a stochastic issue, repeats give pretty tight results, and the overall mean looks good, group size looks good. Any ideas where I've messed up?

from itertools import izip_longest
import numpy as np
import Counter

def groups(iterable, n=3, padvalue=0):
    "groups('abcde', 3, 'x') --> ('a','b','c'), ('d','e','x')"
    return izip_longest(*[iter(iterable)]*n, fillvalue=padvalue)

def event():
    f = 0.1
    r = np.random.random()
    if r < f:  return 1
    return 0

L = [event() for i in range(100000)]
rL = [sum(g) for g in groups(L,n=10)]
print len(rL)
print sum(list(L))

C = Counter.Counter(rL)
for i in range(max(C.keys())+1):
    print str(i).rjust(2), C[i]

$ python script.py 
10000
9949
 0 3509
 1 3845
 2 1971
 3 555
 4 104
 5 15
 6 1
$ python script.py 
10000
10152
 0 3417
 1 3879
 2 1978
 3 599
 4 115
 5 12
Was it helpful?

Solution

I did a combinatorial reality check on your math, and it looks like your results are correct actually. P(0) should not be roughly equivalent to P(1)

.9^10 = 0.34867844 = probability of 0 events
.1 * .9^9 * (10 choose 1) = .1 * .9^9 * 10 = 0.387420489 = probability of 1 event

I wonder if you accidentally did your math thusly:

.1 * .9^10 * (10 choose 1) = 0.34867844 = incorrect probability of 1 event
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top