PAC Learning: Probably Approximately Correct Framework

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

import numpy as np
import matplotlib.pyplot as plt
from typing import Callable
 
def demonstrate_generalization_gap():
    """
    Empirically demonstrate how training error converges to true error
    as sample size increases.
    """
    np.random.seed(42)
    
    # True concept: AND of first 3 variables (x0 ∧ x1 ∧ x2)
    def true_concept(x):
        return int(x[0] == 1 and x[1] == 1 and x[2] == 1)
    
    # Hypothesis: Approximation (AND of first 2 variables)
    def hypothesis(x):
        return int(x[0] == 1 and x[1] == 1)  # Missing x2
    
    # True distribution: Uniform over {0,1}^4
    def sample_from_D():
        return np.random.randint(0, 2, 4)
    
    # True error of hypothesis
    # h(x)=1 when x0=x1=1, c*(x)=1 when x0=x1=x2=1
    # Error when h(x) != c*(x):
    #   - h=1, c*=0: x0=1, x1=1, x2=0 (prob = 1/8)
    # True error = 1/8 = 0.125
    TRUE_ERROR = 0.125
    
    # Measure training error for various sample sizes
    sample_sizes = [10, 25, 50, 100, 250, 500, 1000, 2500, 5000]
    n_trials = 100
    
    results = {m: [] for m in sample_sizes}
    
    for m in sample_sizes:
        for _ in range(n_trials):
            # Generate training set
            training_errors = 0
            for _ in range(m):
                x = sample_from_D()
                if hypothesis(x) != true_concept(x):
                    training_errors += 1
            
            training_error = training_errors / m
            gap = abs(training_error - TRUE_ERROR)
            results[m].append(gap)
    
    # Compute statistics
    mean_gaps = [np.mean(results[m]) for m in sample_sizes]
    std_gaps = [np.std(results[m]) for m in sample_sizes]
    
    # Hoeffding bound for comparison
    hoeffding_bounds = [np.sqrt(np.log(2) / (2 * m)) for m in sample_sizes]
    
    print("Generalization Gap Analysis:")
    print("-" * 60)
    for m, mean, std, bound in zip(sample_sizes, mean_gaps, std_gaps, 
                                   hoeffding_bounds):
        print(f"m={m:5d}: Mean gap={mean:.4f} ± {std:.4f}, "
              f"Hoeffding bound={bound:.4f}")
    
    # Key observations:
    # 1. Gap decreases as m increases (as expected)
    # 2. Empirical gap is typically smaller than Hoeffding bound
    # 3. Gap scales roughly as 1/sqrt(m)
    
    return results
 
# Example output:
# m=   10: Mean gap=0.0875 ± 0.0632, Hoeffding bound=0.1863
# m=  100: Mean gap=0.0275 ± 0.0188, Hoeffding bound=0.0589
# m= 1000: Mean gap=0.0092 ± 0.0058, Hoeffding bound=0.0186
# m= 5000: Mean gap=0.0037 ± 0.0024, Hoeffding bound=0.0083

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

import numpy as np
from math import log, ceil
 
def sample_complexity_finite_class(epsilon: float, delta: float, 
                                    hypothesis_class_size: int) -> int:
    """
    Compute sample complexity for PAC learning a finite hypothesis class.
    
    Args:
        epsilon: Desired accuracy (max generalization error)
        delta: Failure probability
        hypothesis_class_size: |H|, number of hypotheses
    
    Returns:
        m: Number of samples sufficient for PAC learning
    
    Derivation:
        From union bound + Hoeffding:
        Pr[∃h: |err_S(h) - err_D(h)| > ε] ≤ 2|H|exp(-2mε²)
        
        Setting this ≤ δ and solving:
        m ≥ (1/2ε²) ln(2|H|/δ)
    """
    return ceil((1 / (2 * epsilon**2)) * log(2 * hypothesis_class_size / delta))
 
 
def analyze_sample_complexity():
    """Analyze how parameters affect sample complexity."""
    
    print("Sample Complexity Analysis for Finite Hypothesis Classes")
    print("=" * 70)
    
    # Fix δ, vary ε
    print("
1. Effect of ε (accuracy) [|H|=1000, δ=0.05]:")
    for epsilon in [0.20, 0.10, 0.05, 0.02, 0.01]:
        m = sample_complexity_finite_class(epsilon, 0.05, 1000)
        print(f"   ε={epsilon}: m={m:,}")
    
    # Fix ε, vary |H|
    print("
2. Effect of |H| (hypothesis class size) [ε=0.05, δ=0.05]:")
    for H_size in [10, 100, 1000, 10000, 100000, 1000000]:
        m = sample_complexity_finite_class(0.05, 0.05, H_size)
        print(f"   |H|={H_size:>10,}: m={m:,}")
    
    # Fix ε, vary δ
    print("
3. Effect of δ (failure probability) [|H|=1000, ε=0.05]:")
    for delta in [0.20, 0.10, 0.05, 0.01, 0.001]:
        m = sample_complexity_finite_class(0.05, delta, 1000)
        print(f"   δ={delta}: m={m:,}")
    
    # Concrete examples
    print("
4. Concrete Learning Scenarios:")
    
    # Boolean conjunctions: |H| = 3^n
    n = 10
    H_conj = 3 ** n  # Each of n variables: included positive, negative, or absent
    m = sample_complexity_finite_class(0.05, 0.05, H_conj)
    print(f"   Conjunctions (n={n}, |H|={H_conj:,}): m={m:,}")
    
    # Decision stumps: |H| = 2n  
    n = 100
    H_stump = 2 * n  # Each of n features, threshold above or below
    m = sample_complexity_finite_class(0.05, 0.05, H_stump)
    print(f"   Decision stumps (n={n}, |H|={H_stump}): m={m:,}")
    
    # Decision trees depth d: |H| ≈ (2n)^(2^d)
    n, d = 10, 3
    H_tree = (2 * n) ** (2 ** d)
    m = sample_complexity_finite_class(0.05, 0.05, H_tree)
    print(f"   Decision trees (n={n}, d={d}, |H|≈{H_tree:,}): m={m:,}")
 
analyze_sample_complexity()
 
# Sample Output:
# 1. Effect of ε (accuracy):
#    ε=0.20: m=204
#    ε=0.10: m=817
#    ε=0.05: m=3,268
#    ε=0.02: m=20,422
#    ε=0.01: m=81,687
#
# 2. Effect of |H| (logarithmic growth!):
#    |H|=    10: m=1,838
#    |H|=   100: m=2,298
#    |H|= 1,000: m=2,759
#    |H|=10,000: m=3,219
#    |H|=100,000: m=3,679
#    |H|=1,000,000: m=4,139

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

import numpy as np
from typing import List, Tuple, Set, Optional
from itertools import combinations
 
class ERMLearner:
    """
    Empirical Risk Minimization for monotone conjunctions.
    
    For monotone conjunctions, ERM is tractable because we can 
    efficiently find the consistent hypothesis (if one exists).
    """
    
    def __init__(self, n: int):
        """
        Initialize ERM learner for n Boolean variables.
        
        Hypothesis class: All subsets of {0, 1, ..., n-1}
        Each subset S represents the conjunction ∧_{i ∈ S} xᵢ
        """
        self.n = n
        self.hypothesis: Set[int] = set()
    
    def learn(self, samples: List[Tuple[np.ndarray, int]]) -> Set[int]:
        """
        Find the ERM hypothesis: consistent hypothesis with 
        minimum training error.
        
        For conjunctions, the 'tightest' consistent hypothesis is:
        - Start with all variables
        - Remove any variable that's 0 in a positive example
        """
        # Start with all variables
        self.hypothesis = set(range(self.n))
        
        # Remove variables inconsistent with positive examples
        for x, label in samples:
            if label == 1:  # Positive example
                for i in range(self.n):
                    if x[i] == 0 and i in self.hypothesis:
                        self.hypothesis.remove(i)
        
        # Verify: this is the UNIQUE consistent hypothesis (if any exists)
        # It has err_S(h) = 0 if data is consistent with some conjunction
        
        training_error = self.compute_training_error(samples)
        return self.hypothesis
    
    def compute_training_error(self, samples: List[Tuple[np.ndarray, int]]) -> float:
        """Compute training error of current hypothesis."""
        errors = sum(1 for x, label in samples if self.predict(x) != label)
        return errors / len(samples) if samples else 0.0
    
    def predict(self, x: np.ndarray) -> int:
        """Predict using current hypothesis."""
        return int(all(x[i] == 1 for i in self.hypothesis))
 
 
def demonstrate_erm_learning():
    """Demonstrate ERM learning for conjunctions."""
    np.random.seed(42)
    
    n = 6  # 6 Boolean variables
    
    # True target: x0 ∧ x2 ∧ x4 (conjunction of even-indexed vars)
    target = {0, 2, 4}
    
    def true_concept(x):
        return int(all(x[i] == 1 for i in target))
    
    # Generate training data from uniform distribution
    def generate_samples(m: int):
        samples = []
        for _ in range(m):
            x = np.random.randint(0, 2, n)
            label = true_concept(x)
            samples.append((x, label))
        return samples
    
    # Learn with increasing sample sizes
    print("ERM Learning for Conjunctions")
    print("Target: x0 ∧ x2 ∧ x4")
    print("-" * 50)
    
    for m in [5, 10, 20, 50, 100]:
        samples = generate_samples(m)
        learner = ERMLearner(n)
        hypothesis = learner.learn(samples)
        
        # Count positive examples (these determine what we can eliminate)
        n_positive = sum(1 for _, label in samples if label == 1)
        
        # Training error (should be 0 for ERM on realizable data)
        train_err = learner.compute_training_error(samples)
        
        # Is hypothesis correct?
        correct = hypothesis == target
        
        print(f"m={m:3d}: "
              f"Positive examples: {n_positive:2d}, "
              f"Learned: {sorted(hypothesis)}, "
              f"Correct: {correct}")
    
    # With enough positive examples, ERM converges to target
    # Variables NOT in target get eliminated by positive examples
    # Variables IN target never get eliminated
 
demonstrate_erm_learning()
 
# Sample Output:
# m=  5: Positive examples:  0, Learned: [0, 1, 2, 3, 4, 5], Correct: False
# m= 10: Positive examples:  1, Learned: [0, 2, 4], Correct: True
# m= 20: Positive examples:  2, Learned: [0, 2, 4], Correct: True
# m= 50: Positive examples:  4, Learned: [0, 2, 4], Correct: True
# m=100: Positive examples: 11, Learned: [0, 2, 4], Correct: True

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

import numpy as np
from math import log2, ceil
 
def information_theoretic_lower_bound():
    """
    Illustrate the information-theoretic lower bound for PAC learning.
    
    Key insight: To identify the correct hypothesis from |H| options,
    we need log₂(|H|) bits of information. Each labeled example
    provides at most 1 bit.
    
    Factor of 1/ε comes from: hypotheses might differ on only ε
    fraction of the domain. To detect differences with probability
    > 1/2, we need ~1/ε samples in that ε-measure region.
    """
    
    print("Information-Theoretic Lower Bound Analysis")
    print("=" * 60)
    
    # Scenario: Hidden hypothesis among |H| candidates
    for H_size in [16, 256, 4096, 65536]:
        bits_needed = log2(H_size)
        
        # Each example: ~1 bit of information
        # To achieve ε error, need to distinguish hypotheses
        # differing on ε fraction of domain
        
        for epsilon in [0.10, 0.05, 0.01]:
            # Rough lower bound: (1/ε) * log(|H|)
            lower_bound = (1/epsilon) * np.log(H_size)
            
            print(f"|H|={H_size:>6}, ε={epsilon:.2f}: "
                  f"≥{lower_bound:,.0f} samples needed")
        print()
    
    # Adversarial construction intuition
    print("
Adversarial Argument:")
    print("-" * 60)
    print("""
    1. Adversary picks c* uniformly from H after seeing algorithm
    2. Algorithm must output h ≈ c* with high probability
    3. But algorithm's choice depends only on m labeled examples
    4. Each example rules out at most half the hypotheses
    5. After m samples, at least |H|/2^m hypotheses remain consistent
    6. For |H|/2^m < 2, need m ≥ log₂(|H|) - 1
    
    The 1/ε factor: if hypotheses differ by only ε on distribution D,
    probability of a random sample distinguishing them is only ε.
    Need 1/ε samples to see the distinguishing region.
    """)
 
information_theoretic_lower_bound()