A/B Testing in ML: Statistically Rigorous Experiments in Production

import math
from scipy import stats

def required_sample_size(
    baseline_rate: float,
    mde: float,
    alpha: float = 0.05,
    power: float = 0.80
) -> int:
    """
    Compute required sample size PER GROUP for a two-proportion z-test.

    This must be computed BEFORE the experiment starts. Running an
    experiment without pre-computed sample size is guessing, not testing.

    Args:
        baseline_rate: control group metric (e.g., 0.05 for 5% CTR)
        mde: minimum detectable effect as absolute difference (e.g., 0.005
             for 0.5 percentage points)
        alpha: significance level — probability of false positive (default 0.05)
        power: probability of detecting a real effect (default 0.80)

    Returns:
        Required number of users per group (control and treatment each)
    """
    p1 = baseline_rate
    p2 = baseline_rate + mde

    z_alpha = stats.norm.ppf(1 - alpha / 2)  # two-sided test
    z_beta = stats.norm.ppf(power)

    # Variance under each hypothesis
    numerator = (z_alpha + z_beta) ** 2 * (p1 * (1 - p1) + p2 * (1 - p2))
    denominator = (p2 - p1) ** 2

    return math.ceil(numerator / denominator)


def experiment_plan(
    baseline_rate: float,
    mde: float,
    daily_users: int,
    alpha: float = 0.05,
    power: float = 0.80
) -> dict:
    """
    Generate a complete experiment plan with timeline and guardrails.
    """
    n_per_group = required_sample_size(baseline_rate, mde, alpha, power)
    total_users = n_per_group * 2
    min_days_for_power = math.ceil(total_users / daily_users)

    # Enforce minimum duration rules
    min_business_cycles = 14  # 2 full weeks
    novelty_buffer = 7        # 1 additional week
    min_duration = max(min_days_for_power, min_business_cycles + novelty_buffer)

    return {
        "n_per_group": n_per_group,
        "total_users_needed": total_users,
        "min_days_for_statistical_power": min_days_for_power,
        "min_days_for_business_cycles": min_business_cycles,
        "novelty_buffer_days": novelty_buffer,
        "recommended_duration_days": min_duration,
        "alpha": alpha,
        "power": power,
        "baseline_rate": baseline_rate,
        "mde": mde,
        "note": "Do NOT stop early even if p < 0.05 before this duration."
    }


# --- Example: plan an experiment for a 5% baseline CTR ---
plan = experiment_plan(
    baseline_rate=0.05,
    mde=0.005,          # detect a 0.5 percentage point lift
    daily_users=10_000
)
for k, v in plan.items():
    print(f"{k}: {v}")

import hashlib
from collections import Counter

def assign_group(
    user_id: str,
    experiment_id: str,
    split_pct: int = 50
) -> str:
    """
    Deterministic hash-based traffic splitting.

    Properties:
      - Same (user_id, experiment_id) always returns the same group.
      - Different experiment_ids produce independent assignments.
      - Uniform distribution verified empirically on large populations.

    Args:
        user_id: stable server-side user identifier (not cookie)
        experiment_id: unique experiment identifier
        split_pct: percentage of traffic routed to treatment (0-100)

    Returns:
        'treatment' or 'control'
    """
    hash_input = f"{user_id}:{experiment_id}"
    hash_bytes = hashlib.sha256(hash_input.encode('utf-8')).digest()
    # Use first 4 bytes for a 32-bit integer — more than enough entropy
    bucket = int.from_bytes(hash_bytes[:4], 'big') % 100

    return "treatment" if bucket < split_pct else "control"


def validate_split_uniformity(
    experiment_id: str,
    n_users: int = 100_000,
    split_pct: int = 50
) -> dict:
    """
    Empirically verify that hash-based splitting is uniform.

    A non-uniform split means your randomization is biased and
    every experiment result is unreliable. Run this validation
    after any change to the hashing logic.
    """
    counts = Counter(
        assign_group(f"user_{i}", experiment_id, split_pct)
        for i in range(n_users)
    )
    treatment_pct = counts['treatment'] / n_users * 100
    control_pct = counts['control'] / n_users * 100

    # Expected: within 0.5pp of target split for 100K users
    deviation = abs(treatment_pct - split_pct)
    is_uniform = deviation < 1.0  # 1pp tolerance

    return {
        "experiment_id": experiment_id,
        "n_users": n_users,
        "treatment": f"{counts['treatment']} ({treatment_pct:.1f}%)",
        "control": f"{counts['control']} ({control_pct:.1f}%)",
        "deviation_from_target": f"{deviation:.2f}pp",
        "is_uniform": is_uniform
    }


def validate_independence_across_experiments(
    n_users: int = 50_000
) -> dict:
    """
    Verify that assignments across two different experiments are independent.
    A user in treatment for experiment A should have ~50% chance of
    treatment for experiment B.
    """
    both_treatment = 0
    for i in range(n_users):
        uid = f"user_{i}"
        in_A = assign_group(uid, "exp_A") == "treatment"
        in_B = assign_group(uid, "exp_B") == "treatment"
        if in_A and in_B:
            both_treatment += 1

    # Expected: ~25% in both treatment (50% * 50%)
    actual_pct = both_treatment / n_users * 100
    expected_pct = 25.0
    deviation = abs(actual_pct - expected_pct)

    return {
        "both_treatment_pct": f"{actual_pct:.1f}%",
        "expected_pct": f"{expected_pct}%",
        "deviation": f"{deviation:.2f}pp",
        "independent": deviation < 1.0
    }


# Run validations
print("=== Split Uniformity ===")
result = validate_split_uniformity("rec_model_v2")
for k, v in result.items():
    print(f"  {k}: {v}")

print("\n=== Cross-Experiment Independence ===")
result = validate_independence_across_experiments()
for k, v in result.items():
    print(f"  {k}: {v}")

import numpy as np
from scipy import stats

def detect_novelty_effect(
    week1_treatment: float,
    week1_control: float,
    week3_treatment: float,
    week3_control: float,
    threshold: float = 0.50
) -> dict:
    """
    Detect novelty effect by comparing early vs late treatment lift.

    The novelty effect manifests as a positive lift in week 1 that
    decays significantly by week 3. A stable lift across weeks
    indicates genuine improvement; decaying lift indicates curiosity.

    Args:
        week1_treatment: treatment group metric value in week 1
        week1_control: control group metric value in week 1
        week3_treatment: treatment group metric value in week 3
        week3_control: control group metric value in week 3
        threshold: decay fraction above which novelty is flagged (default 0.50)

    Returns:
        dict with detection result, decay percentage, and recommendation
    """
    week1_lift = week1_treatment - week1_control
    week3_lift = week3_treatment - week3_control

    if week1_lift <= 0:
        return {
            "novelty_detected": False,
            "decay_pct": 0.0,
            "week1_lift": round(week1_lift, 4),
            "week3_lift": round(week3_lift, 4),
            "recommendation": "No positive lift in week 1 — novelty not applicable."
        }

    if week3_lift <= 0:
        decay_pct = 100.0
    else:
        decay_pct = (1 - week3_lift / week1_lift) * 100

    novelty_detected = decay_pct > (threshold * 100)

    if novelty_detected:
        recommendation = (
            f"Lift decayed {decay_pct:.0f}% from week 1 to week 3. "
            f"DO NOT SHIP. Extend test to 4+ weeks. "
            f"Segment by new vs returning users. Add post-rollout holdback."
        )
    else:
        recommendation = (
            f"Lift decayed only {decay_pct:.0f}% — appears stable. "
            f"Proceed with caution. Add 5% post-rollout holdback for 2 weeks."
        )

    return {
        "novelty_detected": novelty_detected,
        "decay_pct": round(decay_pct, 1),
        "week1_lift": round(week1_lift, 4),
        "week3_lift": round(week3_lift, 4),
        "recommendation": recommendation
    }


def bonferroni_correction(
    p_values: list[float],
    base_alpha: float = 0.05
) -> list[dict]:
    """
    Apply Bonferroni correction for multiple testing.

    With 20 metrics at alpha=0.05, the family-wise error rate is 64%.
    Bonferroni reduces alpha per metric to keep the overall rate at 5%.
    """
    n = len(p_values)
    corrected_alpha = base_alpha / n

    return [
        {
            "metric_index": i,
            "p_value": round(p, 4),
            "corrected_alpha": round(corrected_alpha, 4),
            "significant_after_correction": p < corrected_alpha
        }
        for i, p in enumerate(p_values)
    ]


# --- Example: detect novelty effect ---
print("=== Novelty Effect Detection ===")
result = detect_novelty_effect(
    week1_treatment=0.085,  # 8.5% CTR in treatment week 1
    week1_control=0.078,    # 7.8% CTR in control week 1
    week3_treatment=0.080,  # 8.0% CTR in treatment week 3
    week3_control=0.079     # 7.9% CTR in control week 3
)
for k, v in result.items():
    print(f"  {k}: {v}")

# --- Example: multiple testing correction ---
print("\n=== Bonferroni Correction ===")
# Simulated p-values from 10 secondary metrics
np.random.seed(42)
p_values = [0.03, 0.12, 0.04, 0.45, 0.72, 0.01, 0.88, 0.06, 0.51, 0.002]
corrected = bonferroni_correction(p_values)
for item in corrected:
    marker = "✓" if item['significant_after_correction'] else "✗"
    print(f"  Metric {item['metric_index']}: p={item['p_value']:.4f} "
          f"corrected_alpha={item['corrected_alpha']:.4f} {marker}")
print(f"\n  Without correction: {sum(1 for p in p_values if p < 0.05)} metrics look significant")
print(f"  With Bonferroni:    {sum(1 for c in corrected if c['significant_after_correction'])} metrics are significant")

package io.thecodeforge.mlops;

import java.nio.charset.StandardCharsets;
import java.security.MessageDigest;
import java.security.NoSuchAlgorithmException;
import java.util.*;
import java.util.concurrent.ConcurrentHashMap;
import java.util.stream.Collectors;

/**
 * Production experiment manager for ML A/B testing.
 *
 * Responsibilities:
 *   - Deterministic hash-based variant assignment
 *   - Thread-safe metric logging with variant tagging
 *   - Automated metric aggregation per variant
 *   - Novelty effect detection across time windows
 *
 * Usage:
 *   1. Create ExperimentManager with experiment ID and split percentage.
 *   2. Call assignVariant(userId) at the start of every request.
 *   3. Call logMetric(userId, metricValue, weekNumber) for every outcome event.
 *   4. Call computeResults() after the pre-committed test duration completes.
 */
public class ExperimentManager {

    private final String experimentId;
    private final int splitPercentage;

    // Per-user, per-week metric storage for novelty detection
    // Key: "variant:userId:week", Value: list of metric observations
    private final ConcurrentHashMap<String, List<Double>> metrics =
        new ConcurrentHashMap<>();

    public ExperimentManager(String experimentId, int splitPercentage) {
        if (splitPercentage < 1 || splitPercentage > 99) {
            throw new IllegalArgumentException(
                "Split percentage must be between 1 and 99, got: " + splitPercentage);
        }
        this.experimentId = experimentId;
        this.splitPercentage = splitPercentage;
    }

    /**
     * Deterministic hash-based variant assignment.
     * Same (userId, experimentId) always returns the same variant.
     * Different experimentIds produce independent assignments.
     */
    public String assignVariant(String userId) {
        try {
            MessageDigest digest = MessageDigest.getInstance("SHA-256");
            String input = userId + ":" + experimentId;
            byte[] hash = digest.digest(input.getBytes(StandardCharsets.UTF_8));
            // Use first 4 bytes for uniform bucket distribution
            int bucket = Math.abs(
                ((hash[0] & 0xFF) << 24) |
                ((hash[1] & 0xFF) << 16) |
                ((hash[2] & 0xFF) << 8)  |
                 (hash[3] & 0xFF)
            ) % 100;
            return bucket < splitPercentage ? "treatment" : "control";
        } catch (NoSuchAlgorithmException e) {
            throw new RuntimeException("SHA-256 unavailable", e);
        }
    }

    /**
     * Log a metric observation tagged with variant, user, and time window.
     * The week parameter enables novelty effect detection by comparing
     * lift in week 1 against lift in week 3.
     */
    public void logMetric(String userId, double metricValue, int week) {
        String variant = assignVariant(userId);
        String key = variant + ":" + userId + ":" + week;
        metrics.computeIfAbsent(key, k -> Collections.synchronizedList(
            new ArrayList<>())).add(metricValue);
    }

    /**
     * Compute per-variant means for a specific week.
     */
    public Map<String, Double> computeWeeklyMeans(int week) {
        double treatmentSum = 0, controlSum = 0;
        int treatmentCount = 0, controlCount = 0;

        for (Map.Entry<String, List<Double>> entry : metrics.entrySet()) {
            String[] parts = entry.getKey().split(":");
            String variant = parts[0];
            int entryWeek = Integer.parseInt(parts[2]);

            if (entryWeek != week) continue;

            double sum = entry.getValue().stream()
                .mapToDouble(Double::doubleValue).sum();
            int count = entry.getValue().size();

            if ("treatment".equals(variant)) {
                treatmentSum += sum;
                treatmentCount += count;
            } else {
                controlSum += sum;
                controlCount += count;
            }
        }

        Map<String, Double> result = new LinkedHashMap<>();
        result.put("treatment_mean",
            treatmentCount > 0 ? treatmentSum / treatmentCount : 0.0);
        result.put("control_mean",
            controlCount > 0 ? controlSum / controlCount : 0.0);
        result.put("lift",
            result.get("treatment_mean") - result.get("control_mean"));
        result.put("treatment_n", (double) treatmentCount);
        result.put("control_n", (double) controlCount);
        return result;
    }

    /**
     * Detect novelty effect by comparing week 1 and week 3 lift.
     */
    public Map<String, Object> detectNovelty() {
        Map<String, Double> week1 = computeWeeklyMeans(1);
        Map<String, Double> week3 = computeWeeklyMeans(3);

        double week1Lift = week1.get("lift");
        double week3Lift = week3.get("lift");

        double decayPct = week1Lift > 0
            ? (1 - week3Lift / week1Lift) * 100
            : 0.0;

        Map<String, Object> result = new LinkedHashMap<>();
        result.put("week1_lift", String.format("%.4f", week1Lift));
        result.put("week3_lift", String.format("%.4f", week3Lift));
        result.put("decay_pct", String.format("%.1f%%", decayPct));
        result.put("novelty_detected", decayPct > 50);
        result.put("recommendation",
            decayPct > 50
                ? "DO NOT SHIP — novelty artifact detected"
                : "Lift appears stable — proceed with holdback");
        return result;
    }

    public static void main(String[] args) {
        ExperimentManager exp = new ExperimentManager("rec_model_v2", 50);
        Random rng = new Random(42);

        // Simulate 3 weeks of data with novelty decay
        for (int week = 1; week <= 3; week++) {
            // Novelty boost decays each week:
            //   week 1: +0.02, week 2: +0.01, week 3: +0.003
            double noveltyBoost = 0.02 / week;

            for (int i = 0; i < 5000; i++) {
                String userId = "user_" + i;
                String variant = exp.assignVariant(userId);
                double baseMetric = 0.05 + rng.nextGaussian() * 0.02;
                double metric = "treatment".equals(variant)
                    ? baseMetric + noveltyBoost
                    : baseMetric;
                exp.logMetric(userId, Math.max(0, metric), week);
            }
        }

        System.out.println("=== Weekly Metric Comparison ===");
        for (int w = 1; w <= 3; w++) {
            Map<String, Double> means = exp.computeWeeklyMeans(w);
            System.out.printf("Week %d: treatment=%.4f control=%.4f lift=%.4f%n",
                w, means.get("treatment_mean"),
                means.get("control_mean"), means.get("lift"));
        }

        System.out.println("\n=== Novelty Effect Detection ===");
        Map<String, Object> novelty = exp.detectNovelty();
        novelty.forEach((k, v) -> System.out.printf("  %s: %s%n", k, v));
    }
}