# Our A/B test hypotheses:
# H₀: p_green = p_blue (no difference)
# H₁: p_green ≠ p_blue (there is a difference)

import numpy as np
from scipy import stats

# Data
blue_visitors = 10000
blue_purchases = 320
blue_rate = blue_purchases / blue_visitors

green_visitors = 10000
green_purchases = 355
green_rate = green_purchases / green_visitors

print(f"Blue conversion rate: {blue_rate:.2%}")
print(f"Green conversion rate: {green_rate:.2%}")
print(f"Observed difference: {green_rate - blue_rate:.2%}")

def two_proportion_z_test(x1, n1, x2, n2):
    """
    Test if two proportions are significantly different.
    
    x1, x2: number of successes
    n1, n2: number of trials
    
    Returns: z-statistic, p-value (two-tailed)
    """
    # Sample proportions
    p1 = x1 / n1
    p2 = x2 / n2
    
    # Pooled proportion (under null hypothesis)
    p_pool = (x1 + x2) / (n1 + n2)
    
    # Standard error under null
    se = np.sqrt(p_pool * (1 - p_pool) * (1/n1 + 1/n2))
    
    # Z-statistic
    z = (p1 - p2) / se
    
    # Two-tailed p-value
    p_value = 2 * (1 - stats.norm.cdf(abs(z)))
    
    return z, p_value

# Run the test
z_stat, p_value = two_proportion_z_test(
    green_purchases, green_visitors,
    blue_purchases, blue_visitors
)

print(f"\nZ-statistic: {z_stat:.3f}")
print(f"P-value: {p_value:.4f}")

if p_value < 0.05:
    print("\nResult: Reject null hypothesis")
    print("The difference is statistically significant at α=0.05")
else:
    print("\nResult: Fail to reject null hypothesis")
    print("The difference is NOT statistically significant at α=0.05")

Blue conversion rate: 3.20%
Green conversion rate: 3.55%
Observed difference: 0.35%

Z-statistic: 1.404
P-value: 0.1603

Result: Fail to reject null hypothesis
The difference is NOT statistically significant at α=0.05

from scipy.stats import chi2_contingency

# Contingency table
#                  Purchased    Not Purchased
# Blue               320          9680
# Green              355          9645

contingency_table = np.array([
    [blue_purchases, blue_visitors - blue_purchases],
    [green_purchases, green_visitors - green_purchases]
])

chi2, p_value, dof, expected = chi2_contingency(contingency_table)

print(f"Chi-square statistic: {chi2:.3f}")
print(f"P-value: {p_value:.4f}")
print(f"Degrees of freedom: {dof}")

Chi-square statistic: 1.972
P-value: 0.1603
Degrees of freedom: 1

# Visualize the tradeoff
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Type I Error
x = np.linspace(-4, 4, 1000)
null_dist = stats.norm.pdf(x, 0, 1)

axes[0].plot(x, null_dist, 'b-', linewidth=2, label='Null Distribution')
axes[0].fill_between(x[x > 1.96], null_dist[x > 1.96], alpha=0.3, color='red', 
                      label=f'Type I Error Region (α/2)')
axes[0].fill_between(x[x < -1.96], null_dist[x < -1.96], alpha=0.3, color='red')
axes[0].axvline(1.96, color='red', linestyle='--')
axes[0].axvline(-1.96, color='red', linestyle='--')
axes[0].set_title('Type I Error (False Positive)')
axes[0].legend()

# Type II Error (with alternative distribution)
alt_dist = stats.norm.pdf(x, 2, 1)  # Effect exists, shifted right

axes[1].plot(x, null_dist, 'b-', linewidth=2, label='Null (No Effect)')
axes[1].plot(x, alt_dist, 'g-', linewidth=2, label='Alternative (Real Effect)')
axes[1].fill_between(x[(x > -1.96) & (x < 1.96)], alt_dist[(x > -1.96) & (x < 1.96)], 
                      alpha=0.3, color='orange', label='Type II Error Region (β)')
axes[1].axvline(1.96, color='red', linestyle='--')
axes[1].axvline(-1.96, color='red', linestyle='--')
axes[1].set_title('Type II Error (False Negative)')
axes[1].legend()

plt.tight_layout()
plt.show()

def power_proportion_test(p1, p2, n, alpha=0.05):
    """
    Calculate power for a two-proportion test.
    
    p1: control proportion
    p2: treatment proportion
    n: sample size per group
    alpha: significance level
    """
    # Effect size
    effect = abs(p2 - p1)
    
    # Pooled standard error under null
    p_pool = (p1 + p2) / 2
    se_null = np.sqrt(2 * p_pool * (1 - p_pool) / n)
    
    # Standard error under alternative
    se_alt = np.sqrt((p1 * (1 - p1) + p2 * (1 - p2)) / n)
    
    # Critical value
    z_crit = stats.norm.ppf(1 - alpha / 2)
    
    # Power
    z_power = (effect - z_crit * se_null) / se_alt
    power = stats.norm.cdf(z_power)
    
    return power

# Our A/B test: 3.20% vs 3.55%, n=10,000 per group
power = power_proportion_test(0.032, 0.0355, 10000)
print(f"Power of our test: {power:.1%}")

# What if we had 50,000 per group?
power_large = power_proportion_test(0.032, 0.0355, 50000)
print(f"Power with n=50,000: {power_large:.1%}")

Power of our test: 27.3%
Power with n=50,000: 75.8%

def sample_size_proportion_test(p1, p2, power=0.80, alpha=0.05):
    """
    Calculate required sample size per group.
    
    p1: expected control proportion
    p2: expected treatment proportion
    power: desired power (typically 0.80)
    alpha: significance level (typically 0.05)
    """
    # Effect size
    effect = abs(p2 - p1)
    
    # Pooled proportion
    p_pool = (p1 + p2) / 2
    
    # Z-scores
    z_alpha = stats.norm.ppf(1 - alpha / 2)
    z_beta = stats.norm.ppf(power)
    
    # Variance terms
    var_null = 2 * p_pool * (1 - p_pool)
    var_alt = p1 * (1 - p1) + p2 * (1 - p2)
    
    # Sample size formula
    n = ((z_alpha * np.sqrt(var_null) + z_beta * np.sqrt(var_alt)) / effect) ** 2
    
    return int(np.ceil(n))

# How many visitors do we need to detect a 0.35% difference?
n_needed = sample_size_proportion_test(0.032, 0.0355)
print(f"Need {n_needed:,} visitors per group to detect 0.35% difference with 80% power")

# For a larger effect (1% improvement)
n_1pct = sample_size_proportion_test(0.032, 0.042)
print(f"Need {n_1pct:,} visitors per group to detect 1.0% difference with 80% power")

Need 48,614 visitors per group to detect 0.35% difference with 80% power
Need 6,038 visitors per group to detect 1.0% difference with 80% power

# Are our website load times different from the 3-second industry standard?
load_times = np.array([2.8, 3.2, 2.9, 3.5, 2.7, 3.1, 2.6, 3.0, 2.9, 3.3])

t_stat, p_value = stats.ttest_1samp(load_times, 3.0)
print(f"Sample mean: {np.mean(load_times):.2f}s")
print(f"t-statistic: {t_stat:.3f}")
print(f"p-value: {p_value:.4f}")

# Do users spend more time on new homepage design?
old_design_time = np.array([45, 52, 38, 61, 42, 55, 48, 50, 44, 58])
new_design_time = np.array([58, 62, 55, 70, 65, 60, 68, 72, 63, 59])

t_stat, p_value = stats.ttest_ind(old_design_time, new_design_time)
print(f"Old design mean: {np.mean(old_design_time):.1f}s")
print(f"New design mean: {np.mean(new_design_time):.1f}s")
print(f"t-statistic: {t_stat:.3f}")
print(f"p-value: {p_value:.4f}")

# Does a training program improve test scores?
before = np.array([65, 72, 58, 80, 75, 62, 70, 68, 74, 78])
after = np.array([70, 78, 62, 85, 82, 68, 75, 72, 80, 82])

t_stat, p_value = stats.ttest_rel(before, after)
print(f"Mean improvement: {np.mean(after - before):.1f} points")
print(f"t-statistic: {t_stat:.3f}")
print(f"p-value: {p_value:.4f}")

# Do three different ad campaigns have different click rates?
campaign_a = np.array([2.1, 2.3, 2.0, 2.4, 2.2])
campaign_b = np.array([2.8, 3.0, 2.9, 3.1, 2.7])
campaign_c = np.array([2.3, 2.5, 2.4, 2.6, 2.2])

f_stat, p_value = stats.f_oneway(campaign_a, campaign_b, campaign_c)
print(f"F-statistic: {f_stat:.3f}")
print(f"p-value: {p_value:.4f}")

import numpy as np
from scipy import stats
from dataclasses import dataclass
from typing import Tuple, Optional

@dataclass
class ABTestResult:
    """Results of an A/B test."""
    control_rate: float
    treatment_rate: float
    relative_lift: float
    absolute_lift: float
    z_statistic: float
    p_value: float
    confidence_interval: Tuple[float, float]
    is_significant: bool
    power: float
    
class ABTestAnalyzer:
    """
    Complete A/B testing framework with proper statistical methodology.
    """
    
    def __init__(self, alpha: float = 0.05, power_threshold: float = 0.80):
        self.alpha = alpha
        self.power_threshold = power_threshold
    
    def run_test(
        self, 
        control_successes: int, 
        control_total: int,
        treatment_successes: int, 
        treatment_total: int
    ) -> ABTestResult:
        """Run a two-proportion z-test."""
        
        # Calculate rates
        p_control = control_successes / control_total
        p_treatment = treatment_successes / treatment_total
        
        # Lifts
        absolute_lift = p_treatment - p_control
        relative_lift = (p_treatment - p_control) / p_control if p_control > 0 else 0
        
        # Pooled proportion
        p_pool = (control_successes + treatment_successes) / (control_total + treatment_total)
        
        # Standard error
        se = np.sqrt(p_pool * (1 - p_pool) * (1/control_total + 1/treatment_total))
        
        # Z-statistic
        z = absolute_lift / se if se > 0 else 0
        
        # P-value (two-tailed)
        p_value = 2 * (1 - stats.norm.cdf(abs(z)))
        
        # Confidence interval for the difference
        se_diff = np.sqrt(
            p_control * (1 - p_control) / control_total +
            p_treatment * (1 - p_treatment) / treatment_total
        )
        z_crit = stats.norm.ppf(1 - self.alpha / 2)
        ci = (absolute_lift - z_crit * se_diff, absolute_lift + z_crit * se_diff)
        
        # Power (approximate)
        power = self._calculate_power(p_control, p_treatment, min(control_total, treatment_total))
        
        return ABTestResult(
            control_rate=p_control,
            treatment_rate=p_treatment,
            relative_lift=relative_lift,
            absolute_lift=absolute_lift,
            z_statistic=z,
            p_value=p_value,
            confidence_interval=ci,
            is_significant=p_value < self.alpha,
            power=power
        )
    
    def _calculate_power(self, p1: float, p2: float, n: int) -> float:
        """Calculate statistical power."""
        effect = abs(p2 - p1)
        if effect == 0:
            return 0
        
        p_pool = (p1 + p2) / 2
        se_null = np.sqrt(2 * p_pool * (1 - p_pool) / n)
        se_alt = np.sqrt((p1 * (1 - p1) + p2 * (1 - p2)) / n)
        
        z_crit = stats.norm.ppf(1 - self.alpha / 2)
        z_power = (effect - z_crit * se_null) / se_alt
        
        return stats.norm.cdf(z_power)
    
    def required_sample_size(
        self, 
        baseline_rate: float, 
        minimum_detectable_effect: float,
        power: float = 0.80
    ) -> int:
        """Calculate required sample size per group."""
        
        p1 = baseline_rate
        p2 = baseline_rate * (1 + minimum_detectable_effect)
        effect = abs(p2 - p1)
        
        p_pool = (p1 + p2) / 2
        z_alpha = stats.norm.ppf(1 - self.alpha / 2)
        z_beta = stats.norm.ppf(power)
        
        var_null = 2 * p_pool * (1 - p_pool)
        var_alt = p1 * (1 - p1) + p2 * (1 - p2)
        
        n = ((z_alpha * np.sqrt(var_null) + z_beta * np.sqrt(var_alt)) / effect) ** 2
        
        return int(np.ceil(n))
    
    def print_report(self, result: ABTestResult, test_name: str = "A/B Test"):
        """Print a formatted test report."""
        
        print("\n" + "=" * 60)
        print(f"A/B TEST REPORT: {test_name}")
        print("=" * 60)
        
        print(f"\nConversion Rates:")
        print(f"  Control:   {result.control_rate:.2%}")
        print(f"  Treatment: {result.treatment_rate:.2%}")
        
        print(f"\nLift:")
        print(f"  Absolute: {result.absolute_lift:+.2%}")
        print(f"  Relative: {result.relative_lift:+.1%}")
        
        print(f"\nStatistical Analysis:")
        print(f"  Z-statistic: {result.z_statistic:.3f}")
        print(f"  P-value: {result.p_value:.4f}")
        print(f"  95% CI for difference: ({result.confidence_interval[0]:+.2%}, {result.confidence_interval[1]:+.2%})")
        
        print(f"\nTest Quality:")
        print(f"  Power: {result.power:.1%}")
        if result.power < self.power_threshold:
            print(f"  Warning: Low power. Consider larger sample size.")
        
        print(f"\nConclusion (α = {self.alpha}):")
        if result.is_significant:
            if result.absolute_lift > 0:
                print("  SIGNIFICANT: Treatment performs BETTER than control")
            else:
                print("  SIGNIFICANT: Treatment performs WORSE than control")
        else:
            print("  NOT SIGNIFICANT: Cannot conclude a difference exists")
            if result.power < self.power_threshold:
                print("  Note: Low power means we might be missing a real effect")
        
        print("=" * 60)


# Usage example
analyzer = ABTestAnalyzer(alpha=0.05)

# Test 1: Original example (not significant)
result1 = analyzer.run_test(
    control_successes=320, control_total=10000,
    treatment_successes=355, treatment_total=10000
)
analyzer.print_report(result1, "Checkout Button Color")

# Test 2: Larger sample (now significant!)
result2 = analyzer.run_test(
    control_successes=3200, control_total=100000,
    treatment_successes=3550, treatment_total=100000
)
analyzer.print_report(result2, "Checkout Button Color (Large Sample)")

# Calculate required sample size
n_required = analyzer.required_sample_size(
    baseline_rate=0.032,
    minimum_detectable_effect=0.10  # 10% relative improvement
)
print(f"\nTo detect 10% relative improvement with 80% power:")
print(f"Need {n_required:,} visitors per group")

# BAD: Stopping as soon as p < 0.05
# This inflates false positive rate to ~30%!

# GOOD: Pre-specify sample size and run to completion

# Testing 20 variations? Some will be "significant" by chance!

# Bonferroni correction
alpha_corrected = 0.05 / 20  # = 0.0025

# Or use False Discovery Rate (FDR) correction