LPC_stats_2017/HW2/p1.py

#!/usr/bin/env python3
""" LPC stats HW2, Problem 1
    Author: Caleb Fangmeier
    Created: Sep. 25, 2017
"""
from numpy import sum, sqrt
from numpy.random import exponential, poisson


def find_relative_frequency(n_experiments, b):
    # First get the mean counts for each experiment
    a = exponential(scale=b, size=n_experiments)

    # Since this is a counting experiment, sample from a Poisson distribution with the
    # previously generated means.
    x = poisson(lam=a)

    # Finally, count for how many experiments the a lies in
    # the range [x-sqrt(x), x+sqrt(x)]
    n_pass = sum((x - sqrt(x) < a) & (a < x + sqrt(x)))
    relative_frequency = n_pass / n_experiments

    # and print the results
    print(f"For b={b}, {n_pass}/{n_experiments} passed, R={relative_frequency*100:4.2f}%")


find_relative_frequency(10000, 5)
find_relative_frequency(10000, 10)

# Example output
# > For b=5, 6012/10000 passed, R=60.12%
# > For b=10, 6372/10000 passed, R=63.72%