Probability Interview Questions

# Prior probabilities
p_disease = 0.01
p_no_disease = 0.99

# Likelihood (test characteristics)
p_pos_given_disease = 0.99      # Sensitivity (True Positive Rate)
p_pos_given_no_disease = 0.05   # False Positive Rate (1 - Specificity)

# Evidence: P(Positive) using Law of Total Probability
p_positive = (p_pos_given_disease * p_disease + 
              p_pos_given_no_disease * p_no_disease)
# = 0.99 * 0.01 + 0.05 * 0.99 = 0.0099 + 0.0495 = 0.0594

# Posterior: Apply Bayes' Theorem
p_disease_given_pos = (p_pos_given_disease * p_disease) / p_positive
# = 0.0099 / 0.0594 ≈ 0.167 or 16.7%

print(f"P(Disease|Positive) = {p_disease_given_pos:.1%}")  # 16.7%

# Prior: 30% of emails are spam
p_spam = 0.30
p_not_spam = 0.70

# Likelihood: P("free" appears | spam/not spam)
p_free_given_spam = 0.80       # 80% of spam contains "free"
p_free_given_not_spam = 0.10   # 10% of legitimate emails contain "free"

# Evidence: P("free" appears)
p_free = (p_free_given_spam * p_spam + 
          p_free_given_not_spam * p_not_spam)
# = 0.80 * 0.30 + 0.10 * 0.70 = 0.24 + 0.07 = 0.31

# Posterior: P(spam | "free" appears)
p_spam_given_free = (p_free_given_spam * p_spam) / p_free
# = 0.24 / 0.31 ≈ 0.774 or 77.4%

print(f"P(Spam|'free') = {p_spam_given_free:.1%}")  # 77.4%

# Complete Bayesian inference example
import numpy as np

def bayes_theorem(prior, likelihood, evidence):
    """
    Calculate posterior probability using Bayes' theorem.

    Args:
        prior: P(H) - prior probability of hypothesis
        likelihood: P(E|H) - probability of evidence given hypothesis
        evidence: P(E) - total probability of evidence

    Returns:
        posterior: P(H|E) - updated probability after seeing evidence
    """
    return (likelihood * prior) / evidence

def calculate_evidence(prior, likelihood, likelihood_complement):
    """Calculate P(E) using law of total probability."""
    return likelihood * prior + likelihood_complement * (1 - prior)

# Example: Updated medical test with sequential testing
prior = 0.01  # Initial disease prevalence

# First positive test
sensitivity = 0.99
false_positive_rate = 0.05

evidence = calculate_evidence(prior, sensitivity, false_positive_rate)
posterior_1 = bayes_theorem(prior, sensitivity, evidence)
print(f"After 1st positive test: {posterior_1:.1%}")  # 16.7%

# Second positive test (prior is now the previous posterior)
evidence_2 = calculate_evidence(posterior_1, sensitivity, false_positive_rate)
posterior_2 = bayes_theorem(posterior_1, sensitivity, evidence_2)
print(f"After 2nd positive test: {posterior_2:.1%}")  # 79.5%

# Third positive test
evidence_3 = calculate_evidence(posterior_2, sensitivity, false_positive_rate)
posterior_3 = bayes_theorem(posterior_2, sensitivity, evidence_3)
print(f"After 3rd positive test: {posterior_3:.1%}")  # 98.7%

┌──────────────────────────────────────────────────────────────┐
│          CONDITIONAL PROBABILITY vs INDEPENDENCE              │
├──────────────────────────────────────────────────────────────┤
│                                                              │
│  CONDITIONAL PROBABILITY                                     │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ "What's P(A) if we KNOW B occurred?"                   │ │
│  │                                                        │ │
│  │ Formula: P(A|B) = P(A∩B) / P(B)                       │ │
│  │                                                        │ │
│  │ Example: P(Rain | Dark Clouds) = 0.8                  │ │
│  │          P(Rain alone) = 0.3                          │ │
│  │          → Knowing clouds CHANGES probability         │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  INDEPENDENCE TEST                                           │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Does knowing B change P(A)?                            │ │
│  │                                                        │ │
│  │ IF P(A|B) = P(A) → INDEPENDENT                        │ │
│  │ IF P(A|B) ≠ P(A) → DEPENDENT                          │ │
│  │                                                        │ │
│  │ Equivalent: P(A∩B) = P(A) × P(B) ✓ Independent       │ │
│  └────────────────────────────────────────────────────────┘ │
│                                                              │
└──────────────────────────────────────────────────────────────┘

import numpy as np
import pandas as pd
from typing import Tuple, Dict, List, Optional
from dataclasses import dataclass
from enum import Enum


class EventType(Enum):
    """Event relationship types."""
    INDEPENDENT = "independent"
    DEPENDENT = "dependent"
    MUTUALLY_EXCLUSIVE = "mutually_exclusive"


@dataclass
class ProbabilityAnalysis:
    """Results from conditional probability analysis."""
    p_a: float
    p_b: float
    p_a_and_b: float
    p_a_given_b: float
    p_b_given_a: float
    event_type: EventType
    independence_score: float  # |P(A|B) - P(A)| / P(A)
    details: Dict[str, any]


class ConditionalProbabilityAnalyzer:
    """
    Production-ready analyzer for conditional probability and independence.

    Used by Netflix for user behavior analysis, Google for ad targeting,
    and Meta for feed ranking independence tests.
    """

    def __init__(self, tolerance: float = 1e-6):
        """
        Initialize analyzer.

        Args:
            tolerance: Numerical tolerance for independence test
        """
        self.tolerance = tolerance

    def analyze_events(
        self,
        data: pd.DataFrame,
        event_a: str,
        event_b: str,
        event_a_condition: any = True,
        event_b_condition: any = True
    ) -> ProbabilityAnalysis:
        """
        Analyze conditional probability and independence from data.

        Args:
            data: DataFrame with event columns
            event_a: Column name for event A
            event_b: Column name for event B
            event_a_condition: Value/condition for event A occurring
            event_b_condition: Value/condition for event B occurring

        Returns:
            ProbabilityAnalysis with full statistical breakdown
        """
        n = len(data)

        # Calculate marginal probabilities
        a_occurs = data[event_a] == event_a_condition
        b_occurs = data[event_b] == event_b_condition
        both_occur = a_occurs & b_occurs

        p_a = a_occurs.sum() / n
        p_b = b_occurs.sum() / n
        p_a_and_b = both_occur.sum() / n

        # Calculate conditional probabilities
        p_a_given_b = p_a_and_b / p_b if p_b > 0 else 0
        p_b_given_a = p_a_and_b / p_a if p_a > 0 else 0

        # Determine event relationship
        event_type, independence_score = self._classify_events(
            p_a, p_b, p_a_and_b, p_a_given_b
        )

        details = {
            'n_samples': n,
            'n_a': a_occurs.sum(),
            'n_b': b_occurs.sum(),
            'n_both': both_occur.sum(),
            'expected_if_independent': p_a * p_b * n,
            'observed': both_occur.sum(),
            'deviation': abs(both_occur.sum() - p_a * p_b * n)
        }

        return ProbabilityAnalysis(
            p_a=p_a,
            p_b=p_b,
            p_a_and_b=p_a_and_b,
            p_a_given_b=p_a_given_b,
            p_b_given_a=p_b_given_a,
            event_type=event_type,
            independence_score=independence_score,
            details=details
        )

    def _classify_events(
        self,
        p_a: float,
        p_b: float,
        p_a_and_b: float,
        p_a_given_b: float
    ) -> Tuple[EventType, float]:
        """Classify event relationship."""
        # Check mutual exclusivity
        if abs(p_a_and_b) < self.tolerance:
            return EventType.MUTUALLY_EXCLUSIVE, 1.0

        # Check independence: P(A∩B) ≈ P(A) × P(B)
        expected_if_independent = p_a * p_b
        independence_deviation = abs(p_a_and_b - expected_if_independent)

        if independence_deviation < self.tolerance:
            return EventType.INDEPENDENT, 0.0

        # Calculate independence score (normalized deviation)
        independence_score = abs(p_a_given_b - p_a) / p_a if p_a > 0 else 1.0

        return EventType.DEPENDENT, independence_score

    def chi_square_independence_test(
        self,
        data: pd.DataFrame,
        event_a: str,
        event_b: str
    ) -> Dict[str, float]:
        """
        Perform chi-square test for independence.

        Used by Google Analytics for feature correlation analysis.
        """
        from scipy.stats import chi2_contingency

        # Create contingency table
        contingency = pd.crosstab(data[event_a], data[event_b])

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

        return {
            'chi_square': chi2,
            'p_value': p_value,
            'degrees_of_freedom': dof,
            'is_independent': p_value > 0.05,
            'effect_size_cramers_v': np.sqrt(chi2 / (contingency.sum().sum() * (min(contingency.shape) - 1)))
        }


# ============================================================================
# EXAMPLE 1: NETFLIX - VIDEO STREAMING QUALITY vs USER RETENTION
# ============================================================================

print("=" * 70)
print("EXAMPLE 1: NETFLIX - Streaming Quality Impact on Retention")
print("=" * 70)

# Simulate Netflix user data (10,000 sessions)
np.random.seed(42)
n_users = 10000

# High quality users have 85% retention, low quality 60% retention
quality_high = np.random.choice([True, False], n_users, p=[0.4, 0.6])

retention = []
for is_high_quality in quality_high:
    if is_high_quality:
        retention.append(np.random.choice([True, False], p=[0.85, 0.15]))
    else:
        retention.append(np.random.choice([True, False], p=[0.60, 0.40]))

netflix_data = pd.DataFrame({
    'high_quality': quality_high,
    'retained': retention
})

analyzer = ConditionalProbabilityAnalyzer()
result = analyzer.analyze_events(
    netflix_data,
    event_a='retained',
    event_b='high_quality'
)

print(f"\nP(Retained) = {result.p_a:.3f}")
print(f"P(High Quality) = {result.p_b:.3f}")
print(f"P(Retained ∩ High Quality) = {result.p_a_and_b:.3f}")
print(f"P(Retained | High Quality) = {result.p_a_given_b:.3f}")
print(f"P(Retained | Low Quality) = {(result.p_a - result.p_a_and_b) / (1 - result.p_b):.3f}")
print(f"\nEvent Type: {result.event_type.value.upper()}")
print(f"Independence Score: {result.independence_score:.2%}")
print(f"\n💡 Insight: Knowing quality {'CHANGES' if result.event_type == EventType.DEPENDENT else 'DOES NOT CHANGE'} retention probability")

chi_result = analyzer.chi_square_independence_test(netflix_data, 'retained', 'high_quality')
print(f"\nChi-square test: χ² = {chi_result['chi_square']:.2f}, p = {chi_result['p_value']:.2e}")
print(f"Statistical conclusion: Events are {'INDEPENDENT' if chi_result['is_independent'] else 'DEPENDENT'}")


# ============================================================================
# EXAMPLE 2: GOOGLE ADS - Click-Through Rate Analysis  
# ============================================================================

print("\n" + "=" * 70)
print("EXAMPLE 2: GOOGLE ADS - Ad Position vs CTR Independence")
print("=" * 70)

# Top positions get more clicks (DEPENDENT)
ad_positions = np.random.choice(['top', 'side', 'bottom'], 5000, p=[0.3, 0.5, 0.2])

clicks = []
for pos in ad_positions:
    if pos == 'top':
        clicks.append(np.random.choice([True, False], p=[0.15, 0.85]))
    elif pos == 'side':
        clicks.append(np.random.choice([True, False], p=[0.05, 0.95]))
    else:
        clicks.append(np.random.choice([True, False], p=[0.02, 0.98]))

google_data = pd.DataFrame({
    'position': ad_positions,
    'clicked': clicks
})

# Analyze top position
google_top = google_data.copy()
google_top['is_top'] = google_top['position'] == 'top'

result_google = analyzer.analyze_events(google_top, 'clicked', 'is_top')

print(f"\nP(Click) = {result_google.p_a:.3f}")
print(f"P(Top Position) = {result_google.p_b:.3f}")
print(f"P(Click | Top Position) = {result_google.p_a_given_b:.3f}")
print(f"P(Click | Not Top) = {(result_google.p_a - result_google.p_a_and_b) / (1 - result_google.p_b):.3f}")
print(f"\nLift from Top Position: {(result_google.p_a_given_b / result_google.p_a - 1) * 100:.1f}%")


# ============================================================================
# EXAMPLE 3: TRUE INDEPENDENCE - COIN FLIPS
# ============================================================================

print("\n" + "=" * 70)
print("EXAMPLE 3: COIN FLIPS - True Independence Verification")
print("=" * 70)

# Two independent coin flips
coin1 = np.random.choice([True, False], 10000, p=[0.5, 0.5])
coin2 = np.random.choice([True, False], 10000, p=[0.5, 0.5])

coin_data = pd.DataFrame({
    'flip1_heads': coin1,
    'flip2_heads': coin2
})

result_coin = analyzer.analyze_events(coin_data, 'flip1_heads', 'flip2_heads')

print(f"\nP(Flip1 = Heads) = {result_coin.p_a:.3f}")
print(f"P(Flip2 = Heads) = {result_coin.p_b:.3f}")
print(f"P(Both Heads) = {result_coin.p_a_and_b:.3f}")
print(f"Expected if independent: {result_coin.p_a * result_coin.p_b:.3f}")
print(f"\nEvent Type: {result_coin.event_type.value.upper()}")
print(f"Independence Score: {result_coin.independence_score:.4f} (≈0 confirms independence)")

┌──────────────────────────────────────────────────────────────┐
│         LAW OF TOTAL PROBABILITY WORKFLOW                     │
├──────────────────────────────────────────────────────────────┤
│                                                              │
│  STEP 1: Identify Complex Event A                           │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Event with unknown probability: P(A) = ?               │ │
│  │ Example: P(Customer Churns)                            │ │
│  └────────────┬───────────────────────────────────────────┘ │
│               ↓                                              │
│  STEP 2: Find Partition {B₁, B₂, ..., Bₙ}                  │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Requirements:                                          │ │
│  │  • Mutually exclusive: Bᵢ ∩ Bⱼ = ∅ for i≠j           │ │
│  │  • Exhaustive: ∪Bᵢ = Sample Space                     │ │
│  │                                                        │ │
│  │ Example: {Premium, Standard, Free} user tiers         │ │
│  └────────────┬───────────────────────────────────────────┘ │
│               ↓                                              │
│  STEP 3: Calculate Conditional P(A|Bᵢ) and Prior P(Bᵢ)     │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ P(Churn|Premium) = 0.05,  P(Premium) = 0.20          │ │
│  │ P(Churn|Standard) = 0.15, P(Standard) = 0.50         │ │
│  │ P(Churn|Free) = 0.30,     P(Free) = 0.30             │ │
│  └────────────┬───────────────────────────────────────────┘ │
│               ↓                                              │
│  STEP 4: Apply Formula (Weighted Average)                   │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ P(A) = Σ P(A|Bᵢ) × P(Bᵢ)                             │ │
│  │      = 0.05×0.20 + 0.15×0.50 + 0.30×0.30            │ │
│  │      = 0.01 + 0.075 + 0.09 = 0.175 or 17.5%         │ │
│  └────────────────────────────────────────────────────────┘ │
│                                                              │
│  BONUS: Bayes' Theorem Inversion                            │
│  P(Bᵢ|A) = P(A|Bᵢ) × P(Bᵢ) / P(A)                          │
│                                                              │
└──────────────────────────────────────────────────────────────┘

import numpy as np
import pandas as pd
from typing import List, Dict, Tuple, Callable, Optional
from dataclasses import dataclass
import matplotlib.pyplot as plt


@dataclass
class PartitionElement:
    """Single element in a probability partition."""
    name: str
    prior_prob: float  # P(Bᵢ)
    conditional_prob: float  # P(A|Bᵢ)

    @property
    def contribution(self) -> float:
        """Contribution to total probability."""
        return self.prior_prob * self.conditional_prob


@dataclass
class TotalProbabilityResult:
    """Results from Law of Total Probability calculation."""
    total_probability: float  # P(A)
    partitions: List[PartitionElement]
    posterior_probs: Dict[str, float]  # P(Bᵢ|A) via Bayes
    contributions: Dict[str, float]  # Each partition's contribution
    dominant_partition: str  # Which Bᵢ contributes most


class TotalProbabilityCalculator:
    """
    Production calculator for Law of Total Probability.

    Used by:
    - Amazon: Customer lifetime value across segments
    - Google: Ad click-through rates across devices
    - Uber: Ride acceptance rates across driver tiers
    - Netflix: Content engagement across user cohorts
    """

    def __init__(self):
        pass

    def calculate(
        self,
        partitions: List[PartitionElement],
        validate: bool = True
    ) -> TotalProbabilityResult:
        """
        Apply Law of Total Probability.

        Args:
            partitions: List of partition elements {Bᵢ, P(Bᵢ), P(A|Bᵢ)}
            validate: Check partition validity (exhaustive, mutually exclusive)

        Returns:
            TotalProbabilityResult with P(A) and Bayesian posteriors
        """
        if validate:
            self._validate_partition(partitions)

        # Calculate P(A) = Σ P(A|Bᵢ) × P(Bᵢ)
        total_prob = sum(p.contribution for p in partitions)

        # Calculate contributions
        contributions = {
            p.name: p.contribution for p in partitions
        }

        # Find dominant partition
        dominant = max(partitions, key=lambda p: p.contribution)

        # Calculate posteriors P(Bᵢ|A) via Bayes' theorem
        posteriors = {}
        for p in partitions:
            if total_prob > 0:
                posteriors[p.name] = p.contribution / total_prob
            else:
                posteriors[p.name] = 0.0

        return TotalProbabilityResult(
            total_probability=total_prob,
            partitions=partitions,
            posterior_probs=posteriors,
            contributions=contributions,
            dominant_partition=dominant.name
        )

    def _validate_partition(self, partitions: List[PartitionElement]):
        """Validate partition properties."""
        total_prior = sum(p.prior_prob for p in partitions)

        if not np.isclose(total_prior, 1.0, atol=1e-6):
            raise ValueError(
                f"Partition not exhaustive: Σ P(Bᵢ) = {total_prior:.4f} ≠ 1.0"
            )

        # Check all probabilities in [0,1]
        for p in partitions:
            if not (0 <= p.prior_prob <= 1 and 0 <= p.conditional_prob <= 1):
                raise ValueError(f"Invalid probability for {p.name}")

    def sensitivity_analysis(
        self,
        base_partitions: List[PartitionElement],
        param_name: str,
        param_range: np.ndarray
    ) -> Dict[str, np.ndarray]:
        """
        Analyze sensitivity of P(A) to parameter changes.

        Used by data scientists to understand which partitions drive outcomes.
        """
        results = {'param_values': param_range, 'total_probs': []}

        for param_val in param_range:
            # Create modified partition (simplified: scale one conditional)
            modified = base_partitions.copy()
            # Implementation depends on specific parameter
            # This is a template
            results['total_probs'].append(param_val)  # Placeholder

        return results


# ============================================================================
# EXAMPLE 1: UBER - RIDE ACCEPTANCE RATES ACROSS DRIVER TIERS
# ============================================================================

print("=" * 70)
print("EXAMPLE 1: UBER - Overall Ride Acceptance Rate")
print("=" * 70)

# Uber has 3 driver tiers with different acceptance rates
uber_partitions = [
    PartitionElement(
        name="Diamond (Top 10%)",
        prior_prob=0.10,  # 10% of drivers
        conditional_prob=0.95  # 95% acceptance rate
    ),
    PartitionElement(
        name="Platinum (Next 30%)",
        prior_prob=0.30,
        conditional_prob=0.85  # 85% acceptance rate
    ),
    PartitionElement(
        name="Standard (60%)",
        prior_prob=0.60,
        conditional_prob=0.65  # 65% acceptance rate
    )
]

calc = TotalProbabilityCalculator()
uber_result = calc.calculate(uber_partitions)

print(f"\nOverall acceptance rate: {uber_result.total_probability:.1%}")
print(f"\nContributions by tier:")
for name, contrib in uber_result.contributions.items():
    pct_of_total = contrib / uber_result.total_probability * 100
    print(f"  {name}: {contrib:.3f} ({pct_of_total:.1f}% of total)")

print(f"\nDominant tier: {uber_result.dominant_partition}")

print(f"\nIf ride accepted, which tier? (Bayesian posterior):")
for name, post_prob in uber_result.posterior_probs.items():
    print(f"  P({name} | Accepted) = {post_prob:.1%}")

print(f"\n💡 Insight: Standard tier drivers (60%) contribute ")
print(f"   {uber_result.contributions['Standard (60%)'] / uber_result.total_probability:.1%} of acceptances")


# ============================================================================
# EXAMPLE 2: AMAZON - PRODUCT DEFECT RATE ACROSS SUPPLIERS
# ============================================================================

print("\n" + "=" * 70)
print("EXAMPLE 2: AMAZON - Product Defect Rate Analysis")
print("=" * 70)

# Three suppliers with different defect rates
amazon_partitions = [
    PartitionElement("Supplier A (High Volume)", 0.55, 0.02),
    PartitionElement("Supplier B (Medium Volume)", 0.30, 0.035),
    PartitionElement("Supplier C (Low Volume)", 0.15, 0.08)
]

amazon_result = calc.calculate(amazon_partitions)

print(f"\nOverall defect rate: {amazon_result.total_probability:.2%}")
print(f"Expected defective units per 10,000: {amazon_result.total_probability * 10000:.0f}")

print(f"\nIf product is defective, which supplier?")
for name, post_prob in amazon_result.posterior_probs.items():
    print(f"  {name}: {post_prob:.1%}")

print(f"\n🎯 Action: Focus QA on {amazon_result.dominant_partition}")


# ============================================================================
# EXAMPLE 3: GOOGLE ADS - CLICK-THROUGH RATE ACROSS DEVICES
# ============================================================================

print("\n" + "=" * 70)
print("EXAMPLE 3: GOOGLE ADS - Overall CTR Across Devices")
print("=" * 70)

google_partitions = [
    PartitionElement("Desktop", 0.35, 0.08),  # 35% traffic, 8% CTR
    PartitionElement("Mobile", 0.55, 0.04),   # 55% traffic, 4% CTR
    PartitionElement("Tablet", 0.10, 0.06)    # 10% traffic, 6% CTR
]

google_result = calc.calculate(google_partitions)

print(f"\nOverall CTR: {google_result.total_probability:.2%}")
print(f"Revenue if $2 per click on 1M impressions: ${google_result.total_probability * 1e6 * 2:,.0f}")

print(f"\nDevice mix optimization:")
for p in google_partitions:
    current_revenue = p.contribution * 1e6 * 2
    print(f"  {p.name}: ${current_revenue:,.0f} ({p.prior_prob:.0%} traffic × {p.conditional_prob:.1%} CTR)")

# What if we increase mobile CTR by 1%?
google_partitions_improved = [
    PartitionElement("Desktop", 0.35, 0.08),
    PartitionElement("Mobile", 0.55, 0.05),  # 4% → 5%
    PartitionElement("Tablet", 0.10, 0.06)
]
google_result_improved = calc.calculate(google_partitions_improved)

revenue_lift = (google_result_improved.total_probability - google_result.total_probability) * 1e6 * 2
print(f"\n🚀 If mobile CTR improves 4% → 5%: +${revenue_lift:,.0f} revenue")

┌──────────────────────────────────────────────────────────────┐
│        EXPECTED VALUE PROPERTIES HIERARCHY                    │
├──────────────────────────────────────────────────────────────┤
│                                                              │
│  PROPERTY 1: LINEARITY (Most Important!)                     │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ E[aX + bY + c] = aE[X] + bE[Y] + c                     │ │
│  │                                                        │ │
│  │ ✅ Works for ANY X, Y (even dependent!)                │ │
│  │ ✅ Extends to any linear combination                   │ │
│  │ ✅ Foundation of portfolio theory, ML loss functions   │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PROPERTY 2: PRODUCT (Requires Independence)                 │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ E[XY] = E[X] · E[Y]   IFF X ⊥ Y                       │ │
│  │                                                        │ │
│  │ ⚠️  Only if X and Y are independent                    │ │
│  │ ⚠️  Otherwise: E[XY] = E[X]E[Y] + Cov(X,Y)            │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PROPERTY 3: MONOTONICITY                                    │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ If X ≤ Y, then E[X] ≤ E[Y]                            │ │
│  │                                                        │ │
│  │ Preserves ordering                                     │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PROPERTY 4: LAW OF ITERATED EXPECTATIONS                    │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ E[X] = E[E[X|Y]]                                       │ │
│  │                                                        │ │
│  │ "Expectation of conditional expectation = expectation"│ │
│  │ Used in hierarchical models, Bayesian inference       │ │
│  └────────────────────────────────────────────────────────┘ │
│                                                              │
└──────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────┐
│         VARIANCE PROPERTIES (CRITICAL FOR INTERVIEWS)         │
├──────────────────────────────────────────────────────────────┤
│                                                              │
│  PROPERTY 1: Variance of Constant = 0                        │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Var(c) = 0                                             │ │
│  │                                                        │ │
│  │ No randomness → no variance                           │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PROPERTY 2: Scaling (QUADRATIC!)                            │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Var(aX) = a² · Var(X)                                 │ │
│  │                                                        │ │
│  │ ⚠️  SQUARES the constant (unlike expectation)          │ │
│  │ Example: Var(2X) = 4·Var(X), not 2·Var(X)            │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PROPERTY 3: Translation Invariance                          │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Var(X + b) = Var(X)                                    │ │
│  │                                                        │ │
│  │ Shifting all values doesn't change spread             │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PROPERTY 4: Sum of Independent Variables                    │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Var(X + Y) = Var(X) + Var(Y)   IFF X ⊥ Y             │ │
│  │                                                        │ │
│  │ Variances ADD for independent variables               │ │
│  │ ⚠️  Requires independence (unlike expectation)         │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PROPERTY 5: Sum with Covariance (General Case)              │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Var(X + Y) = Var(X) + Var(Y) + 2·Cov(X,Y)            │ │
│  │                                                        │ │
│  │ Covariance term captures dependence                   │ │
│  │ Cov(X,Y) > 0 → more variance (positive correlation)   │ │
│  │ Cov(X,Y) < 0 → less variance (hedging effect)         │ │
│  └────────────────────────────────────────────────────────┘ │
│                                                              │
└──────────────────────────────────────────────────────────────┘

import numpy as np

# Roll of fair die
outcomes = np.array([1, 2, 3, 4, 5, 6])
probs = np.array([1/6] * 6)

# E[X] = 3.5
E_X = np.sum(outcomes * probs)

# E[X²] = 1²·(1/6) + 2²·(1/6) + ... + 6²·(1/6) = 91/6 ≈ 15.167
E_X2 = np.sum(outcomes**2 * probs)

# Var(X) = E[X²] - (E[X])²  = 91/6 - (7/2)² = 91/6 - 49/4 ≈ 2.917
variance = E_X2 - E_X**2
std_dev = np.sqrt(variance)  # σ ≈ 1.708

print(f"E[X] = {E_X:.3f}")
print(f"Var(X) = {variance:.3f}")
print(f"SD(X) = {std_dev:.3f}")

# Scaling demonstration: Var(2X) = 4·Var(X)
variance_2x = np.var(2 * outcomes * np.repeat(1/6, 6))
print(f"\nVar(2X) = {4 * variance:.3f}  (4 times larger!)")

┌──────────────────────────────────────────────────────────────┐
│         CENTRAL LIMIT THEOREM MAGIC EXPLAINED                 │
├──────────────────────────────────────────────────────────────┤
│                                                              │
│  STEP 1: Start with ANY Population Distribution              │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Exponential, Uniform, Binomial, Even Bimodal!          │ │
│  │                                                        │ │
│  │ Only requirement: Finite variance σ²                  │ │
│  │ Population: μ (mean), σ² (variance)                   │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  STEP 2: Draw Samples of Size n                              │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Take n observations: X₁, X₂, ..., Xₙ                  │ │
│  │                                                        │ │
│  │ Calculate sample mean: X̄ = (ΣXᵢ) / n                 │ │
│  │                                                        │ │
│  │ Repeat many times → get distribution of X̄            │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  STEP 3: Observe the Miracle! 🎉                            │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Distribution of X̄ becomes NORMAL as n increases!     │ │
│  │                                                        │ │
│  │ • Mean: E[X̄] = μ (same as population)                │ │
│  │ • Variance: Var(X̄) = σ²/n (decreases with n)        │ │
│  │ • Std Dev (SE): SD(X̄) = σ/√n ("Standard Error")     │ │
│  │                                                        │ │
│  │ For n≥30: X̄ ~ N(μ, σ²/n) approximately              │ │
│  └────────────────────────────────────────────────────────┘ │
│               ↓                                              │
│  PRACTICAL CONSEQUENCE                                       │
│  ┌────────────────────────────────────────────────────────┐ │
│  │ Can use z-tests, t-tests, confidence intervals         │ │
│  │ WITHOUT assuming population is normal!                 │ │
│  │                                                        │ │
│  │ This is the foundation of A/B testing!                │ │
│  └────────────────────────────────────────────────────────┘ │
│                                                              │
└──────────────────────────────────────────────────────────────┘

import numpy as np
import pandas as pd
from typing import Callable, List, Tuple
from scipy import stats
import matplotlib.pyplot as plt
from dataclasses import dataclass


@dataclass
class CLTDemonstration:
    """Results from CLT simulation."""
    population_mean: float
    population_std: float
    sample_size: int
    num_samples: int
    sample_means: np.ndarray
    theoretical_mean: float
    theoretical_std_error: float
    empirical_mean: float
    empirical_std_error: float
    normality_test_pvalue: float


class CentralLimitTheoremAnalyzer:
    """
    Demonstrate and apply Central Limit Theorem.

    Used by:
    - Google: A/B test sample size calculations
    - Netflix: Confidence intervals for engagement metrics
    - Amazon: Quality control (defect rate estimation)
    - Uber: Trip duration confidence intervals
    """

    def demonstrate_clt(
        self,
        population_dist: Callable,
        sample_size: int,
        num_samples: int = 10000,
        population_mean: Optional[float] = None,
        population_std: Optional[float] = None
    ) -> CLTDemonstration:
        """
        Demonstrate CLT by simulation.

        Args:
            population_dist: Function that generates n samples from population
            sample_size: Size of each sample (n)
            num_samples: Number of sample means to generate
            population_mean: True population mean (if known)
            population_std: True population std dev (if known)

        Returns:
            CLTDemonstration with empirical and theoretical statistics
        """
        # Generate many sample means
        sample_means = np.array([
            np.mean(population_dist(sample_size))
            for _ in range(num_samples)
        ])

        # Empirical statistics from simulation
        empirical_mean = np.mean(sample_means)
        empirical_std_error = np.std(sample_means, ddof=1)

        # Theoretical statistics (if population params known)
        if population_mean is None or population_std is None:
            # Estimate from large sample
            large_sample = population_dist(100000)
            population_mean = np.mean(large_sample)
            population_std = np.std(large_sample, ddof=1)

        theoretical_std_error = population_std / np.sqrt(sample_size)

        # Test normality (Shapiro-Wilk)
        _, normality_pvalue = stats.shapiro(
            sample_means[:5000]  # Shapiro-Wilk limit
        )

        return CLTDemonstration(
            population_mean=population_mean,
            population_std=population_std,
            sample_size=sample_size,
            num_samples=num_samples,
            sample_means=sample_means,
            theoretical_mean=population_mean,
            theoretical_std_error=theoretical_std_error,
            empirical_mean=empirical_mean,
            empirical_std_error=empirical_std_error,
            normality_test_pvalue=normality_pvalue
        )

    def minimum_sample_size(
        self,
        population_std: float,
        margin_of_error: float,
        confidence_level: float = 0.95
    ) -> int:
        """
        Calculate minimum sample size for desired precision.

        Based on CLT: n = (z*σ / E)²
        where E = margin of error

        Used by data scientists for experiment design.
        """
        z_score = stats.norm.ppf((1 + confidence_level) / 2)
        n = (z_score * population_std / margin_of_error) ** 2
        return int(np.ceil(n))


# ============================================================================
# EXAMPLE 1: EXPONENTIAL DISTRIBUTION → NORMAL SAMPLE MEANS
# ============================================================================

print("=" * 70)
print("EXAMPLE 1: CLT with Exponential Distribution (Highly Skewed)")
print("=" * 70)

analyzer = CentralLimitTheoremAnalyzer()

# Exponential(λ=1): Mean=1, Var=1, Highly right-skewed
exponential_dist = lambda n: np.random.exponential(scale=1.0, size=n)

# Small sample size (n=5) - CLT weak
result_n5 = analyzer.demonstrate_clt(
    exponential_dist,
    sample_size=5,
    population_mean=1.0,
    population_std=1.0
)

print(f"\nSample size n=5:")
print(f"  Theoretical SE: {result_n5.theoretical_std_error:.3f}")
print(f"  Empirical SE: {result_n5.empirical_std_error:.3f}")
print(f"  Normality test p-value: {result_n5.normality_test_pvalue:.4f}")
print(f"  Normal? {result_n5.normality_test_pvalue > 0.05}")

# Medium sample size (n=30) - CLT kicks in!
result_n30 = analyzer.demonstrate_clt(
    exponential_dist,
    sample_size=30,
    population_mean=1.0,
    population_std=1.0
)

print(f"\nSample size n=30:")
print(f"  Theoretical SE: {result_n30.theoretical_std_error:.3f}")
print(f"  Empirical SE: {result_n30.empirical_std_error:.3f}")
print(f"  Normality test p-value: {result_n30.normality_test_pvalue:.4f}")
print(f"  Normal? {result_n30.normality_test_pvalue > 0.05}")

# Large sample size (n=100) - Strongly normal
result_n100 = analyzer.demonstrate_clt(
    exponential_dist,
    sample_size=100,
    population_mean=1.0,
    population_std=1.0
)

print(f"\nSample size n=100:")
print(f"  Theoretical SE: {result_n100.theoretical_std_error:.3f}")
print(f"  Empirical SE: {result_n100.empirical_std_error:.3f}")
print(f"  Normality test p-value: {result_n100.normality_test_pvalue:.4f}")
print(f"  Normal? {result_n100.normality_test_pvalue > 0.05}")

print(f"\n🎯 As n increases, SE decreases (√n): {1/np.sqrt(5):.3f} → {1/np.sqrt(30):.3f} → {1/np.sqrt(100):.3f}")


# ============================================================================
# EXAMPLE 2: NETFLIX - CONFIDENCE INTERVAL FOR AVG WATCH TIME
# ============================================================================

print("\n" + "=" * 70)
print("EXAMPLE 2: NETFLIX - Watch Time Confidence Interval (CLT)")
print("=" * 70)

# Sample of 500 users
# Population: Unknown distribution (probably right-skewed)
# But CLT lets us use normal inference!

np.random.seed(42)
watch_times = np.random.gamma(shape=2, scale=2.5, size=500)  # Skewed data

n = len(watch_times)
sample_mean = np.mean(watch_times)
sample_std = np.std(watch_times, ddof=1)

# Standard error (by CLT)
se = sample_std / np.sqrt(n)

# 95% confidence interval
z_95 = 1.96
ci_95 = (sample_mean - z_95 * se, sample_mean + z_95 * se)

print(f"\nSample size: {n} users")
print(f"Sample mean: {sample_mean:.2f} hours")
print(f"Sample std dev: {sample_std:.2f} hours")
print(f"Standard error: {se:.3f} hours")
print(f"\n95% CI: ({ci_95[0]:.2f}, {ci_95[1]:.2f}) hours")
print(f"\n🎬 We're 95% confident true mean watch time is in [{ci_95[0]:.2f}, {ci_95[1]:.2f}]")
print(f"   (Thanks to CLT, even though data is skewed!)")


# ============================================================================
# EXAMPLE 3: GOOGLE - A/B TEST SAMPLE SIZE CALCULATION
# ============================================================================

print("\n" + "=" * 70)
print("EXAMPLE 3: GOOGLE - Sample Size for A/B Test (CLT-based)")
print("=" * 70)

# Google wants to detect 2% CTR improvement
# Control CTR: 5% (σ ≈ √(0.05 × 0.95) ≈ 0.218 for binary outcome)

baseline_ctr = 0.05
population_std = np.sqrt(baseline_ctr * (1 - baseline_ctr))

# Want margin of error = 0.005 (0.5%) at 95% confidence
margin_of_error = 0.005

n_required = analyzer.minimum_sample_size(
    population_std=population_std,
    margin_of_error=margin_of_error,
    confidence_level=0.95
)

print(f"\nBaseline CTR: {baseline_ctr:.1%}")
print(f"Population std: {population_std:.3f}")
print(f"Desired margin of error: {margin_of_error:.2%}")
print(f"\nRequired sample size per variant: {n_required:,}")
print(f"Total experiment size: {2 * n_required:,}")

# At 1M daily users, how long to run?
daily_users = 1_000_000
users_per_variant = daily_users / 2
days_needed = np.ceil(n_required / users_per_variant)

print(f"\n📊 With {daily_users:,} daily users (50/50 split):")
print(f"   Need to run experiment for {int(days_needed)} days")

┌────────────────────────────────────────────────────────────┐
│         BINOMIAL DISTRIBUTION CONDITIONS (BINS)            │
├────────────────────────────────────────────────────────────┤
│                                                            │
│  B - Binary outcomes:    Each trial has 2 outcomes        │
│                           (Success/Failure, Yes/No)        │
│                                                            │
│  I - Independent trials: Outcome of one trial doesn't     │
│                           affect others                    │
│                                                            │
│  N - Number fixed:       n trials determined in advance   │
│                           (not stopping after X successes) │
│                                                            │
│  S - Same probability:   p constant across all trials     │
│                           (homogeneous)                    │
│                                                            │
│  \ud83d\udcca Distribution Shape:                                    │
│                                                            │
│   p=0.1 (skewed right)     p=0.5 (symmetric)              │
│   █                        ████                            │
│   ██                      ██████                           │
│   ███                    ████████                          │
│   ████                  ██████████                         │
│   ─────────────────     ────────────────                  │
│   0 1 2 3 4 5...        0  1  2  3  4  5                  │
│                                                            │
└────────────────────────────────────────────────────────────┘

from scipy.stats import poisson

# 4 customers per hour on average
lambda_rate = 4

# P(exactly 6 customers)?
p_6 = poisson.pmf(k=6, mu=4)  # ≈ 0.104

# P(at most 2 customers)?
p_le_2 = poisson.cdf(k=2, mu=4)  # ≈ 0.238

# P(more than 5)?
p_gt_5 = 1 - poisson.cdf(k=5, mu=4)  # ≈ 0.215

from scipy.stats import expon

# λ = 0.1, scale = 1/λ = 10
wait_time = expon(scale=10)

# P(wait < 5 minutes)?
p_lt_5 = wait_time.cdf(5)  # ≈ 0.393

# P(wait > 15 minutes)?
p_gt_15 = 1 - wait_time.cdf(15)  # ≈ 0.223

# Mean wait time
mean_wait = wait_time.mean()  # 10 minutes

from scipy.stats import geom

p = 0.3

# P(pass on exactly 3rd interview)?
p_3rd = geom.pmf(k=3, p=0.3)
# = (0.7)^2 * 0.3 = 0.147

# Expected interviews until first pass?
expected = 1 / 0.3  # ≈ 3.33 interviews

# P(need more than 5 interviews)?
p_gt_5 = 1 - geom.cdf(k=5, p=0.3)
# = (0.7)^5 ≈ 0.168

def birthday_probability(n):
    """P(at least 2 share birthday in group of n)"""
    p_all_different = 1.0
    for i in range(n):
        p_all_different *= (365 - i) / 365
    return 1 - p_all_different

# Results:
# n=23: 50.7% (famous result!)
# n=50: 97.0%
# n=70: 99.9%

for n in [10, 23, 30, 50, 70]:
    print(f"n={n}: {birthday_probability(n):.1%}")

Initial pick: P(Car) = 1/3
Other doors:  P(Car) = 2/3

After host reveals goat:
- Your door still has P = 1/3
- Remaining door gets all 2/3

import random

def monty_hall(switch, n_simulations=100000):
    wins = 0
    for _ in range(n_simulations):
        car = random.randint(0, 2)
        choice = random.randint(0, 2)

        # Host opens a goat door (not your choice, not car)
        goat_doors = [i for i in range(3) if i != choice and i != car]
        host_opens = random.choice(goat_doors)

        if switch:
            # Switch to remaining door
            choice = [i for i in range(3) if i != choice and i != host_opens][0]

        if choice == car:
            wins += 1

    return wins / n_simulations

print(f"Stay:   {monty_hall(switch=False):.1%}")   # ~33.3%
print(f"Switch: {monty_hall(switch=True):.1%}")    # ~66.7%

# Covariance
Cov(X, X) = Var(X)
Cov(X, Y) = Cov(Y, X)  # Symmetric
Cov(aX + b, Y) = a·Cov(X, Y)

# Correlation
Corr(aX + b, Y) = sign(a) · Corr(X, Y)  # Unaffected by linear transform

import numpy as np

x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 4, 5, 4, 5])

cov_matrix = np.cov(x, y)
cov_xy = cov_matrix[0, 1]  # Covariance

corr_matrix = np.corrcoef(x, y)
corr_xy = corr_matrix[0, 1]  # Correlation

import numpy as np
import matplotlib.pyplot as plt

# Fair coin: P(Heads) = 0.5
np.random.seed(42)

flips = np.random.binomial(1, 0.5, 10000)
running_mean = np.cumsum(flips) / np.arange(1, 10001)

plt.plot(running_mean)
plt.axhline(y=0.5, color='r', linestyle='--', label='True Mean')
plt.xlabel('Number of Flips')
plt.ylabel('Running Mean')
plt.title('Law of Large Numbers: Coin Flips')

from scipy import stats
import numpy as np

# Discrete: Binomial PMF and CDF
x_discrete = np.arange(0, 11)
pmf = stats.binom.pmf(x_discrete, n=10, p=0.5)
cdf_discrete = stats.binom.cdf(x_discrete, n=10, p=0.5)

# Continuous: Normal PDF and CDF
x_continuous = np.linspace(-4, 4, 100)
pdf = stats.norm.pdf(x_continuous)
cdf_continuous = stats.norm.cdf(x_continuous)

from scipy import stats
import numpy as np

data = [23, 25, 28, 22, 26, 27, 24, 29, 25, 26]

# 95% CI for mean
mean = np.mean(data)
se = stats.sem(data)  # Standard error
ci = stats.t.interval(0.95, len(data)-1, loc=mean, scale=se)

print(f"95% CI: ({ci[0]:.2f}, {ci[1]:.2f})")

from scipy import stats

# Control: 100 conversions out of 1000
# Treatment: 120 conversions out of 1000

control_conv = 100
control_n = 1000
treatment_conv = 120
treatment_n = 1000

# H₀: p1 = p2 (no difference)
# H₁: p1 ≠ p2 (difference exists)

# Two-proportion z-test
from statsmodels.stats.proportion import proportions_ztest

stat, pvalue = proportions_ztest(
    [control_conv, treatment_conv],
    [control_n, treatment_n]
)

print(f"p-value: {pvalue:.4f}")
# If p < 0.05, reject H₀ → significant difference

from statsmodels.stats.power import TTestIndPower

# Parameters
effect_size = 0.5  # Cohen's d (medium effect)
alpha = 0.05
power = 0.80

# Calculate required sample size
analysis = TTestIndPower()
n = analysis.solve_power(
    effect_size=effect_size,
    alpha=alpha,
    power=power,
    ratio=1  # Equal group sizes
)

print(f"Required n per group: {n:.0f}")
# ~64 per group for medium effect

from math import factorial, comb, perm

# 5 people, select 3 for positions (President, VP, Secretary)
# Order matters → Permutation
positions = perm(5, 3)  # = 5 × 4 × 3 = 60

# 5 people, select 3 for a committee
# Order doesn't matter → Combination
committee = comb(5, 3)  # = 10

# Relationship
assert perm(5, 3) == comb(5, 3) * factorial(3)
# 60 = 10 × 6

# 4-digit PIN (0-9): 10^4 = 10000
# Choose 3 scoops from 5 flavors (repeats OK): C(5+3-1, 3) = C(7,3) = 35

from scipy.stats import nbinom

r, p = 3, 0.8

# Expected trials until 3 successes
expected_trials = r / p  # = 3 / 0.8 = 3.75

# P(need exactly 5 trials)?
# 5 trials, 3 successes, 2 failures
p_5 = nbinom.pmf(k=2, n=3, p=0.8)  # k = failures
# = C(4,2) * 0.8^3 * 0.2^2 = 0.0512

# P(need at most 4 trials)?
p_le_4 = nbinom.cdf(k=1, n=3, p=0.8)  # ≤1 failure