LPC_stats_2017/HW2/p2.py

#!/usr/bin/env python3
""" LPC stats HW2, Problem 2
    Author: Caleb Fangmeier
    Created: Sep. 27, 2017
"""
from scipy.stats import beta, poisson
from scipy.optimize import root

# Declare the experimental results
N = 25
M = 353
k = 37.6
x = 1/(1+k)

p = 0.683


def DL(z):
    """
        This is just the DL function as defined in the problem statement. It is the
        cumulative distribution function for p(D|s).
    """
    num = 0
    den = 0
    for r in range(N+1):
        b = beta.pdf(x, r+1, M)
        num += b*poisson.cdf(N-r, z)
        den += b
    return num/den


def DR(z):
    """
       Same as DL, but for the opposite side.
    """
    return 1 - DL(z)


# Scipy's root function finds a zero crossing for the given function.
# The second argument is just the starting point for the algorithm.
root_DL = root(lambda z: DL(z) - (1-p)/2, 10).x[0]
root_DR = root(lambda z: DR(z) - (1-p)/2, 10).x[0]
print(f"The interval is: [{root_DR:.2f}, {root_DL:.2f}]")

# Output
# > The interval is: [11.49, 21.68]