Bias-Variance Tradeoff in Cross-Validation

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

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import learning_curve, cross_val_score
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
 
# Generate data
np.random.seed(42)
X, y = make_classification(n_samples=1000, n_features=20, 
                           n_informative=10, random_state=42)
 
# Compute learning curve
train_sizes, train_scores, val_scores = learning_curve(
    RandomForestClassifier(n_estimators=50, random_state=42),
    X, y,
    train_sizes=np.linspace(0.1, 1.0, 10),
    cv=5,
    scoring='accuracy',
    random_state=42
)
 
# Plot learning curve
plt.figure(figsize=(10, 6))
plt.plot(train_sizes, train_scores.mean(axis=1), 'o-', label='Training score')
plt.plot(train_sizes, val_scores.mean(axis=1), 'o-', label='Validation score')
plt.fill_between(train_sizes, 
                 val_scores.mean(axis=1) - val_scores.std(axis=1),
                 val_scores.mean(axis=1) + val_scores.std(axis=1), 
                 alpha=0.2)
 
# Mark what different k values estimate
for k in [2, 5, 10]:
    train_fraction = (k-1)/k
    train_n = int(1000 * train_fraction)
    idx = np.argmin(np.abs(train_sizes - train_n))
    plt.axvline(train_n, color='red', linestyle='--', alpha=0.5)
    plt.annotate(f'k={k}
({train_fraction:.0%})', 
                xy=(train_n, val_scores.mean(axis=1)[idx]),
                xytext=(train_n-80, val_scores.mean(axis=1)[idx]+0.02))
 
plt.xlabel('Training Set Size')
plt.ylabel('Score')
plt.title('Learning Curve: What Different k Values Estimate')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
 
# Demonstrate bias empirically
print("Empirical CV estimates for different k:")
for k in [2, 3, 5, 10, 20]:
    cv_score = cross_val_score(
        RandomForestClassifier(n_estimators=50, random_state=42),
        X, y, cv=k, scoring='accuracy'
    ).mean()
    train_fraction = (k-1)/k
    print(f"k={k:2d}: CV accuracy = {cv_score:.4f} "
          f"(trains on {train_fraction:.1%} of data)")

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
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.datasets import make_classification
 
def estimate_true_error(X, y, model_class, n_test=10000, n_trials=50):
    """
    Estimate true generalization error by training on all data
    and testing on fresh samples.
    """
    errors = []
    n_samples = len(X)
    
    for _ in range(n_trials):
        # Generate fresh test data from same distribution
        X_test, y_test = make_classification(
            n_samples=n_test,
            n_features=X.shape[1],
            n_informative=10,
            random_state=np.random.randint(0, 10000)
        )
        
        # Train on all available data
        model = model_class()
        model.fit(X, y)
        
        # Evaluate on fresh data
        error = 1 - model.score(X_test, y_test)
        errors.append(error)
    
    return np.mean(errors), np.std(errors)
 
def compute_cv_bias(X, y, model_class, k_values=[2, 5, 10, 20]):
    """
    Compute CV estimates and compare to true error.
    """
    # Estimate true error (training on all n samples)
    true_error, true_std = estimate_true_error(X, y, model_class)
    print(f"True generalization error (trained on n={len(X)}): "
          f"{true_error:.4f} ± {true_std:.4f}")
    print()
    
    results = []
    for k in k_values:
        cv_scores = cross_val_score(
            model_class(), X, y, cv=k, scoring='accuracy'
        )
        cv_error = 1 - cv_scores.mean()
        cv_std = cv_scores.std()
        
        bias = cv_error - true_error
        train_fraction = (k-1)/k
        
        print(f"k={k:2d}: CV error = {cv_error:.4f} ± {cv_std:.4f}, "
              f"Bias = {bias:+.4f} ({train_fraction:.0%} train)")
        
        results.append({
            'k': k,
            'cv_error': cv_error,
            'cv_std': cv_std,
            'bias': bias,
            'train_fraction': train_fraction
        })
    
    return results
 
# Run simulation
np.random.seed(42)
X, y = make_classification(n_samples=500, n_features=20,
                           n_informative=10, random_state=42)
 
print("=" * 60)
print("Bias Analysis: Logistic Regression")
print("=" * 60)
results = compute_cv_bias(X, y, LogisticRegression)

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

import numpy as np
from sklearn.model_selection import cross_val_score, RepeatedKFold
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification
 
def analyze_cv_variance(X, y, model_class, k_values=[2, 5, 10, 20], 
                        n_repeats=100, random_seed=42):
    """
    Analyze variance of CV estimates across different k values.
    
    We run CV multiple times with different random partitions to
    estimate the full variance of the CV procedure.
    """
    np.random.seed(random_seed)
    
    results = {}
    
    for k in k_values:
        cv_means = []
        fold_scores_all = []
        
        for i in range(n_repeats):
            # Different random partition each time
            scores = cross_val_score(
                model_class(), X, y, cv=k, 
                scoring='accuracy'
            )
            cv_means.append(scores.mean())
            fold_scores_all.append(scores)
        
        cv_means = np.array(cv_means)
        fold_scores_all = np.array(fold_scores_all)  # (n_repeats, k)
        
        # Mean fold variance (within a single CV run)
        within_cv_var = np.mean([s.var() for s in fold_scores_all])
        
        # Variance of the CV estimate across repeats
        cv_estimate_var = cv_means.var()
        
        # Estimate correlation between fold scores
        # Average pairwise correlation within each CV run
        correlations = []
        for scores in fold_scores_all[:20]:  # Sample for efficiency
            if len(scores) >= 2:
                for i in range(len(scores)):
                    for j in range(i+1, len(scores)):
                        # This is not quite right - we need paired observations
                        pass
        
        results[k] = {
            'cv_estimate_mean': cv_means.mean(),
            'cv_estimate_std': cv_means.std(),
            'cv_estimate_var': cv_estimate_var,
            'within_cv_std': np.sqrt(within_cv_var),
            'coefficient_of_variation': cv_means.std() / cv_means.mean()
        }
        
        print(f"k={k:2d}: CV estimate = {cv_means.mean():.4f} ± {cv_means.std():.4f}")
        print(f"       Within-run fold std = {np.sqrt(within_cv_var):.4f}")
        print(f"       Coefficient of variation = {results[k]['coefficient_of_variation']:.4f}")
        print()
    
    return results
 
# Compare models with different stability
np.random.seed(42)
X, y = make_classification(n_samples=500, n_features=20,
                           n_informative=10, random_state=42)
 
print("=" * 60)
print("Variance Analysis: Logistic Regression (stable model)")
print("=" * 60)
lr_results = analyze_cv_variance(X, y, LogisticRegression)
 
print("=" * 60)
print("Variance Analysis: Random Forest (less stable)")  
print("=" * 60)
rf_results = analyze_cv_variance(
    X, y, 
    lambda: RandomForestClassifier(n_estimators=50, random_state=None)
)
 
# Compare fold overlap
print("=" * 60)
print("Training Set Overlap Analysis")
print("=" * 60)
n = len(X)
for k in [2, 5, 10, 20]:
    train_size = n * (k-1) / k
    overlap = (k-2) / (k-1) if k > 1 else 0  # Fraction of overlap between two folds
    print(f"k={k:2d}: Training size = {train_size:.0f} ({(k-1)/k:.0%}), "
          f"Pairwise training overlap = {overlap:.1%}")

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

import numpy as np
import time
from sklearn.model_selection import (
    LeaveOneOut, KFold, cross_val_score, cross_val_predict
)
from sklearn.linear_model import Ridge
from sklearn.tree import DecisionTreeRegressor
from sklearn.datasets import make_regression
 
# Generate regression data
np.random.seed(42)
X, y = make_regression(n_samples=200, n_features=10, noise=10, random_state=42)
 
def compare_cv_strategies(X, y, model_class, model_name):
    """Compare LOOCV and k-fold CV."""
    print(f"
{model_name}")
    print("-" * 40)
    
    # LOOCV
    loo = LeaveOneOut()
    start = time.time()
    loo_scores = cross_val_score(model_class(), X, y, cv=loo, scoring='r2')
    loo_time = time.time() - start
    
    print(f"LOOCV (k={len(X)}): R² = {loo_scores.mean():.4f} ± {loo_scores.std():.4f}")
    print(f"  Time: {loo_time:.2f}s, Iterations: {len(loo_scores)}")
    
    # Various k-fold
    for k in [5, 10, 20, 50]:
        kfold = KFold(n_splits=k, shuffle=True, random_state=42)
        start = time.time()
        kf_scores = cross_val_score(model_class(), X, y, cv=kfold, scoring='r2')
        kf_time = time.time() - start
        
        print(f"{k:2d}-Fold: R² = {kf_scores.mean():.4f} ± {kf_scores.std():.4f}")
        print(f"  Time: {kf_time:.2f}s, Iterations: {k}")
    
    return loo_scores, loo_time
 
# Compare stable vs unstable models
print("=" * 60)
print("LOOCV vs K-Fold Comparison")
print("=" * 60)
 
# Stable model (Ridge regression)
compare_cv_strategies(X, y, lambda: Ridge(alpha=1.0), "Ridge Regression (Stable)")
 
# Unstable model (Deep decision tree)
compare_cv_strategies(X, y, 
                     lambda: DecisionTreeRegressor(max_depth=None, random_state=None),
                     "Decision Tree (Unstable)")
 
# Demonstrate variance across different random partitions
print("
" + "=" * 60)
print("Variance Across Random Partitions (10-fold)")
print("=" * 60)
 
cv_means = []
for seed in range(50):
    kfold = KFold(n_splits=10, shuffle=True, random_state=seed)
    scores = cross_val_score(
        DecisionTreeRegressor(max_depth=3), X, y, cv=kfold, scoring='r2'
    )
    cv_means.append(scores.mean())
 
print(f"Mean of CV estimates: {np.mean(cv_means):.4f}")
print(f"Std of CV estimates: {np.std(cv_means):.4f}")
print(f"Range: [{np.min(cv_means):.4f}, {np.max(cv_means):.4f}]")
print("
Note: LOOCV would give the same estimate every time (deterministic)")

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

def recommend_cv_strategy(n_samples, n_features, model_type, goal,
                            compute_budget='medium'):
    """
    Recommend cross-validation strategy based on dataset and goals.
    
    Parameters:
    -----------
    n_samples : int
        Number of training samples
    n_features : int
        Number of features
    model_type : str
        'linear', 'tree', 'ensemble', 'neural_net', 'stable', 'unstable'
    goal : str
        'model_selection', 'performance_estimation', 'feature_selection'
    compute_budget : str
        'low', 'medium', 'high'
        
    Returns:
    --------
    dict with recommended strategy
    """
    recommendation = {
        'k': None,
        'n_repeats': 1,
        'reasoning': []
    }
    
    # Base k selection by sample size
    if n_samples < 50:
        recommendation['k'] = n_samples  # LOOCV
        recommendation['reasoning'].append(
            f"Very small dataset (n={n_samples}): Use LOOCV to maximize training data"
        )
    elif n_samples < 200:
        recommendation['k'] = 10
        recommendation['reasoning'].append(
            f"Small dataset (n={n_samples}): Use 10-fold to balance bias/variance"
        )
    elif n_samples < 5000:
        recommendation['k'] = 5 if compute_budget == 'low' else 10
        recommendation['reasoning'].append(
            f"Medium dataset (n={n_samples}): k=5 or 10 both reasonable"
        )
    else:
        recommendation['k'] = 5
        recommendation['reasoning'].append(
            f"Large dataset (n={n_samples}): k=5 sufficient, bias minimal"
        )
    
    # Adjust for model stability
    unstable_models = ['tree', 'neural_net', 'unstable']
    if model_type in unstable_models:
        recommendation['n_repeats'] = max(3, recommendation['n_repeats'])
        recommendation['reasoning'].append(
            f"Unstable model ({model_type}): Use repeated CV to reduce variance"
        )
    
    # Adjust for goal
    if goal == 'performance_estimation':
        if recommendation['k'] < 10:
            recommendation['k'] = min(10, n_samples // 10)
            recommendation['reasoning'].append(
                "Performance estimation: Prefer larger k to reduce pessimistic bias"
            )
    elif goal == 'feature_selection':
        recommendation['n_repeats'] = max(5, recommendation['n_repeats'])
        recommendation['reasoning'].append(
            "Feature selection: Multiple repeats for stable importance rankings"
        )
    
    # Compute budget constraints
    if compute_budget == 'low':
        recommendation['k'] = min(recommendation['k'], 5)
        recommendation['n_repeats'] = 1
        recommendation['reasoning'].append(
            "Low compute budget: Limiting to k=5, single repeat"
        )
    elif compute_budget == 'high':
        recommendation['n_repeats'] = max(5, recommendation['n_repeats'])
        recommendation['reasoning'].append(
            "High compute budget: Using multiple repeats for robustness"
        )
    
    return recommendation
 
# Example usage
scenarios = [
    {'n_samples': 50, 'n_features': 10, 'model_type': 'linear', 
     'goal': 'model_selection'},
    {'n_samples': 500, 'n_features': 50, 'model_type': 'ensemble', 
     'goal': 'performance_estimation'},
    {'n_samples': 200, 'n_features': 100, 'model_type': 'tree', 
     'goal': 'feature_selection'},
    {'n_samples': 10000, 'n_features': 20, 'model_type': 'neural_net', 
     'goal': 'model_selection', 'compute_budget': 'low'},
]
 
for scenario in scenarios:
    print("Scenario:", scenario)
    rec = recommend_cv_strategy(**scenario)
    print(f"  Recommended: {rec['k']}-fold with {rec['n_repeats']} repeats")
    for reason in rec['reasoning']:
        print(f"    - {reason}")
    print()