import math, random

# Find mu,sigma for Normal Distr that approximate Bin(n,p)
#
#  mu = n*p
#  sigma = sqrt( n*p*(1-p) )

def normal_approximation_to_binomial(n, p):
    """finds mu and sigma corresponding to a Binomial(n, p)"""
    mu = p * n
    sigma = math.sqrt(p * (1 - p) * n)
    return mu, sigma

###################################################################
# normal_pdf(x, mu, sigma) = value of bell curve for x=x

def normal_pdf(x: float, mu: float = 0, sigma: float = 1) -> float:
    return (math.exp(-(x-mu) ** 2 / 2 / sigma ** 2) / (SQRT_TWO_PI * sigma))


#######################################################################
# normal_cdf(x, mu, sigma) = p = P[X <= x] for Normal(mu,sigma) var X

def normal_cdf(x: float, mu: float = 0, sigma: float = 1) -> float:
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2


#######################################################################
# inverse_normal_cdf(p, mu, sigma) = x for which
#
#     P[X <= x] = p for Normal(mu,sigma) var X
#
# I.e.: it's a probability LOOKUP function

def inverse_normal_cdf(p: float,
                       mu: float = 0,
                       sigma: float = 1,
                       tolerance: float = 0.00001) -> float:
    """Find approximate inverse using binary search"""

    # if not standard, compute standard and rescale
    if mu != 0 or sigma != 1:
        return mu + sigma * inverse_normal_cdf(p, tolerance=tolerance)

    low_z = -10.0                      # normal_cdf(-10) is (very close to) 0
    hi_z  =  10.0                      # normal_cdf(10)  is (very close to) 1
    while hi_z - low_z > tolerance:
        mid_z = (low_z + hi_z) / 2     # Consider the midpoint
        mid_p = normal_cdf(mid_z)      # and the cdf's value there
        if mid_p < p:
            low_z = mid_z              # Midpoint too low, search above it
        else:
            hi_z = mid_z               # Midpoint too high, search below it

    return mid_z



# ####################################################
# P[Z <= p]  ** This is just the CDF

normal_probability_below = normal_cdf

# ####################################################
# P[Z >= p] = 1 - P[Z <= p]

def normal_probability_above(lo, mu=0, sigma=1):
    return 1 - normal_cdf(lo, mu, sigma)

# ####################################################
# P[ lo <= Z <= hi] = P[Z <= hi] - P[Z <= lo]

def normal_probability_between(lo, hi, mu=0, sigma=1):
    return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)

# ####################################################
# P[ Z <= lo || Z >= hi] = 1 - (P[Z <= hi] - P[Z <= lo])

def normal_probability_outside(lo, hi, mu=0, sigma=1):
    return 1 - normal_probability_between(lo, hi, mu, sigma)



# ######################################################
# normal_upper_bound(p) = z where P[ Z <= z ] = p

def normal_upper_bound(probability, mu=0, sigma=1):
    """
    returns the z for which P(Z <= z) = probability
    """
    return inverse_normal_cdf(probability, mu, sigma)

# ######################################################
# normal_upper_bound(p) = z where P[ Z >= z ] = p

def normal_lower_bound(probability, mu=0, sigma=1):
    """
    returns the z for which P(Z >= z) = probability
    """
    return inverse_normal_cdf(1 - probability, mu, sigma)

# ######################################################
# normal_two_sided_bounds(p) = a,b where P[ a <= Z <= b ] = p
#
# a, b are equi-distance from mean (mu)

def normal_two_sided_bounds(probability, mu=0, sigma=1):
    """
    returns the symmetric (about the mean) bounds
    that contain the specified probability
    """

    tail_probability = (1 - probability) / 2
    # upper bound should have tail_probability above it
    upper_bound = normal_lower_bound(tail_probability, mu, sigma)

    # lower bound should have tail_probability below it
    lower_bound = normal_upper_bound(tail_probability, mu, sigma)
    return lower_bound, upper_bound