Tutorial 4: Classical Probability Tests & Null Models

"""Tutorial 4: Test if the lottery is behaving randomly using chi-squared tests"""

import pandas as pd
import numpy as np
from scipy.stats import chisquare

# --- Load feature table ---

print("Loading feature table from Tutorial 2...")
features = pd.read_parquet('data/processed/features_powerball.parquet')
print(f"Loaded {len(features)} draws\n")

# --- Test 1: Odd Count Distribution ---

print("Test 1: Odd Count Distribution")
print("-" * 40)

# Count how many draws have 0, 1, 2, 3, 4, 5 odd balls
observed_odd_counts = features['odd_count'].value_counts().sort_index()
observed_frequencies = [observed_odd_counts.get(i, 0) for i in range(6)]

print("\nObserved frequencies:")
for i, count in enumerate(observed_frequencies):
    print(f"  {i} odd balls: {count} draws")

# Calculate expected frequencies under the null hypothesis
# (binomial distribution with n=5, p=0.5)
total_draws = len(features)
expected_probabilities = [0.029, 0.153, 0.318, 0.330, 0.170, 0.034]
expected_frequencies = [prob * total_draws for prob in expected_probabilities]

print("\nExpected frequencies (if random):")
for i, count in enumerate(expected_frequencies):
    print(f"  {i} odd balls: {count:.1f} draws")

# Perform chi-squared test
chi2_stat, p_value = chisquare(observed_frequencies, expected_frequencies)

print(f"\nChi-squared statistic: {chi2_stat:.4f}")
print(f"P-value: {p_value:.4f}")
print(f"Degrees of freedom: {len(observed_frequencies) - 1}")

# Interpret the result
alpha = 0.05
if p_value < alpha:
    print(f"\nConclusion: REJECT null hypothesis (p < {alpha})")
    print("The odd count distribution does NOT match random expectation.")
else:
    print(f"\nConclusion: FAIL TO REJECT null hypothesis (p >= {alpha})")
    print("The odd count distribution is consistent with randomness.")

# --- Test 2: Decade Distribution ---

print("\n\nTest 2: Decade Distribution (Uniformity)")
print("-" * 40)

decade_columns = ['decade_00_09', 'decade_10_19', 'decade_20_29', 
                  'decade_30_39', 'decade_40_49', 'decade_50_59', 'decade_60_69']

# Count total hits in each decade
observed_decade_counts = [features[col].sum() for col in decade_columns]

print("\nObserved decade counts:")
decades = ['0-9', '10-19', '20-29', '30-39', '40-49', '50-59', '60-69']
for decade, count in zip(decades, observed_decade_counts):
    print(f"  {decade}: {count} hits")

# Under uniformity, each decade should get roughly the same number of hits
# But decade 0-9 only has 9 numbers vs 10 for others, so adjust expected values
decade_sizes = [9, 10, 10, 10, 10, 10, 10]  # Number of balls in each decade
total_size = sum(decade_sizes)
total_hits = sum(observed_decade_counts)

expected_decade_counts = [(size / total_size) * total_hits for size in decade_sizes]

print("\nExpected decade counts (if uniform):")
for decade, count in zip(decades, expected_decade_counts):
    print(f"  {decade}: {count:.1f} hits")

# Perform chi-squared test
chi2_stat, p_value = chisquare(observed_decade_counts, expected_decade_counts)

print(f"\nChi-squared statistic: {chi2_stat:.4f}")
print(f"P-value: {p_value:.4f}")
print(f"Degrees of freedom: {len(observed_decade_counts) - 1}")

# Interpret the result
if p_value < alpha:
    print(f"\nConclusion: REJECT null hypothesis (p < {alpha})")
    print("The decade distribution is NOT uniform.")
else:
    print(f"\nConclusion: FAIL TO REJECT null hypothesis (p >= {alpha})")
    print("The decade distribution is consistent with uniformity.")

# --- Summary ---

print("\n\nSummary")
print("-" * 40)
print("Both tests used alpha = 0.05 (5% significance level)")
print("If p-value < 0.05, we reject the null hypothesis")
print("If p-value >= 0.05, we fail to reject (data looks random)\n")