poisson distribution on sagemath

poisson_sagemath

# this code can be run on cloud.sagemath.org


per_day = 1000000
per_sec = per_day/(24 * 3600)
cache = {}

def poisson(mean_per_second, x):
    return e ** (-mean_per_second) * mean_per_second ** x / factorial(x)
  
def poisson_cached(mean_per_second, x):
    if mean_per_second in cache:
        if x in cache[mean_per_second]:
            res = cache[mean_per_second][x]
        else:
            res = poisson(mean_per_second, x)
            cache[mean_per_second][x] = res
    else:
        res = poisson(mean_per_second, x)
        cache[mean_per_second] = {x: res}
    return res

def prob_range(mean_per_second, start=0, end=0):
    return sum(poisson_cached(mean_per_second, x) for x in range(start, end))
  
f = lambda v: float(1-prob_range(110, end=v))
line_poins = [[xv, f(xv)] for xv in  range(150, 190)]

plot(line(res), axes_labels=("more hits  than", "probability"), frame=True, axes=True)

poisson_using_decimal.py

# pure python implementation

from math import e, factorial
from decimal import Decimal


per_day = 1000000
per_sec = per_day/(24 * 3600)
cache = {}

def poisson(mean_per_second, x):
    return Decimal(e) ** Decimal(-mean_per_second) * mean_per_second ** x / factorial(x)
  
def poisson_cached(mean_per_second, x):
    if mean_per_second in cache:
        if x in cache[mean_per_second]:
            res = cache[mean_per_second][x]
        else:
            res = poisson(mean_per_second, x)
            cache[mean_per_second][x] = res
    else:
        res = poisson(mean_per_second, x)
        cache[mean_per_second] = {x: res}
    return res

def prob_range(mean_per_second, start=0, end=0):
    return sum(poisson_cached(mean_per_second, x) for x in range(start, end))
  
f = lambda v: float(1-prob_range(110, end=v))
line_poins = [[xv, f(xv)] for xv in  range(150, 190)]