A/B Testing Interview Questions

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.proportion import proportions_ztest
from statsmodels.stats.weightstats import ttest_ind as welch_ttest
import matplotlib.pyplot as plt

# Production: Comprehensive Hypothesis Testing Pipeline for A/B Tests

# ===== 1. Generate realistic A/B test data =====
np.random.seed(42)

# Example 1: Conversion rate test (proportions)
n_control = 10000
n_treatment = 10000
baseline_cvr = 0.10  # 10% baseline conversion
treatment_lift = 0.02  # +2pp absolute lift (20% relative)

control_conversions = np.random.binomial(n_control, baseline_cvr)
treatment_conversions = np.random.binomial(n_treatment, baseline_cvr + treatment_lift)

control_cvr = control_conversions / n_control
treatment_cvr = treatment_conversions / n_treatment

print("===== Conversion Rate Test (Proportions) =====")
print(f"Control: {control_conversions}/{n_control} = {control_cvr:.4f}")
print(f"Treatment: {treatment_conversions}/{n_treatment} = {treatment_cvr:.4f}")
print(f"Absolute lift: {(treatment_cvr - control_cvr):.4f}")
print(f"Relative lift: {100*(treatment_cvr/control_cvr - 1):.2f}%\n")

# ===== 2. Z-test for proportions (most common in A/B testing) =====
# Use when: Large samples (n*p > 5 and n*(1-p) > 5), binary outcome

counts = np.array([treatment_conversions, control_conversions])
nobs = np.array([n_treatment, n_control])

z_stat, p_value_two_sided = proportions_ztest(counts, nobs, alternative='two-sided')
_, p_value_one_sided = proportions_ztest(counts, nobs, alternative='larger')

print("--- Z-test for Proportions ---")
print(f"Z-statistic: {z_stat:.4f}")
print(f"p-value (two-sided): {p_value_two_sided:.6f}")
print(f"p-value (one-sided): {p_value_one_sided:.6f}")
print(f"Result: {'REJECT H0 (significant)' if p_value_two_sided < 0.05 else 'FAIL TO REJECT H0'}\n")

# ===== 3. Example 2: Revenue per user (continuous metric) =====
control_revenue = np.random.exponential(scale=20, size=5000)  # Right-skewed
treatment_revenue = np.random.exponential(scale=22, size=5000)  # +10% lift

print("===== Revenue Per User Test (Continuous) =====")
print(f"Control mean: ${control_revenue.mean():.2f} (median: ${np.median(control_revenue):.2f})")
print(f"Treatment mean: ${treatment_revenue.mean():.2f} (median: ${np.median(treatment_revenue):.2f})")
print(f"Absolute lift: ${(treatment_revenue.mean() - control_revenue.mean()):.2f}")
print(f"Relative lift: {100*(treatment_revenue.mean()/control_revenue.mean() - 1):.2f}%\n")

# ===== 4. T-test for continuous metrics =====
# Standard t-test (assumes equal variance)
t_stat_standard, p_value_standard = stats.ttest_ind(treatment_revenue, control_revenue)

# Welch's t-test (does NOT assume equal variance - PREFERRED)
t_stat_welch, p_value_welch, _ = welch_ttest(treatment_revenue, control_revenue, usevar='unequal')

print("--- T-test (Standard, assumes equal variance) ---")
print(f"t-statistic: {t_stat_standard:.4f}")
print(f"p-value: {p_value_standard:.6f}")

print("\n--- Welch's T-test (Unequal variance, PREFERRED) ---")
print(f"t-statistic: {t_stat_welch:.4f}")
print(f"p-value: {p_value_welch:.6f}")
print(f"Result: {'REJECT H0 (significant)' if p_value_welch < 0.05 else 'FAIL TO REJECT H0'}\n")

# ===== 5. Mann-Whitney U test (non-parametric alternative) =====
# Use when: Data is heavily skewed or violates normality assumptions
u_stat, p_value_mann = stats.mannwhitneyu(treatment_revenue, control_revenue, alternative='two-sided')

print("--- Mann-Whitney U Test (Non-parametric) ---")
print(f"U-statistic: {u_stat:.0f}")
print(f"p-value: {p_value_mann:.6f}")
print(f"Result: {'REJECT H0 (significant)' if p_value_mann < 0.05 else 'FAIL TO REJECT H0'}\n")

# ===== 6. Bootstrap confidence interval (robust alternative) =====
# Useful when: Unsure about distributions, want robust CI
def bootstrap_diff(control, treatment, n_iterations=10000):
    """Bootstrap confidence interval for difference in means"""
    diffs = []
    for _ in range(n_iterations):
        control_sample = np.random.choice(control, size=len(control), replace=True)
        treatment_sample = np.random.choice(treatment, size=len(treatment), replace=True)
        diffs.append(treatment_sample.mean() - control_sample.mean())
    return np.array(diffs)

bootstrap_diffs = bootstrap_diff(control_revenue, treatment_revenue)
ci_lower = np.percentile(bootstrap_diffs, 2.5)
ci_upper = np.percentile(bootstrap_diffs, 97.5)

print("--- Bootstrap 95% CI for Difference in Means ---")
print(f"95% CI: [${ci_lower:.2f}, ${ci_upper:.2f}]")
print(f"Observed diff: ${(treatment_revenue.mean() - control_revenue.mean()):.2f}")
print(f"Result: {'REJECT H0' if ci_lower > 0 else 'FAIL TO REJECT H0'}\n")

# ===== 7. Effect size calculation (practical significance) =====
# Cohen's d for continuous metrics
pooled_std = np.sqrt((control_revenue.var() + treatment_revenue.var()) / 2)
cohens_d = (treatment_revenue.mean() - control_revenue.mean()) / pooled_std

print("--- Effect Size (Cohen's d) ---")
print(f"Cohen's d: {cohens_d:.4f}")
print(f"Interpretation: ", end="")
if abs(cohens_d) < 0.2:
    print("Small effect")
elif abs(cohens_d) < 0.5:
    print("Medium effect")
else:
    print("Large effect")
print()

# ===== 8. Confidence interval for proportions =====
from statsmodels.stats.proportion import proportion_confint

ci_control = proportion_confint(control_conversions, n_control, alpha=0.05, method='wilson')
ci_treatment = proportion_confint(treatment_conversions, n_treatment, alpha=0.05, method='wilson')

print("--- 95% Confidence Intervals (Wilson method) ---")
print(f"Control CVR: {control_cvr:.4f} [{ci_control[0]:.4f}, {ci_control[1]:.4f}]")
print(f"Treatment CVR: {treatment_cvr:.4f} [{ci_treatment[0]:.4f}, {ci_treatment[1]:.4f}]")

# Confidence interval for difference in proportions
se_diff = np.sqrt(control_cvr*(1-control_cvr)/n_control + treatment_cvr*(1-treatment_cvr)/n_treatment)
diff = treatment_cvr - control_cvr
ci_diff_lower = diff - 1.96 * se_diff
ci_diff_upper = diff + 1.96 * se_diff

print(f"Difference: {diff:.4f} [{ci_diff_lower:.4f}, {ci_diff_upper:.4f}]")

import numpy as np
import pandas as pd
from statsmodels.stats.power import TTestIndPower, zt_ind_solve_power
from statsmodels.stats.proportion import proportion_effectsize
import matplotlib.pyplot as plt

# Production: Comprehensive Sample Size Calculator for A/B Tests

# ===== 1. Sample size for proportions (conversion rate) =====
def sample_size_proportions(baseline_rate, mde_relative, alpha=0.05, power=0.8):
    """
    Calculate sample size for proportion test

    Args:
        baseline_rate: Baseline conversion rate (e.g., 0.10 for 10%)
        mde_relative: Minimum detectable effect as relative lift (e.g., 0.05 for 5%)
        alpha: Significance level (Type I error rate)
        power: Statistical power (1 - Type II error rate)

    Returns:
        Sample size per group
    """
    # Convert relative MDE to absolute
    treatment_rate = baseline_rate * (1 + mde_relative)

    # Calculate effect size (Cohen's h for proportions)
    effect_size = proportion_effectsize(baseline_rate, treatment_rate)

    # Solve for sample size
    analysis = TTestIndPower()
    sample_size = analysis.solve_power(
        effect_size=effect_size,
        power=power,
        alpha=alpha,
        ratio=1.0,
        alternative='two-sided'
    )

    return int(np.ceil(sample_size))

# Example 1: Conversion rate test
baseline_cvr = 0.10  # 10% baseline
mde_relative = 0.05  # 5% relative lift (detect 10% → 10.5%)

n_per_group = sample_size_proportions(baseline_cvr, mde_relative)

print("===== Sample Size for Conversion Rate Test =====")
print(f"Baseline CVR: {baseline_cvr:.2%}")
print(f"MDE (relative): {mde_relative:.1%}")
print(f"MDE (absolute): {baseline_cvr * mde_relative:.4f}")
print(f"Significance (α): 0.05")
print(f"Power (1-β): 0.80")
print(f"Sample size per group: {n_per_group:,}")
print(f"Total sample size: {n_per_group * 2:,}\n")

# ===== 2. Sample size for continuous metrics (revenue) =====
def sample_size_continuous(mean, std, mde_relative, alpha=0.05, power=0.8):
    """
    Calculate sample size for continuous metric (e.g., revenue)

    Args:
        mean: Baseline mean (e.g., average revenue per user)
        std: Standard deviation
        mde_relative: Minimum detectable effect as relative lift
        alpha: Significance level
        power: Statistical power

    Returns:
        Sample size per group
    """
    # Absolute MDE
    mde_absolute = mean * mde_relative

    # Effect size (Cohen's d)
    effect_size = mde_absolute / std

    # Solve for sample size
    analysis = TTestIndPower()
    sample_size = analysis.solve_power(
        effect_size=effect_size,
        power=power,
        alpha=alpha,
        ratio=1.0,
        alternative='two-sided'
    )

    return int(np.ceil(sample_size))

# Example 2: Revenue per user test
avg_revenue = 50.0  # $50 average revenue
revenue_std = 30.0  # $30 standard deviation (CV = 60%, high variance)
mde_relative_revenue = 0.10  # Detect 10% lift

n_revenue = sample_size_continuous(avg_revenue, revenue_std, mde_relative_revenue)

print("===== Sample Size for Revenue Test =====")
print(f"Baseline revenue: ${avg_revenue:.2f}")
print(f"Std deviation: ${revenue_std:.2f} (CV = {revenue_std/avg_revenue:.1%})")
print(f"MDE (relative): {mde_relative_revenue:.1%}")
print(f"MDE (absolute): ${avg_revenue * mde_relative_revenue:.2f}")
print(f"Sample size per group: {n_revenue:,}\n")

# ===== 3. Sensitivity analysis: MDE vs sample size =====
baseline_cvr = 0.10
mde_values = [0.01, 0.02, 0.03, 0.05, 0.10, 0.15, 0.20]  # Relative MDE
sample_sizes = [sample_size_proportions(baseline_cvr, mde) for mde in mde_values]

print("===== Sensitivity Analysis: MDE vs Sample Size =====")
print(f"{'MDE (relative)':<15} {'MDE (absolute)':<15} {'Sample Size/Group':<20} {'Test Duration (days)':<20}")
print("-" * 70)

daily_traffic = 10000  # Assume 10k daily users

for mde, n in zip(mde_values, sample_sizes):
    mde_abs = baseline_cvr * mde
    duration = (n * 2) / daily_traffic
    print(f"{mde:>6.1%}          {mde_abs:>8.4f}          {n:>12,}          {duration:>10.1f}")

print()

# ===== 4. Tradeoff: Power vs Sample Size =====
power_values = [0.70, 0.75, 0.80, 0.85, 0.90, 0.95]
mde_fixed = 0.05  # 5% relative MDE

print("===== Power vs Sample Size Tradeoff (MDE = 5%) =====")
print(f"{'Power (1-β)':<12} {'Sample Size/Group':<20} {'% Increase vs 80%':<20}")
print("-" * 50)

baseline_n = sample_size_proportions(baseline_cvr, mde_fixed, power=0.80)

for pwr in power_values:
    n = sample_size_proportions(baseline_cvr, mde_fixed, power=pwr)
    increase = 100 * (n / baseline_n - 1)
    print(f"{pwr:>6.0%}       {n:>12,}          {increase:>+8.1f}%")

print()

# ===== 5. Unequal allocation (90/10 split) =====
# Sometimes you want 90% treatment, 10% control (faster rollout)
def sample_size_unequal(baseline_rate, mde_relative, ratio=9.0, alpha=0.05, power=0.8):
    """
    Calculate sample size with unequal allocation

    Args:
        ratio: Treatment size / Control size (e.g., 9.0 for 90/10 split)

    Returns:
        (control_size, treatment_size)
    """
    treatment_rate = baseline_rate * (1 + mde_relative)
    effect_size = proportion_effectsize(baseline_rate, treatment_rate)

    analysis = TTestIndPower()
    n_control = analysis.solve_power(
        effect_size=effect_size,
        power=power,
        alpha=alpha,
        ratio=ratio,
        alternative='two-sided'
    )

    n_treatment = n_control * ratio

    return int(np.ceil(n_control)), int(np.ceil(n_treatment))

n_control_90_10, n_treatment_90_10 = sample_size_unequal(baseline_cvr, mde_fixed, ratio=9.0)
n_equal = sample_size_proportions(baseline_cvr, mde_fixed)

print("===== Equal vs Unequal Allocation =====")
print(f"Equal (50/50): {n_equal:,} per group, {n_equal*2:,} total")
print(f"Unequal (10/90): Control = {n_control_90_10:,}, Treatment = {n_treatment_90_10:,}, Total = {n_control_90_10 + n_treatment_90_10:,}")
print(f"Cost of unequal allocation: {100*((n_control_90_10 + n_treatment_90_10)/(n_equal*2) - 1):.1f}% more users\n")

# ===== 6. Sequential testing adjustment =====
# If using sequential testing (early stopping), need ~10-20% more samples
n_sequential = int(n_per_group * 1.15)  # 15% inflation for sequential testing

print("===== Sequential Testing Adjustment =====")
print(f"Fixed-horizon sample size: {n_per_group:,}")
print(f"Sequential testing sample size: {n_sequential:,} (+15% inflation)")
print(f"Benefit: Can stop early if strong signal, but need buffer for Type I error control\n")

import numpy as np
import pandas as pd
from statsmodels.stats.power import TTestIndPower
from statsmodels.stats.proportion import proportion_effectsize
import matplotlib.pyplot as plt
from scipy import stats

# Production: Statistical Power Analysis and Visualization

# ===== 1. Power calculation for given sample size =====
def calculate_power(baseline_rate, treatment_rate, n_per_group, alpha=0.05):
    """
    Calculate statistical power for a proportion test

    Args:
        baseline_rate: Control conversion rate
        treatment_rate: Treatment conversion rate
        n_per_group: Sample size per group
        alpha: Significance level

    Returns:
        Statistical power (0 to 1)
    """
    effect_size = proportion_effectsize(baseline_rate, treatment_rate)
    analysis = TTestIndPower()
    power = analysis.solve_power(
        effect_size=effect_size,
        nobs1=n_per_group,
        alpha=alpha,
        ratio=1.0,
        alternative='two-sided'
    )
    return power

# Example: Power for different sample sizes
baseline_cvr = 0.10
treatment_cvr = 0.105  # 5% relative lift

print("===== Power Analysis: Effect of Sample Size =====")
print(f"Baseline: {baseline_cvr:.1%}, Treatment: {treatment_cvr:.1%} (5% relative lift)")
print(f"{'Sample Size/Group':<20} {'Power':<10} {'Interpretation':<30}")
print("-" * 60)

sample_sizes = [1000, 5000, 10000, 20000, 50000, 100000]
for n in sample_sizes:
    pwr = calculate_power(baseline_cvr, treatment_cvr, n)
    interp = "Too low" if pwr < 0.5 else "Low" if pwr < 0.8 else "Good" if pwr < 0.9 else "Excellent"
    print(f"{n:>12,}        {pwr:>6.1%}     {interp}")

print()

# ===== 2. Power curves: Visualize tradeoffs =====
# Power vs sample size for different effect sizes
sample_range = np.linspace(1000, 50000, 100)
effect_sizes = [0.01, 0.02, 0.05, 0.10]  # Relative lifts: 1%, 2%, 5%, 10%

print("===== Power Curves: Sample Size vs Power =====")
print("For different effect sizes (relative lift):")

for rel_lift in effect_sizes:
    treatment = baseline_cvr * (1 + rel_lift)
    powers = [calculate_power(baseline_cvr, treatment, int(n)) for n in sample_range]
    print(f"  {rel_lift:.0%} lift: Power reaches 80% at n = {sample_range[np.argmax(np.array(powers) >= 0.80)]:.0f}")

print()

# ===== 3. Minimum Detectable Effect (MDE) at fixed power =====
def calculate_mde(baseline_rate, n_per_group, power=0.8, alpha=0.05):
    """
    Calculate minimum detectable effect for given sample size and power

    Returns:
        MDE in absolute terms
    """
    analysis = TTestIndPower()

    # Binary search for MDE
    mde_low, mde_high = 0.0001, 0.50
    while mde_high - mde_low > 0.0001:
        mde_mid = (mde_low + mde_high) / 2
        treatment_rate = baseline_rate + mde_mid
        effect_size = proportion_effectsize(baseline_rate, treatment_rate)

        current_power = analysis.solve_power(
            effect_size=effect_size,
            nobs1=n_per_group,
            alpha=alpha,
            ratio=1.0
        )

        if current_power < power:
            mde_low = mde_mid
        else:
            mde_high = mde_mid

    return mde_high

print("===== MDE Analysis: What Can You Detect? =====")
print(f"Baseline CVR: {baseline_cvr:.1%}, Power: 80%, α: 0.05")
print(f"{'Sample Size/Group':<20} {'MDE (absolute)':<15} {'MDE (relative)':<15}")
print("-" * 50)

for n in [5000, 10000, 20000, 50000]:
    mde_abs = calculate_mde(baseline_cvr, n)
    mde_rel = mde_abs / baseline_cvr
    print(f"{n:>12,}        {mde_abs:>8.4f}        {mde_rel:>8.1%}")

print()

# ===== 4. Power vs significance level tradeoff =====
print("===== Power vs Significance Level Tradeoff =====")
print("Sample size: 20,000 per group, MDE: 5% relative lift")
print(f"{'Alpha (α)':<10} {'Power (1-β)':<12} {'Type I Error':<15} {'Type II Error':<15}")
print("-" * 50)

n_fixed = 20000
alpha_levels = [0.01, 0.025, 0.05, 0.10, 0.15]

for alpha in alpha_levels:
    pwr = calculate_power(baseline_cvr, baseline_cvr * 1.05, n_fixed, alpha=alpha)
    type1_error = alpha
    type2_error = 1 - pwr
    print(f"{alpha:>6.2f}      {pwr:>6.1%}        {type1_error:>6.1%}          {type2_error:>6.1%}")

print()

# ===== 5. Variance reduction impact on power =====
# CUPED and stratification can reduce variance by 20-50%
print("===== Impact of Variance Reduction on Power =====")
print("Sample size: 15,000 per group, MDE: 5% relative lift")
print(f"{'Variance Reduction':<25} {'Effective Sample Size':<25} {'Power':<10}")
print("-" * 60)

n_base = 15000
variance_reductions = [0.0, 0.20, 0.30, 0.40, 0.50]

for var_red in variance_reductions:
    # Variance reduction effectively increases sample size
    effective_n = n_base / (1 - var_red) if var_red < 1 else n_base * 2
    pwr = calculate_power(baseline_cvr, baseline_cvr * 1.05, int(effective_n))
    print(f"{var_red:>6.0%}                  {effective_n:>15,.0f}           {pwr:>6.1%}")

print()

# ===== 6. Post-hoc power analysis (retrospective) =====
# After experiment, calculate achieved power
def post_hoc_power(control_data, treatment_data, alpha=0.05):
    """
    Calculate observed power from experiment data

    Args:
        control_data: Array of control group data
        treatment_data: Array of treatment group data

    Returns:
        Observed effect size and power
    """
    # Calculate observed effect
    control_rate = control_data.mean()
    treatment_rate = treatment_data.mean()

    effect_size = proportion_effectsize(control_rate, treatment_rate)
    n_control = len(control_data)
    n_treatment = len(treatment_data)

    analysis = TTestIndPower()
    power = analysis.solve_power(
        effect_size=effect_size,
        nobs1=n_control,
        alpha=alpha,
        ratio=n_treatment / n_control
    )

    return {
        'control_rate': control_rate,
        'treatment_rate': treatment_rate,
        'effect_size': effect_size,
        'power': power,
        'n_control': n_control,
        'n_treatment': n_treatment
    }

# Simulate experiment data
np.random.seed(42)
control_sample = np.random.binomial(1, 0.10, size=8000)
treatment_sample = np.random.binomial(1, 0.105, size=8000)

post_hoc = post_hoc_power(control_sample, treatment_sample)

print("===== Post-Hoc Power Analysis (After Experiment) =====")
print(f"Control rate: {post_hoc['control_rate']:.4f}")
print(f"Treatment rate: {post_hoc['treatment_rate']:.4f}")
print(f"Observed lift: {100*(post_hoc['treatment_rate']/post_hoc['control_rate'] - 1):.2f}%")
print(f"Effect size (Cohen's h): {post_hoc['effect_size']:.4f}")
print(f"Sample size: {post_hoc['n_control']:,} per group")
print(f"Achieved power: {post_hoc['power']:.1%}")
print(f"Interpretation: {'Under-powered' if post_hoc['power'] < 0.8 else 'Adequately powered'}")

import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import chi2_contingency, chisquare
import matplotlib.pyplot as plt

# Production: SRM Detection and Debugging Pipeline

# ===== 1. Basic SRM Check (Chi-square test) =====
def check_srm(n_control, n_treatment, expected_ratio=0.5, alpha=0.001):
    """
    Check for Sample Ratio Mismatch using chi-square test

    Args:
        n_control: Number of users in control
        n_treatment: Number of users in treatment
        expected_ratio: Expected proportion for control (default 0.5 for 50/50)
        alpha: Significance level (default 0.001, very strict)

    Returns:
        Dictionary with SRM results
    """
    total = n_control + n_treatment
    expected_control = total * expected_ratio
    expected_treatment = total * (1 - expected_ratio)

    # Chi-square test
    observed = np.array([n_control, n_treatment])
    expected = np.array([expected_control, expected_treatment])

    chi2_stat = np.sum((observed - expected)**2 / expected)
    p_value = 1 - stats.chi2.cdf(chi2_stat, df=1)

    # Observed ratio
    observed_ratio = n_control / total

    return {
        'n_control': n_control,
        'n_treatment': n_treatment,
        'total': total,
        'expected_ratio': expected_ratio,
        'observed_ratio': observed_ratio,
        'chi2_stat': chi2_stat,
        'p_value': p_value,
        'srm_detected': p_value < alpha,
        'deviation_pct': 100 * abs(observed_ratio - expected_ratio) / expected_ratio
    }

# Example 1: No SRM (healthy experiment)
result_healthy = check_srm(n_control=10050, n_treatment=9950, expected_ratio=0.5)

print("===== Example 1: Healthy Experiment (No SRM) =====")
print(f"Control: {result_healthy['n_control']:,}, Treatment: {result_healthy['n_treatment']:,}")
print(f"Expected ratio: {result_healthy['expected_ratio']:.2f}")
print(f"Observed ratio: {result_healthy['observed_ratio']:.4f}")
print(f"Deviation: {result_healthy['deviation_pct']:.2f}%")
print(f"Chi-square statistic: {result_healthy['chi2_stat']:.4f}")
print(f"p-value: {result_healthy['p_value']:.6f}")
print(f"SRM Detected: {result_healthy['srm_detected']}\n")

# Example 2: SRM Detected (problematic experiment)
result_srm = check_srm(n_control=10500, n_treatment=9500, expected_ratio=0.5)

print("===== Example 2: SRM Detected (STOP EXPERIMENT) =====")
print(f"Control: {result_srm['n_control']:,}, Treatment: {result_srm['n_treatment']:,}")
print(f"Expected ratio: {result_srm['expected_ratio']:.2f}")
print(f"Observed ratio: {result_srm['observed_ratio']:.4f}")
print(f"Deviation: {result_srm['deviation_pct']:.2f}%")
print(f"Chi-square statistic: {result_srm['chi2_stat']:.4f}")
print(f"p-value: {result_srm['p_value']:.6f}")
print(f"SRM Detected: {'⚠️ YES - STOP AND DEBUG' if result_srm['srm_detected'] else 'No'}\n")

# ===== 2. SRM Check by Segment (debugging) =====
# Check SRM within different user segments to isolate root cause
def check_srm_by_segment(df, segment_col, expected_ratio=0.5):
    """
    Check SRM for each segment to identify where mismatch occurs

    Args:
        df: DataFrame with columns [segment_col, 'group']
        segment_col: Column name for segmentation (e.g., 'device', 'country')

    Returns:
        DataFrame with SRM results per segment
    """
    results = []

    for segment in df[segment_col].unique():
        segment_data = df[df[segment_col] == segment]
        n_control = (segment_data['group'] == 'control').sum()
        n_treatment = (segment_data['group'] == 'treatment').sum()

        srm_result = check_srm(n_control, n_treatment, expected_ratio)
        srm_result['segment'] = segment
        results.append(srm_result)

    return pd.DataFrame(results)

# Simulate experiment data with SRM in mobile users
np.random.seed(42)
n_users = 20000

# Desktop users: healthy 50/50 split
desktop = pd.DataFrame({
    'user_id': range(12000),
    'device': 'desktop',
    'group': np.random.choice(['control', 'treatment'], size=12000, p=[0.50, 0.50])
})

# Mobile users: SRM (bot filtering affected treatment more)
mobile = pd.DataFrame({
    'user_id': range(12000, 20000),
    'device': 'mobile',
    'group': np.random.choice(['control', 'treatment'], size=8000, p=[0.55, 0.45])  # SRM!
})

df_experiment = pd.concat([desktop, mobile], ignore_index=True)

print("===== SRM Check by Device Segment =====")
srm_by_device = check_srm_by_segment(df_experiment, 'device')
print(srm_by_device[['segment', 'n_control', 'n_treatment', 'observed_ratio', 'deviation_pct', 'p_value', 'srm_detected']].to_string(index=False))
print("\nDiagnosis: SRM detected in MOBILE users only → likely bot filtering or device compatibility issue\n")

# ===== 3. SRM Check Over Time (detect when SRM started) =====
def check_srm_over_time(df_time, date_col='date'):
    """
    Check SRM for each day to identify when mismatch began

    Args:
        df_time: DataFrame with columns [date_col, 'group']

    Returns:
        DataFrame with daily SRM results
    """
    results = []

    for date in sorted(df_time[date_col].unique()):
        day_data = df_time[df_time[date_col] == date]
        n_control = (day_data['group'] == 'control').sum()
        n_treatment = (day_data['group'] == 'treatment').sum()

        srm_result = check_srm(n_control, n_treatment)
        srm_result['date'] = date
        results.append(srm_result)

    return pd.DataFrame(results)

# Simulate 7-day experiment where SRM appears on day 5
dates = pd.date_range('2024-01-01', periods=7, freq='D')
df_daily = []

for i, date in enumerate(dates):
    if i < 4:
        # Days 1-4: healthy
        daily_users = pd.DataFrame({
            'user_id': range(i*3000, (i+1)*3000),
            'date': date,
            'group': np.random.choice(['control', 'treatment'], size=3000, p=[0.50, 0.50])
        })
    else:
        # Days 5-7: SRM (code change broke randomization)
        daily_users = pd.DataFrame({
            'user_id': range(i*3000, (i+1)*3000),
            'date': date,
            'group': np.random.choice(['control', 'treatment'], size=3000, p=[0.45, 0.55])
        })
    df_daily.append(daily_users)

df_daily = pd.concat(df_daily, ignore_index=True)

print("===== SRM Check Over Time (Daily) =====")
srm_daily = check_srm_over_time(df_daily)
print(srm_daily[['date', 'n_control', 'n_treatment', 'observed_ratio', 'p_value', 'srm_detected']].to_string(index=False))
print("\nDiagnosis: SRM started on 2024-01-05 → likely code deployment or config change on that date\n")

# ===== 4. Continuous SRM Monitoring (production dashboard) =====
def srm_alert_threshold(total_users, expected_ratio=0.5, alpha=0.001):
    """
    Calculate acceptable range for control users (no SRM)

    Returns:
        (lower_bound, upper_bound) for number of control users
    """
    expected_control = total_users * expected_ratio

    # Using normal approximation for large n
    se = np.sqrt(total_users * expected_ratio * (1 - expected_ratio))
    z = stats.norm.ppf(1 - alpha/2)  # Two-tailed

    lower = expected_control - z * se
    upper = expected_control + z * se

    return int(np.floor(lower)), int(np.ceil(upper))

print("===== SRM Alert Thresholds (Production Monitoring) =====")
print(f"{'Total Users':<15} {'Expected Control':<18} {'Acceptable Range (99.9% CI)':<35}")
print("-" * 70)

for total in [10000, 50000, 100000, 500000]:
    expected_ctrl = total * 0.5
    lower, upper = srm_alert_threshold(total)
    print(f"{total:>10,}      {expected_ctrl:>10,.0f}        [{lower:,} - {upper:,}]")

print()

# ===== 5. Common SRM Root Causes and Debugging =====
print("===== Common SRM Root Causes =====")
print("""
1. RANDOMIZATION BUGS
   - Hash collision (user ID % 2 not perfectly random)
   - Cookie issues (treatment cookie not set correctly)
   - Server-side vs client-side mismatch

2. BOT FILTERING
   - Bots filtered more in one group (e.g., bot detection flag added mid-experiment)
   - Anti-fraud rules differ by variant

3. BROWSER/DEVICE COMPATIBILITY
   - Treatment code crashes on certain browsers → users drop out
   - Mobile app version incompatibility

4. EXPERIMENT INTERACTION
   - Another experiment running affects assignment
   - Triggered experiments (only some users eligible)

5. DATA PIPELINE ISSUES
   - Logging inconsistency between groups
   - Data loss in one variant

DEBUGGING STEPS:
1. Check SRM by segment (device, browser, country, day)
2. Review code changes around SRM start date
3. Verify randomization logic (run A/A test)
4. Check bot filtering and fraud rules
5. Validate experiment eligibility criteria
""")

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.power import TTestIndPower
from statsmodels.stats.proportion import proportions_ztest, proportion_effectsize
import matplotlib.pyplot as plt

# Production: Type I/II Error Analysis and Simulation

# ===== 1. Simulate Type I Errors (False Positives) =====
# Run 1000 A/A tests (no real effect) and count false positives
def simulate_type1_errors(n_simulations=1000, n_per_group=5000, alpha=0.05):
    """
    Simulate Type I error rate with A/A tests (H0 is true)

    Returns:
        Empirical Type I error rate
    """
    np.random.seed(42)
    false_positives = 0
    baseline_rate = 0.10

    for _ in range(n_simulations):
        # Both groups same (H0 true)
        control = np.random.binomial(1, baseline_rate, size=n_per_group)
        treatment = np.random.binomial(1, baseline_rate, size=n_per_group)

        # Test
        _, p_value = proportions_ztest(
            [treatment.sum(), control.sum()],
            [len(treatment), len(control)]
        )

        if p_value < alpha:
            false_positives += 1

    return false_positives / n_simulations

type1_rate = simulate_type1_errors(n_simulations=1000)

print("===== Type I Error Simulation (A/A Tests) =====")
print(f"Expected Type I error rate (α): 5.0%")
print(f"Observed Type I error rate: {100*type1_rate:.1f}%")
print(f"Result: {'✓ Matches expected' if abs(type1_rate - 0.05) < 0.01 else '✗ Calibration issue'}")
print(f"\nInterpretation: In {1000} A/A tests, we falsely rejected H0 {int(type1_rate*1000)} times\n")

# ===== 2. Simulate Type II Errors (False Negatives) =====
# Run tests with real effect, count how often we miss it
def simulate_type2_errors(n_simulations=1000, n_per_group=5000,
                          baseline_rate=0.10, effect_size_rel=0.05, alpha=0.05):
    """
    Simulate Type II error rate with A/B tests (H0 is false, H1 is true)

    Returns:
        Empirical Type II error rate and power
    """
    np.random.seed(42)
    true_positives = 0
    treatment_rate = baseline_rate * (1 + effect_size_rel)

    for _ in range(n_simulations):
        # Real effect exists (H1 true)
        control = np.random.binomial(1, baseline_rate, size=n_per_group)
        treatment = np.random.binomial(1, treatment_rate, size=n_per_group)

        # Test
        _, p_value = proportions_ztest(
            [treatment.sum(), control.sum()],
            [len(treatment), len(control)]
        )

        if p_value < alpha:
            true_positives += 1

    power = true_positives / n_simulations
    type2_rate = 1 - power

    return type2_rate, power

type2_rate, observed_power = simulate_type2_errors(
    n_simulations=1000,
    n_per_group=5000,
    effect_size_rel=0.05
)

print("===== Type II Error Simulation (A/B Tests with 5% Lift) =====")
print(f"Sample size: 5,000 per group")
print(f"Effect: 10% → 10.5% (5% relative lift)")
print(f"Observed power (1-β): {100*observed_power:.1f}%")
print(f"Observed Type II error rate (β): {100*type2_rate:.1f}%")
print(f"\nInterpretation: In {1000} tests with real 5% lift, we detected it {int(observed_power*1000)} times\n")

# ===== 3. Type I vs Type II Tradeoff =====
# Show how changing α affects β (for fixed sample size)
print("===== Type I vs Type II Tradeoff =====")
print(f"Sample size: 10,000 per group, True effect: 5% relative lift")
print(f"{'Alpha (α)':<12} {'Type I Error':<15} {'Type II Error (β)':<20} {'Power (1-β)':<15}")
print("-" * 65)

n_fixed = 10000
baseline = 0.10
treatment = baseline * 1.05

for alpha in [0.01, 0.025, 0.05, 0.10, 0.15, 0.20]:
    effect = proportion_effectsize(baseline, treatment)
    analysis = TTestIndPower()
    power = analysis.solve_power(
        effect_size=effect,
        nobs1=n_fixed,
        alpha=alpha,
        ratio=1.0
    )
    beta = 1 - power
    print(f"{alpha:>6.2f}       {alpha:>6.1%}          {beta:>8.1%}             {power:>8.1%}")

print("\nKey insight: Lower α (fewer false positives) → Higher β (more false negatives)\n")

# ===== 4. Business Cost Analysis =====
# Calculate expected cost of Type I and Type II errors
def calculate_error_costs(
    baseline_cvr=0.10,
    effect_size=0.05,
    n_per_group=10000,
    alpha=0.05,
    cost_type1=100000,  # Cost of shipping bad feature
    cost_type2=50000,   # Opportunity cost of missing good feature
    tests_per_year=100
):
    """
    Calculate annual cost of Type I and Type II errors

    Returns:
        Dictionary with cost analysis
    """
    # Type I errors (false positives)
    # Occur in tests where H0 is true (no real effect)
    # Assume 50% of tests have no effect
    prob_h0_true = 0.50
    expected_type1_per_year = tests_per_year * prob_h0_true * alpha

    # Type II errors (false negatives)
    # Occur in tests where H1 is true (real effect exists)
    treatment_cvr = baseline_cvr * (1 + effect_size)
    effect = proportion_effectsize(baseline_cvr, treatment_cvr)
    analysis = TTestIndPower()
    power = analysis.solve_power(
        effect_size=effect,
        nobs1=n_per_group,
        alpha=alpha,
        ratio=1.0
    )
    beta = 1 - power

    prob_h1_true = 0.50  # Assume 50% of tests have real effect
    expected_type2_per_year = tests_per_year * prob_h1_true * beta

    # Annual costs
    annual_type1_cost = expected_type1_per_year * cost_type1
    annual_type2_cost = expected_type2_per_year * cost_type2
    total_cost = annual_type1_cost + annual_type2_cost

    return {
        'tests_per_year': tests_per_year,
        'expected_type1_errors': expected_type1_per_year,
        'expected_type2_errors': expected_type2_per_year,
        'annual_type1_cost': annual_type1_cost,
        'annual_type2_cost': annual_type2_cost,
        'total_annual_cost': total_cost,
        'power': power,
        'alpha': alpha
    }

cost_analysis = calculate_error_costs(
    n_per_group=10000,
    alpha=0.05,
    cost_type1=100000,
    cost_type2=50000
)

print("===== Business Cost of Type I and Type II Errors =====")
print(f"Tests per year: {cost_analysis['tests_per_year']}")
print(f"Sample size: 10,000 per group")
print(f"Power: {100*cost_analysis['power']:.1f}%\n")

print(f"Expected Type I errors/year: {cost_analysis['expected_type1_errors']:.1f}")
print(f"Cost per Type I error: ${cost_analysis['annual_type1_cost']/cost_analysis['expected_type1_errors']:,.0f}")
print(f"Annual Type I cost: ${cost_analysis['annual_type1_cost']:,.0f}\n")

print(f"Expected Type II errors/year: {cost_analysis['expected_type2_errors']:.1f}")
print(f"Cost per Type II error: ${cost_analysis['annual_type2_cost']/cost_analysis['expected_type2_errors']:,.0f}")
print(f"Annual Type II cost: ${cost_analysis['annual_type2_cost']:,.0f}\n")

print(f"Total annual cost: ${cost_analysis['total_annual_cost']:,.0f}\n")

# ===== 5. Optimal alpha selection based on costs =====
print("===== Optimal Alpha Based on Business Costs =====")
print(f"{'Alpha':<8} {'Type I/year':<12} {'Type II/year':<12} {'Total Cost':<15} {'Power':<8}")
print("-" * 60)

for alpha_test in [0.01, 0.05, 0.10, 0.15]:
    result = calculate_error_costs(
        n_per_group=10000,
        alpha=alpha_test,
        cost_type1=100000,
        cost_type2=50000
    )
    print(f"{alpha_test:>6.2f}   {result['expected_type1_errors']:>6.1f}       "
          f"{result['expected_type2_errors']:>6.1f}       ${result['total_annual_cost']:>10,.0f}     {100*result['power']:>5.1f}%")

print("\nKey insight: If Type I cost >> Type II cost, use lower α (0.01)")
print("            If Type II cost >> Type I cost, use higher α (0.10)\n")

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.proportion import proportions_ztest
import matplotlib.pyplot as plt

# Production: Peeking Problem Analysis and Sequential Testing

# ===== 1. Simulate the peeking problem =====
def simulate_peeking_inflation(
    n_simulations=1000,
    total_users_per_group=10000,
    n_peeks=14,  # Check daily for 14 days
    baseline_rate=0.10,
    alpha=0.05
):
    """
    Simulate how peeking inflates false positive rate

    Args:
        n_simulations: Number of A/A tests to run
        total_users_per_group: Final sample size per group
        n_peeks: Number of times to check p-value
        baseline_rate: True conversion rate (same for both groups in A/A test)

    Returns:
        Empirical false positive rate with peeking
    """
    np.random.seed(42)
    false_positives_peeking = 0
    false_positives_no_peek = 0

    users_per_peek = total_users_per_group // n_peeks

    for sim in range(n_simulations):
        # Generate data (A/A test, both groups identical)
        control_all = np.random.binomial(1, baseline_rate, size=total_users_per_group)
        treatment_all = np.random.binomial(1, baseline_rate, size=total_users_per_group)

        # Strategy 1: PEEKING (check at each peek, stop if significant)
        stopped_early = False
        for peek in range(1, n_peeks + 1):
            n_so_far = peek * users_per_peek
            control_peek = control_all[:n_so_far]
            treatment_peek = treatment_all[:n_so_far]

            _, p_value = proportions_ztest(
                [treatment_peek.sum(), control_peek.sum()],
                [len(treatment_peek), len(control_peek)]
            )

            if p_value < alpha:
                false_positives_peeking += 1
                stopped_early = True
                break

        # Strategy 2: NO PEEKING (only check at end)
        _, p_value_final = proportions_ztest(
            [treatment_all.sum(), control_all.sum()],
            [len(treatment_all), len(control_all)]
            )

        if p_value_final < alpha:
            false_positives_no_peek += 1

    return {
        'peeking_fpr': false_positives_peeking / n_simulations,
        'no_peeking_fpr': false_positives_no_peek / n_simulations,
        'inflation_factor': (false_positives_peeking / n_simulations) / alpha
    }

peeking_result = simulate_peeking_inflation(n_simulations=1000, n_peeks=14)

print("===== Peeking Problem: False Positive Rate Inflation =====")
print(f"Expected FPR (α): 5.0%")
print(f"FPR without peeking: {100*peeking_result['no_peeking_fpr']:.1f}%")
print(f"FPR with peeking (14 checks): {100*peeking_result['peeking_fpr']:.1f}%")
print(f"Inflation factor: {peeking_result['inflation_factor']:.2f}x")
print(f"\n⚠️ Peeking inflates FPR by {100*(peeking_result['peeking_fpr']/0.05 - 1):.0f}%!\n")

# ===== 2. Peeking vs number of looks =====
print("===== Impact of Number of Peeks =====")
print(f"{'Number of Peeks':<18} {'FPR (peeking)':<15} {'Inflation':<12}")
print("-" * 45)

for n_peeks in [1, 3, 7, 14, 30]:
    result = simulate_peeking_inflation(n_simulations=500, n_peeks=n_peeks)
    print(f"{n_peeks:>10}         {100*result['peeking_fpr']:>8.1f}%        {result['inflation_factor']:>6.2f}x")

print()

# ===== 3. Solution 1: Sequential Testing with Alpha Spending =====
# O'Brien-Fleming boundaries (conservative early stopping)
def obrien_fleming_boundary(n_looks, alpha=0.05):
    """
    Calculate O'Brien-Fleming alpha spending boundaries

    Args:
        n_looks: Number of planned interim looks
        alpha: Overall significance level

    Returns:
        List of critical z-values for each look
    """
    # O'Brien-Fleming: spend very little alpha early, most at the end
    # z_k = z_alpha * sqrt(K/k) where K = total looks, k = current look
    z_alpha = stats.norm.ppf(1 - alpha/2)  # Two-sided

    boundaries = []
    for k in range(1, n_looks + 1):
        z_boundary = z_alpha * np.sqrt(n_looks / k)
        boundaries.append(z_boundary)

    return boundaries

# Example: 4 interim looks + final analysis (5 total)
n_looks = 5
of_boundaries = obrien_fleming_boundary(n_looks, alpha=0.05)

print("===== O'Brien-Fleming Boundaries (Sequential Testing) =====")
print(f"Overall α: 0.05, Number of looks: {n_looks}")
print(f"{'Look':<8} {'% Complete':<15} {'Critical Z':<12} {'Equivalent α':<15}")
print("-" * 50)

for i, z_bound in enumerate(of_boundaries, 1):
    pct_complete = (i / n_looks) * 100
    equivalent_alpha = 2 * (1 - stats.norm.cdf(z_bound))  # Two-sided
    print(f"{i:>4}     {pct_complete:>6.0f}%          {z_bound:>6.3f}        {equivalent_alpha:>8.5f}")

print("\nInterpretation: Early looks require very strong evidence (Z > 4.0)")
print("                Final look uses standard threshold (Z > 1.96)\n")

# ===== 4. Solution 2: Pocock Boundaries (less conservative) =====
def pocock_boundary(n_looks, alpha=0.05):
    """
    Calculate Pocock alpha spending boundaries (constant across looks)

    Returns:
        Constant critical z-value for all looks
    """
    # Pocock: equal alpha spending at each look
    # Approximate critical value (exact requires numerical integration)
    if n_looks == 2:
        z_boundary = 2.178
    elif n_looks == 3:
        z_boundary = 2.289
    elif n_looks == 4:
        z_boundary = 2.361
    elif n_looks == 5:
        z_boundary = 2.413
    else:
        # Approximation for other values
        z_boundary = stats.norm.ppf(1 - alpha/(2*n_looks)) * 1.2

    return z_boundary

pocock_z = pocock_boundary(n_looks=5)

print("===== Pocock Boundaries (Uniform Alpha Spending) =====")
print(f"Constant critical Z for all {n_looks} looks: {pocock_z:.3f}")
print(f"Equivalent α per look: {2*(1 - stats.norm.cdf(pocock_z)):.5f}")
print(f"\nComparison: O'Brien-Fleming is more conservative early, Pocock allows earlier stopping\n")

# ===== 5. Always-Valid Inference (Anytime p-values) =====
# Using mSPRT (mixture Sequential Probability Ratio Test)
def always_valid_confidence_sequence(
    control_successes, control_total,
    treatment_successes, treatment_total,
    alpha=0.05
):
    """
    Calculate always-valid confidence interval (can peek anytime)

    Uses sequential testing theory to provide valid inference at any time
    """
    # Simple approach: adjust alpha using a conservative bound
    # (Production systems use more sophisticated methods like mSPRT)
    adjusted_alpha = alpha / np.sqrt(control_total + treatment_total)

    control_rate = control_successes / control_total
    treatment_rate = treatment_successes / treatment_total

    # Wider confidence intervals to maintain validity
    se = np.sqrt(
        control_rate * (1 - control_rate) / control_total +
        treatment_rate * (1 - treatment_rate) / treatment_total
    )

    z_adjusted = stats.norm.ppf(1 - adjusted_alpha/2)
    diff = treatment_rate - control_rate
    ci_lower = diff - z_adjusted * se
    ci_upper = diff + z_adjusted * se

    return {
        'diff': diff,
        'ci_lower': ci_lower,
        'ci_upper': ci_upper,
        'is_significant': ci_lower > 0 or ci_upper < 0
    }

# Example: Check after 5000 users per group
av_result = always_valid_confidence_sequence(
    control_successes=500,
    control_total=5000,
    treatment_successes=525,
    treatment_total=5000
)

print("===== Always-Valid Inference (Anytime p-values) =====")
print(f"Control: 500/5000 = 10.00%")
print(f"Treatment: 525/5000 = 10.50%")
print(f"Difference: {100*av_result['diff']:.2f}%")
print(f"Always-valid 95% CI: [{100*av_result['ci_lower']:.2f}%, {100*av_result['ci_upper']:.2f}%]")
print(f"Significant: {av_result['is_significant']}")
print(f"\nBenefit: Can check progress anytime without inflating FPR\n")

# ===== 6. Production Recommendations =====
print("===== Production Guidelines for Handling Peeking =====")
print("""
PROBLEM:
- Checking p-value daily for 14 days inflates FPR from 5% → 29%
- Stopping when p < 0.05 invalidates statistical guarantees

SOLUTIONS:

1. FIXED-HORIZON (Simplest)
   - Pre-commit to sample size (e.g., 20,000 per group)
   - Don't look at results until reached
   - Check p-value once at end
   - Pros: Simple, maintains α = 0.05
   - Cons: Can't stop early, slow learning

2. SEQUENTIAL TESTING (O'Brien-Fleming)
   - Plan interim looks (e.g., 25%, 50%, 75%, 100%)
   - Use adjusted alpha boundaries (e.g., p < 0.0001 at 25%, p < 0.05 at 100%)
   - Stop early only if boundary crossed
   - Pros: Valid early stopping, maintains α = 0.05
   - Cons: Requires planning, implementation complexity

3. BAYESIAN MONITORING
   - Use posterior probability: P(treatment better | data)
   - Stop when P(treatment > control) > 95%
   - No alpha inflation issue
   - Pros: Continuous monitoring OK, intuitive
   - Cons: Requires priors, different interpretation

4. ALWAYS-VALID INFERENCE (Advanced)
   - Use confidence sequences (valid at any time)
   - Can peek freely without penalty
   - Pros: Maximum flexibility
   - Cons: Wider confidence intervals, cutting-edge research

RECOMMENDATION BY TEAM SIZE:
- Small team (<10 eng): Fixed-horizon (wait 2 weeks)
- Medium team (10-100): Sequential testing with 2-3 interim looks
- Large team (100+): Bayesian or always-valid with real-time dashboards
""")

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

# Production: CUPED Variance Reduction

# ===== 1. Generate experiment data =====
np.random.seed(42)
n = 5000

# Pre-experiment revenue (7-30 days before)
pre_control = np.random.exponential(50, n)
pre_treatment = np.random.exponential(50, n)

# Experiment revenue (correlated r ~ 0.7)
exp_control = 0.7*pre_control + np.random.normal(30, 20, n)
exp_treatment = 0.7*pre_treatment + np.random.normal(35, 20, n)  # +$5 lift

# ===== 2. Standard analysis (no CUPED) =====
mean_diff = exp_treatment.mean() - exp_control.mean()
se_standard = np.sqrt(exp_control.var()/n + exp_treatment.var()/n)
t_standard, p_standard = stats.ttest_ind(exp_treatment, exp_control)

print(f"Standard: Diff=${mean_diff:.2f}, SE=${se_standard:.2f}, p={p_standard:.4f}")

# ===== 3. CUPED analysis =====
X = np.concatenate([pre_control, pre_treatment])
Y = np.concatenate([exp_control, exp_treatment])

theta = np.cov(X, Y)[0,1] / np.var(X)
Y_cuped = Y - theta * (X - X.mean())

exp_control_cuped = Y_cuped[:n]
exp_treatment_cuped = Y_cuped[n:]

se_cuped = np.sqrt(exp_control_cuped.var()/n + exp_treatment_cuped.var()/n)
t_cuped, p_cuped = stats.ttest_ind(exp_treatment_cuped, exp_control_cuped)

var_reduction = 1 - (np.var(Y_cuped) / np.var(Y))
print(f"CUPED: Diff=${mean_diff:.2f}, SE=${se_cuped:.2f}, p={p_cuped:.4f}")
print(f"Variance reduction: {100*var_reduction:.1f}% (SE reduced {100*(1-se_cuped/se_standard):.1f}%)")

import pymc3 as pm

with pm.Model():
    p_control = pm.Beta('p_C', 1, 1)
    p_treatment = pm.Beta('p_T', 1, 1)

    obs_C = pm.Binomial('obs_C', n=n_C, p=p_control, observed=success_C)
    obs_T = pm.Binomial('obs_T', n=n_T, p=p_treatment, observed=success_T)

    trace = pm.sample(1000)

┌────────────────────────────────────────────────────────┐
│           SWITCHBACK EXPERIMENT DESIGN                 │
├────────────────────────────────────────────────────────┤
│                                                         │
│  TIME                                                  │
│   0hr    1hr    2hr    3hr    4hr    5hr    6hr       │
│   │      │      │      │      │      │      │         │
│   ↓      ↓      ↓      ↓      ↓      ↓      ↓         │
│  ┌───┐  ┌───┐  ┌───┐  ┌───┐  ┌───┐  ┌───┐  ┌───┐    │
│  │ C │  │ T │  │ C │  │ T │  │ C │  │ T │  │ C │    │
│  └───┘  └───┘  └───┘  └───┘  └───┘  └───┘  └───┘    │
│                                                         │
│  ALL users in market see same treatment at same time  │
│                                                         │
│  ANALYSIS:                                            │
│  - Compare T periods vs C periods                     │
│  - Account for time-of-day effects                    │
│  - Account for serial correlation                     │
│  - Washout period between switches                    │
│                                                         │
└────────────────────────────────────────────────────────┘

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant
from datetime import datetime, timedelta

class SwitchbackAnalyzer:
    \"\"\"Analyze switchback experiments with time-based randomization.\"\"\"

    def __init__(self, alpha=0.05):
        self.alpha = alpha

    def generate_switchback_schedule(self, start_time, n_periods, period_hours=1):
        \"\"\"
        Generate alternating treatment schedule.
        \"\"\"
        schedule = []
        current = start_time

        for i in range(n_periods):
            treatment = i % 2  # Alternate 0, 1, 0, 1...
            schedule.append({
                'period': i,
                'start_time': current,
                'end_time': current + timedelta(hours=period_hours),
                'treatment': treatment
            })
            current += timedelta(hours=period_hours)

        return pd.DataFrame(schedule)

    def analyze_switchback(self, df, outcome_col, treatment_col, period_col):
        \"\"\"
        Analyze switchback with period-level aggregation.
        Accounts for time-of-day effects.
        \"\"\"
        # Aggregate to period level
        period_summary = df.groupby([period_col, treatment_col])[outcome_col].agg([
            'mean', 'std', 'count'
        ]).reset_index()

        control_periods = period_summary[period_summary[treatment_col]==0]
        treatment_periods = period_summary[period_summary[treatment_col]==1]

        # Paired t-test (periods alternate)
        effect = treatment_periods['mean'].mean() - control_periods['mean'].mean()

        # Use period-level data
        t_stat, p_value = stats.ttest_ind(
            treatment_periods['mean'], 
            control_periods['mean']
        )

        se = effect / t_stat if t_stat != 0 else np.nan

        return {
            'effect': effect,
            'se': se,
            'pvalue': p_value,
            'n_periods_control': len(control_periods),
            'n_periods_treatment': len(treatment_periods)
        }

    def adjust_for_time_effects(self, df, outcome_col, treatment_col, hour_col):
        \"\"\"
        Regression adjustment for hour-of-day effects.
        Y = β0 + β1*Treatment + β2*Hour + ε
        \"\"\"
        df = df.copy()
        df = pd.get_dummies(df, columns=[hour_col], prefix='hour', drop_first=True)

        # Model with time controls
        hour_cols = [c for c in df.columns if c.startswith('hour_')]
        X = add_constant(df[[treatment_col] + hour_cols])
        y = df[outcome_col]

        model = OLS(y, X).fit(cov_type='HC3')  # Robust SE

        return {
            'effect': model.params[treatment_col],
            'se': model.bse[treatment_col],
            'pvalue': model.pvalues[treatment_col],
            'ci': model.conf_int().loc[treatment_col].values
        }

# Example: Uber surge pricing switchback
np.random.seed(42)

print("="*70)
print("UBER - SURGE PRICING SWITCHBACK EXPERIMENT")
print("="*70)

analyzer = SwitchbackAnalyzer()

# Generate 48 hour schedule (hourly switches)
start = datetime(2025, 6, 1, 0, 0)
schedule = analyzer.generate_switchback_schedule(start, n_periods=48, period_hours=1)

print(f"\\nExperiment setup:")
print(f"  Duration: 48 hours (2 days)")
print(f"  Switch frequency: Every 1 hour")
print(f"  Total periods: {len(schedule)}")
print(f"  Control periods: {(schedule['treatment']==0).sum()}")
print(f"  Treatment periods: {(schedule['treatment']==1).sum()}")

# Simulate outcomes
# Baseline varies by hour (peak hours have more rides)
# Treatment effect: +12% rides (surge increases supply)

data = []
for _, period in schedule.iterrows():
    hour = period['start_time'].hour

    # Hour-of-day effect (peak hours)
    if 7 <= hour <= 9 or 17 <= hour <= 19:
        baseline = 500  # Peak
    elif 22 <= hour or hour <= 5:
        baseline = 100  # Late night
    else:
        baseline = 300  # Normal

    # Treatment effect
    if period['treatment'] == 1:
        mean_rides = baseline * 1.12  # +12%
    else:
        mean_rides = baseline

    # Simulate rides in this period
    n_rides = int(np.random.normal(mean_rides, mean_rides * 0.15))

    data.append({
        'period': period['period'],
        'hour': hour,
        'treatment': period['treatment'],
        'rides': n_rides
    })

df = pd.DataFrame(data)

# Naive analysis
naive_control = df[df['treatment']==0]['rides'].mean()
naive_treatment = df[df['treatment']==1]['rides'].mean()
naive_effect = naive_treatment - naive_control
naive_lift = naive_effect / naive_control * 100

print(f"\\nNaive analysis:")
print(f"  Control mean: {naive_control:.1f} rides/hour")
print(f"  Treatment mean: {naive_treatment:.1f} rides/hour")
print(f"  Effect: {naive_effect:+.1f} rides/hour")
print(f"  Lift: {naive_lift:+.1f}%")

# Switchback analysis (period-level)
switchback_result = analyzer.analyze_switchback(df, 'rides', 'treatment', 'period')

print(f"\\nSwitchback analysis (period-level):")
print(f"  Effect: {switchback_result['effect']:+.1f} rides/hour")
print(f"  SE: {switchback_result['se']:.1f}")
print(f"  P-value: {switchback_result['pvalue']:.4f}")
print(f"  Periods: {switchback_result['n_periods_control']} control, {switchback_result['n_periods_treatment']} treatment")

# Adjusted for hour-of-day
adjusted_result = analyzer.adjust_for_time_effects(df, 'rides', 'treatment', 'hour')

print(f"\\nHour-adjusted analysis:")
print(f"  Effect: {adjusted_result['effect']:+.1f} rides/hour")
print(f"  SE: {adjusted_result['se']:.1f}")
print(f"  P-value: {adjusted_result['pvalue']:.4f}")
print(f"  95% CI: ({adjusted_result['ci'][0]:.1f}, {adjusted_result['ci'][1]:.1f})")

print(f"\\nTrue effect (simulation): +12.0% lift")
print(f"Measured lift: {(adjusted_result['effect']/naive_control)*100:+.1f}%")

print("="*70)

# Output:
# ======================================================================
# UBER - SURGE PRICING SWITCHBACK EXPERIMENT
# ======================================================================
# 
# Experiment setup:
#   Duration: 48 hours (2 days)
#   Switch frequency: Every 1 hour
#   Total periods: 48
#   Control periods: 24
#   Treatment periods: 24
# 
# Naive analysis:
#   Control mean: 288.5 rides/hour
#   Treatment mean: 322.5 rides/hour
#   Effect: +34.0 rides/hour
#   Lift: +11.8%
# 
# Switchback analysis (period-level):
#   Effect: +34.0 rides/hour
#   SE: 17.2
#   P-value: 0.0537
# 
# Hour-adjusted analysis:
#   Effect: +34.1 rides/hour
#   SE: 6.8
#   P-value: 0.0000
#   95% CI: (20.7, 47.5)
# 
# True effect (simulation): +12.0% lift
# Measured lift: +11.8%
# ======================================================================