###############################################################################
# Code needed to should "Confidence interval"
###############################################################################

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

from scratch.probability import normal_cdf, inverse_normal_cdf

# ####################################################
# 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





###############################################################################
# Confidence interval
###############################################################################

# We’ve been testing hypotheses about the value of the heads probability p, 
# which is a **parameter** of the unknown “heads” distribution. 
#
# When this is the case, a third approach is to construct 
# a confidence interval around the observed value of the parameter.

# For example, we can **estimate the probability** of the **unfair** coin 
# by looking at the average value of the Bernoulli variables corresponding 
# to each flip  1 if heads, 0 if tails. 

# E.g.: If we observe 525 heads out of 1,000 flips, 
# then we estimate p equals 0.525.

#################################################
# How confident can we be about this estimate?
#################################################

# IF we knew the exact value of p, the he Central Limit Theorem tells us 
# that the average of those Bernoulli variables should be approximately 
# normal, with:
#
#         mean                =   p 		and 
#         standard deviation: sqrt(n*p*(1-p))

# We don't know p, so instead, we use the ESTIMATE:
#
#     p = 0.525
#
# Our estimate is:

p_hat = 525 / 1000				# Our estimate for p

mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000) # 0.0158

print("mu = ", mu, ", sigma = ", sigma)

# Using the normal approximation, we conclude that we are “95% confident” 
# that the following interval contains the true parameter p:

conf_lo, conf_hi = normal_two_sided_bounds(0.95, mu, sigma)  # [0.4940, 0.5560]
print("Observed 525 heads out of 1000 tosses")
print("95% confidence interval = ", conf_lo, conf_hi)

# Meaning: if you were to repeat the experiment many times, 
#          95% of the time the “true” parameter p
#          would lie within the confidence interval

# NOTE: WE CANNOT conclude from 525 heads out of 1000 that
#       the coin is FAIR because:
#       95% confidence interval = [0.4940, 0.5560]
#       And:
# 		0.5 is INSIDE [0.4940, 0.5560] !!!

##################################################################
# How to use confidence interval
##################################################################

# Confidence interval technique can REJECT the statement that a coin is fair

# Suppose toss a coin 1000 and see 540 heads:

# Now: we assume p = 0.54

p_hat = 540 / 1000
mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000) # 0.0158

# Compute the 95% confidence interval for the ACTUAL p value:
conf_lo, conf_hi = normal_two_sided_bounds(0.95, mu, sigma) # [0.5091, 0.5709]

print()
print("Observed 525 heads out of 1000 tosses")
print("95% confidence interval = ", conf_lo, conf_hi)   # [0.5091, 0.5709]

# So: “fair coin” (0.5) doesn’t lie in the confidence interval.
# (The “fair coin” hypothesis doesn’t pass a test that you’d expect 
# it to pass 95% of the time if it were true.)