# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Instance for tossing a coin

# +
import math
import random


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

2
def inverse_normal_cdf(p, mu = 0, sigma = 1, tolerance = 0.00001):
    # 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


# -

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


# The normal cdf _is_ the probability the variable is below a threshold
normal_probability_below = normal_cdf


# It's above the threshold if it's not below the threshold
def normal_probability_above(lo, mu =0, sigma = 1):
    # the probability that a N(mu, sigma) is greater than lo
    return 1 - normal_cdf(lo, mu, sigma)


# It's between if it's less than hi, but not less than lo
def normal_probability_between(lo, hi, mu = 0, sigma = 1):
    # the probability that a N(mu, sigma) is between lo and hi
    return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)


# It's outside if it's not between
def normal_probability_outside(lo, hi, mu = 0, sigma = 1):
    # the probability that a N(mu, sigma) is not between lo and hi
    return 1 - normal_probability_between(lo, hi, mu, sigma)


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)


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)


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


# +
mu_0, sigma_0 = normal_approximation_to_binomial(1000, 0.5)

lower_bound, upper_bound = normal_two_sided_bounds(0.95, mu_0, sigma_0)  # (469, 531)

# 95% bounds based on assumption p is 0.5
lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)

# actual mu and sigma based on p = 0.55
mu_1, sigma_1 = normal_approximation_to_binomial(1000, 0.55)

# a type 2 error means we fail to reject the null hypothesis
# which will happen when X is still in our original interval
type_2_probability = normal_probability_between(lo, hi, mu_1, sigma_1)
power = 1 - type_2_probability      # 0.887

# +
hi = normal_upper_bound(0.95, mu_0, sigma_0)
# is 526 (< 531, since we need more probability in the upper tail)

type_2_probability = normal_probability_below(hi, mu_1, sigma_1)

power = 1 - type_2_probability      # 0.936


# -

# # P-values

# +
def two_sided_p_value(x, mu = 0, sigma = 1):
    # how likely are we to see a value at least as extreme as x (in either
    # direction) if our values are from a N(mu, sigma)?
    if x >= mu:
        # x is greater than the mean, so the tail is everything greater than x
        return 2 * normal_probability_above(x, mu, sigma)
    else:
        # x is less than the mean, so the tail is everything less than x
        return 2 * normal_probability_below(x, mu, sigma)

two_sided_p_value(529.5, mu_0, sigma_0)   # 0.062

# +
extreme_value_count = 0
for _ in range(1000):
    num_heads = sum(1 if random.random() < 0.5 else 0    # Count # of heads
                    for _ in range(1000))                # in 1000 flips,
    if num_heads >= 530 or num_heads <= 470:             # and count how often
        extreme_value_count += 1                         # the # is 'extreme'

print(extreme_value_count / 100000)
# -

two_sided_p_value(531.5, mu_0, sigma_0)   # 0.0463

# +
upper_p_value = normal_probability_above
lower_p_value = normal_probability_below

upper_p_value(524.5, mu_0, sigma_0) # 0.061

upper_p_value(526.5, mu_0, sigma_0) # 0.047
# -

# # Confidence interval 

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

normal_two_sided_bounds(0.95, mu, sigma)        # [0.4940, 0.5560]

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

normal_two_sided_bounds(0.95, mu, sigma) # [0.5091, 0.5709]


# -

# # P-value hacking

def run_experiment():
    # flips a fair coin 1000 times, True = heads, False = tails
    return [random.random() < 0.5 for _ in range(1000)]


def reject_fairness(experiment):
    # using the 5% significance levels
    num_heads = len([flip for flip in experiment if flip])
    return num_heads < 469 or num_heads > 531


# +
random.seed(0)
experiments = [run_experiment() for _ in range(1000)]
num_rejections = len([experiment
                      for experiment in experiments
                      if reject_fairness(experiment)])


num_rejections  # 46


# -

# # A/B tests

def estimated_parameters(N, n):
    p = n / N
    sigma = math.sqrt(p * (1 - p) / N)
    return p, sigma


# +
def a_b_test_statistic(N_A, n_A, N_B, n_B):
    p_A, sigma_A = estimated_parameters(N_A, n_A)
    p_B, sigma_B = estimated_parameters(N_B, n_B)
    return (p_B - p_A) / math.sqrt(sigma_A ** 2 + sigma_B ** 2)

z = a_b_test_statistic(1000, 200, 1000, 180)    # -1.14
z
# -

two_sided_p_value(z)

z = a_b_test_statistic(1000, 200, 1000, 150)    # -2.94
two_sided_p_value(z)                            # 0.003


# # Bayesian statistic inference

def B(alpha, beta):
    # a normalizing constant so that the total probability is 1
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)


def beta_pdf(x, alpha, beta):
    if x <= 0 or x >= 1:          # no weight outside of [0, 1]
        return 0
    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)



StarcoderdataJupyterScriptsDedupFiltered