A/B Testing

graph TD
    A[A/B Testing] --> B[Experiment Design]
    A --> C[Randomization]
    A --> D[Data Collection]
    A --> E[Statistical Analysis]
    A --> F[Decision Making]

    B --> B1[Define hypothesis]
    B --> B2[Select metrics]
    B --> B3[Determine sample size]

    C --> C1[Random assignment]
    C --> C2[Control vs treatment]

    D --> D1[Track user behavior]
    D --> D2[Collect performance data]

    E --> E1[Statistical significance]
    E --> E2[Confidence intervals]
    E --> E3[Effect size]

    F --> F1[Implement winner]
    F --> F2[Iterate and improve]

    style A fill:#f9f,stroke:#333
    style B fill:#cfc,stroke:#333
    style E fill:#fcc,stroke:#333

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

# Simulate A/B test data
np.random.seed(42)
n_samples = 10000
conversion_rate_a = 0.12  # Control group
conversion_rate_b = 0.14  # Treatment group

# Generate data
group_a = np.random.binomial(1, conversion_rate_a, n_samples)
group_b = np.random.binomial(1, conversion_rate_b, n_samples)

# Create DataFrame
data = pd.DataFrame({
    'group': ['A'] * n_samples + ['B'] * n_samples,
    'converted': np.concatenate([group_a, group_b])
})

def analyze_ab_test(data):
    """Analyze A/B test results"""
    # Calculate conversion rates
    conversion_rates = data.groupby('group')['converted'].agg(['mean', 'count', 'sum'])
    conversion_rates.columns = ['conversion_rate', 'sample_size', 'conversions']

    # Calculate difference
    rate_a = conversion_rates.loc['A', 'conversion_rate']
    rate_b = conversion_rates.loc['B', 'conversion_rate']
    difference = rate_b - rate_a
    relative_diff = (difference / rate_a) * 100

    # Statistical test (z-test for proportions)
    successes = conversion_rates['conversions']
    samples = conversion_rates['sample_size']
    z_stat, p_value = proportions_ztest(successes, samples)

    # Confidence intervals
    se_a = np.sqrt(rate_a * (1 - rate_a) / conversion_rates.loc['A', 'sample_size'])
    se_b = np.sqrt(rate_b * (1 - rate_b) / conversion_rates.loc['B', 'sample_size'])
    se_diff = np.sqrt(se_a**2 + se_b**2)

    ci_lower = difference - 1.96 * se_diff
    ci_upper = difference + 1.96 * se_diff

    # Results
    results = {
        'conversion_rates': conversion_rates,
        'difference': difference,
        'relative_difference': relative_diff,
        'p_value': p_value,
        'z_statistic': z_stat,
        'confidence_interval': (ci_lower, ci_upper),
        'significant': p_value < 0.05
    }

    # Print results
    print("A/B Test Results:")
    print(f"Group A: {rate_a:.4f} conversion rate")
    print(f"Group B: {rate_b:.4f} conversion rate")
    print(f"Absolute Difference: {difference:.4f}")
    print(f"Relative Difference: {relative_diff:.2f}%")
    print(f"P-value: {p_value:.4f}")
    print(f"95% Confidence Interval: ({ci_lower:.4f}, {ci_upper:.4f})")
    print(f"Statistically Significant: {'Yes' if results['significant'] else 'No'}")

    # Visualization
    plt.figure(figsize=(10, 6))

    # Bar plot
    plt.subplot(1, 2, 1)
    conversion_rates['conversion_rate'].plot(kind='bar', color=['skyblue', 'salmon'])
    plt.title('Conversion Rates by Group')
    plt.ylabel('Conversion Rate')
    plt.ylim(0, max(rate_a, rate_b) * 1.2)

    # Add value labels
    for i, v in enumerate(conversion_rates['conversion_rate']):
        plt.text(i, v + 0.005, f"{v:.2%}", ha='center')

    # Difference plot
    plt.subplot(1, 2, 2)
    plt.bar(['Difference'], [difference], color='lightgreen')
    plt.axhline(y=0, color='black', linestyle='--')
    plt.title('Conversion Rate Difference')
    plt.ylabel('Difference')
    plt.ylim(min(ci_lower, 0) * 1.2, max(ci_upper, difference) * 1.2)

    # Add confidence interval
    plt.errorbar(['Difference'], [difference], yerr=[[difference - ci_lower], [ci_upper - difference]],
                fmt='o', color='red', capsize=5)

    plt.tight_layout()
    plt.show()

    return results

# Example usage
results = analyze_ab_test(data)

from statsmodels.stats.power import NormalIndPower

def calculate_sample_size(baseline_rate, min_detectable_effect, alpha=0.05, power=0.8):
    """Calculate required sample size for A/B test"""
    # Calculate effect size (Cohen's h)
    h = 2 * (np.arcsin(np.sqrt(baseline_rate + min_detectable_effect)) -
             np.arcsin(np.sqrt(baseline_rate)))

    # Calculate sample size
    analysis = NormalIndPower()
    sample_size = analysis.solve_power(
        effect_size=h,
        nobs1=None,
        alpha=alpha,
        power=power,
        ratio=1.0
    )

    # Round up to nearest integer
    sample_size = int(np.ceil(sample_size))

    print(f"Sample Size Calculation:")
    print(f"Baseline conversion rate: {baseline_rate:.2%}")
    print(f"Minimum detectable effect: {min_detectable_effect:.2%}")
    print(f"Significance level (α): {alpha:.2f}")
    print(f"Power (1-β): {power:.2f}")
    print(f"Required sample size per group: {sample_size:,}")

    return sample_size

# Example usage
sample_size = calculate_sample_size(
    baseline_rate=0.12,
    min_detectable_effect=0.02,  # 2% absolute increase
    alpha=0.05,
    power=0.8
)

class SequentialABTest:
    """Sequential A/B testing with early stopping"""
    def __init__(self, baseline_rate, min_effect=0.02, alpha=0.05, beta=0.2):
        self.baseline_rate = baseline_rate
        self.min_effect = min_effect
        self.alpha = alpha
        self.beta = beta

        # Calculate thresholds
        self.upper_threshold = np.log((1 - self.beta) / self.alpha)
        self.lower_threshold = np.log(self.beta / (1 - self.alpha))

        # Initialize
        self.log_likelihood_ratio = 0
        self.results = {'A': {'conversions': 0, 'samples': 0},
                        'B': {'conversions': 0, 'samples': 0}}

    def update(self, group, converted):
        """Update test with new data"""
        # Update counts
        self.results[group]['conversions'] += converted
        self.results[group]['samples'] += 1

        # Calculate current rates
        rate_a = self.results['A']['conversions'] / self.results['A']['samples'] if self.results['A']['samples'] > 0 else 0
        rate_b = self.results['B']['conversions'] / self.results['B']['samples'] if self.results['B']['samples'] > 0 else 0

        # Calculate log likelihood ratio
        if rate_a > 0 and rate_b > 0 and rate_a < 1 and rate_b < 1:
            llr = (self.results['B']['conversions'] * np.log(rate_b / rate_a) +
                   (self.results['B']['samples'] - self.results['B']['conversions']) * np.log((1 - rate_b) / (1 - rate_a)) +
                   self.results['A']['conversions'] * np.log(rate_a / rate_b) +
                   (self.results['A']['samples'] - self.results['A']['conversions']) * np.log((1 - rate_a) / (1 - rate_b)))

            self.log_likelihood_ratio = llr

        # Check stopping conditions
        if self.log_likelihood_ratio >= self.upper_threshold:
            return 'B_wins'
        elif self.log_likelihood_ratio <= self.lower_threshold:
            return 'A_wins'
        elif self.log_likelihood_ratio == 0 and self.results['A']['samples'] > 1000 and self.results['B']['samples'] > 1000:
            return 'no_difference'
        else:
            return 'continue'

    def get_results(self):
        """Get current test results"""
        rate_a = self.results['A']['conversions'] / self.results['A']['samples'] if self.results['A']['samples'] > 0 else 0
        rate_b = self.results['B']['conversions'] / self.results['B']['samples'] if self.results['B']['samples'] > 0 else 0

        return {
            'conversion_rates': {'A': rate_a, 'B': rate_b},
            'samples': {'A': self.results['A']['samples'], 'B': self.results['B']['samples']},
            'conversions': {'A': self.results['A']['conversions'], 'B': self.results['B']['conversions']},
            'log_likelihood_ratio': self.log_likelihood_ratio,
            'upper_threshold': self.upper_threshold,
            'lower_threshold': self.lower_threshold
        }

    def plot_progress(self):
        """Plot test progress"""
        results = self.get_results()

        plt.figure(figsize=(12, 6))

        # Conversion rates over time
        plt.subplot(1, 2, 1)
        plt.plot(np.cumsum([1] * self.results['A']['samples']), np.cumsum([1] * self.results['A']['conversions']) / np.arange(1, self.results['A']['samples'] + 1),
                label='Group A', color='skyblue')
        plt.plot(np.cumsum([1] * self.results['B']['samples']), np.cumsum([1] * self.results['B']['conversions']) / np.arange(1, self.results['B']['samples'] + 1),
                label='Group B', color='salmon')
        plt.xlabel('Number of Samples')
        plt.ylabel('Conversion Rate')
        plt.title('Conversion Rates Over Time')
        plt.legend()
        plt.grid(True)

        # Log likelihood ratio
        plt.subplot(1, 2, 2)
        plt.plot(np.arange(1, len(self.log_likelihood_ratios) + 1), self.log_likelihood_ratios,
                label='Log Likelihood Ratio', color='green')
        plt.axhline(y=self.upper_threshold, color='red', linestyle='--', label='Upper Threshold')
        plt.axhline(y=self.lower_threshold, color='blue', linestyle='--', label='Lower Threshold')
        plt.xlabel('Number of Samples')
        plt.ylabel('Log Likelihood Ratio')
        plt.title('Sequential Testing Progress')
        plt.legend()
        plt.grid(True)

        plt.tight_layout()
        plt.show()

# Example usage
def simulate_sequential_test(true_rate_a=0.12, true_rate_b=0.14, max_samples=10000):
    """Simulate sequential A/B test"""
    test = SequentialABTest(baseline_rate=true_rate_a, min_effect=0.01)

    # Store log likelihood ratios for plotting
    log_likelihood_ratios = []

    for i in range(max_samples):
        # Randomly assign to group
        group = np.random.choice(['A', 'B'])

        # Simulate conversion
        if group == 'A':
            converted = np.random.binomial(1, true_rate_a)
        else:
            converted = np.random.binomial(1, true_rate_b)

        # Update test
        result = test.update(group, converted)
        current_results = test.get_results()
        log_likelihood_ratios.append(current_results['log_likelihood_ratio'])

        # Check if test should stop
        if result != 'continue':
            print(f"Test stopped after {i+1} samples: {result}")
            print(f"Final results: {current_results}")
            break
    else:
        print(f"Test completed {max_samples} samples without conclusive result")
        print(f"Final results: {current_results}")

    # Store log likelihood ratios for plotting
    test.log_likelihood_ratios = log_likelihood_ratios

    return test

# Run simulation
sequential_test = simulate_sequential_test()
sequential_test.plot_progress()

class EpsilonGreedyBandit:
    """Epsilon-Greedy Multi-Armed Bandit for A/B/n testing"""
    def __init__(self, arms, epsilon=0.1):
        self.arms = arms
        self.epsilon = epsilon
        self.counts = {arm: 0 for arm in arms}
        self.values = {arm: 0.0 for arm in arms}

    def select_arm(self):
        """Select an arm using epsilon-greedy strategy"""
        if np.random.random() < self.epsilon:
            # Explore: choose random arm
            return np.random.choice(self.arms)
        else:
            # Exploit: choose best arm
            return max(self.values.items(), key=lambda x: x[1])[0]

    def update(self, arm, reward):
        """Update arm value with observed reward"""
        self.counts[arm] += 1
        n = self.counts[arm]

        # Update value using incremental average
        self.values[arm] = ((n - 1) / n) * self.values[arm] + (1 / n) * reward

    def get_results(self):
        """Get current bandit results"""
        return {
            'counts': self.counts,
            'values': self.values,
            'best_arm': max(self.values.items(), key=lambda x: x[1])[0]
        }

def simulate_bandit(true_rates, n_trials=10000, epsilon=0.1):
    """Simulate epsilon-greedy bandit"""
    arms = list(true_rates.keys())
    bandit = EpsilonGreedyBandit(arms, epsilon)

    # Store results for analysis
    results = {arm: {'rewards': [], 'counts': []} for arm in arms}

    for _ in range(n_trials):
        # Select arm
        arm = bandit.select_arm()

        # Simulate reward
        reward = np.random.binomial(1, true_rates[arm])

        # Update bandit
        bandit.update(arm, reward)

        # Store results
        for a in arms:
            results[a]['rewards'].append(bandit.values[a])
            results[a]['counts'].append(bandit.counts[a])

    # Plot results
    plt.figure(figsize=(12, 6))

    # Reward estimates over time
    plt.subplot(1, 2, 1)
    for arm in arms:
        plt.plot(results[arm]['rewards'], label=f'Arm {arm}')
    plt.xlabel('Trials')
    plt.ylabel('Estimated Conversion Rate')
    plt.title('Estimated Conversion Rates Over Time')
    plt.legend()
    plt.grid(True)

    # Arm selection counts
    plt.subplot(1, 2, 2)
    counts = [bandit.counts[arm] for arm in arms]
    plt.bar(arms, counts, color=['skyblue', 'salmon', 'lightgreen', 'gold'])
    plt.xlabel('Arms')
    plt.ylabel('Number of Selections')
    plt.title('Arm Selection Counts')
    for i, v in enumerate(counts):
        plt.text(i, v + 50, f"{v:,}", ha='center')

    plt.tight_layout()
    plt.show()

    return bandit

# Example usage
true_rates = {'A': 0.12, 'B': 0.14, 'C': 0.13, 'D': 0.15}
bandit = simulate_bandit(true_rates, n_trials=20000, epsilon=0.1)
print(f"Final results: {bandit.get_results()}")

import pymc3 as pm
import arviz as az

def bayesian_ab_test(data):
    """Bayesian A/B test analysis"""
    # Prepare data
    group_a = data[data['group'] == 'A']['converted'].values
    group_b = data[data['group'] == 'B']['converted'].values

    # Bayesian model
    with pm.Model() as model:
        # Priors
        p_a = pm.Beta('p_a', alpha=1, beta=1)
        p_b = pm.Beta('p_b', alpha=1, beta=1)

        # Likelihood
        obs_a = pm.Bernoulli('obs_a', p=p_a, observed=group_a)
        obs_b = pm.Bernoulli('obs_b', p=p_b, observed=group_b)

        # Difference
        delta = pm.Deterministic('delta', p_b - p_a)

        # Sample
        trace = pm.sample(2000, tune=1000, cores=1)

    # Analyze results
    summary = az.summary(trace)

    # Calculate probability that B is better than A
    prob_b_better = (trace['delta'] > 0).mean()

    # Plot results
    plt.figure(figsize=(15, 10))

    # Posterior distributions
    plt.subplot(2, 2, 1)
    az.plot_posterior(trace, var_names=['p_a', 'p_b'], ref_val=0)
    plt.title('Posterior Distributions')

    # Difference distribution
    plt.subplot(2, 2, 2)
    az.plot_posterior(trace, var_names=['delta'], ref_val=0)
    plt.title('Difference Distribution')

    # Trace plots
    plt.subplot(2, 2, 3)
    az.plot_trace(trace, var_names=['p_a', 'p_b'])
    plt.title('Trace Plots')

    # Forest plot
    plt.subplot(2, 2, 4)
    az.plot_forest(trace, var_names=['p_a', 'p_b', 'delta'])
    plt.title('Forest Plot')

    plt.tight_layout()
    plt.show()

    # Print results
    print("Bayesian A/B Test Results:")
    print(f"Group A conversion rate: {summary.loc['p_a', 'mean']:.4f} (94% HDI: {summary.loc['p_a', 'hdi_3%']:.4f}-{summary.loc['p_a', 'hdi_97%']:.4f})")
    print(f"Group B conversion rate: {summary.loc['p_b', 'mean']:.4f} (94% HDI: {summary.loc['p_b', 'hdi_3%']:.4f}-{summary.loc['p_b', 'hdi_97%']:.4f})")
    print(f"Difference (B - A): {summary.loc['delta', 'mean']:.4f} (94% HDI: {summary.loc['delta', 'hdi_3%']:.4f}-{summary.loc['delta', 'hdi_97%']:.4f})")
    print(f"Probability B is better than A: {prob_b_better:.2%}")

    return {
        'summary': summary,
        'trace': trace,
        'prob_b_better': prob_b_better
    }

# Example usage
bayesian_results = bayesian_ab_test(data)