Limitations

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

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt
 
def analyze_holdout_variance(
    X, y,
    n_splits: int = 100,
    test_size: float = 0.2,
    model_class=RandomForestClassifier
):
    """
    Empirically measure the variance of holdout validation.
    Runs many random splits and measures result variability.
    """
    accuracies = []
    
    for seed in range(n_splits):
        # Different random split each time
        X_train, X_test, y_train, y_test = train_test_split(
            X, y, test_size=test_size, random_state=seed, stratify=y
        )
        
        # Train and evaluate
        model = model_class(random_state=42, n_estimators=50)
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        accuracies.append(accuracy_score(y_test, y_pred))
    
    accuracies = np.array(accuracies)
    
    # Statistics
    results = {
        'mean': np.mean(accuracies),
        'std': np.std(accuracies),
        'min': np.min(accuracies),
        'max': np.max(accuracies),
        'range': np.max(accuracies) - np.min(accuracies),
        '95_ci': (np.percentile(accuracies, 2.5), np.percentile(accuracies, 97.5)),
        'all_accuracies': accuracies
    }
    
    print(f"Holdout Variance Analysis ({n_splits} random splits)")
    print(f"=" * 50)
    print(f"Mean accuracy:     {results['mean']:.4f}")
    print(f"Std deviation:     {results['std']:.4f}")
    print(f"Range:             {results['min']:.4f} - {results['max']:.4f}")
    print(f"Total spread:      {results['range']:.4f}")
    print(f"95% CI:            [{results['95_ci'][0]:.4f}, {results['95_ci'][1]:.4f}]")
    
    return results
 
 
# Example on synthetic data
from sklearn.datasets import make_classification
 
X, y = make_classification(
    n_samples=500,
    n_features=20,
    n_informative=10,
    n_classes=2,
    random_state=42
)
 
results = analyze_holdout_variance(X, y, n_splits=100)
 
# The spread shows how unreliable a single holdout is
print(f"
A single holdout could give any value in [{results['min']:.1%}, {results['max']:.1%}]")
print(f"That's a {results['range']:.1%} spread - too wide for reliable comparison!")

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

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from sklearn.datasets import make_classification
 
def compare_models_sensitivity(X, y, n_splits=50):
    """
    Demonstrate how model rankings change across random splits.
    """
    models = {
        'Logistic Regression': LogisticRegression(max_iter=1000, random_state=42),
        'Random Forest': RandomForestClassifier(n_estimators=50, random_state=42),
        'Gradient Boosting': GradientBoostingClassifier(n_estimators=50, random_state=42)
    }
    
    results = {name: [] for name in models}
    rankings = []  # Who won each split
    
    for seed in range(n_splits):
        X_train, X_test, y_train, y_test = train_test_split(
            X, y, test_size=0.2, random_state=seed, stratify=y
        )
        
        split_scores = {}
        for name, model in models.items():
            model_clone = model.__class__(**model.get_params())
            model_clone.fit(X_train, y_train)
            score = accuracy_score(y_test, model_clone.predict(X_test))
            results[name].append(score)
            split_scores[name] = score
        
        # Rank models for this split
        ranking = sorted(split_scores.items(), key=lambda x: -x[1])
        rankings.append([name for name, _ in ranking])
    
    # Analyze ranking stability
    print("Model Comparison Across Random Splits")
    print("=" * 60)
    
    for name in models:
        scores = results[name]
        print(f"
{name}:")
        print(f"  Mean: {np.mean(scores):.4f} ± {np.std(scores):.4f}")
        print(f"  Range: [{np.min(scores):.4f}, {np.max(scores):.4f}]")
    
    # Count how often each model ranked first
    first_place_counts = {}
    for ranking in rankings:
        winner = ranking[0]
        first_place_counts[winner] = first_place_counts.get(winner, 0) + 1
    
    print(f"
First Place Frequency (out of {n_splits} splits):")
    for name, count in sorted(first_place_counts.items(), key=lambda x: -x[1]):
        print(f"  {name}: {count} ({count/n_splits:.1%})")
    
    # How often do rankings agree with the average-based ranking?
    avg_ranking = sorted(
        [(name, np.mean(scores)) for name, scores in results.items()],
        key=lambda x: -x[1]
    )
    expected_order = [name for name, _ in avg_ranking]
    
    matches = sum(1 for r in rankings if r == expected_order)
    print(f"
Rankings match average-based order: {matches}/{n_splits} ({matches/n_splits:.1%})")
    
    return results, rankings
 
 
# Generate data
X, y = make_classification(
    n_samples=500, n_features=20, n_informative=10,
    n_redundant=5, n_classes=2, random_state=42
)
 
results, rankings = compare_models_sensitivity(X, y, n_splits=50)

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

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score
from scipy import stats
 
def adversarial_validation(X_train, X_test, n_estimators=100):
    """
    Check if train and test distributions are distinguishable.
    
    Returns AUC close to 0.5 if distributions are similar.
    Returns AUC >> 0.5 if distributions are different.
    """
    # Create labels: 0 = train, 1 = test
    y_combined = np.array([0] * len(X_train) + [1] * len(X_test))
    X_combined = np.vstack([X_train, X_test])
    
    # Split for validation
    X_adv_train, X_adv_val, y_adv_train, y_adv_val = train_test_split(
        X_combined, y_combined, test_size=0.2, random_state=42, stratify=y_combined
    )
    
    # Train classifier to distinguish train from test
    clf = RandomForestClassifier(n_estimators=n_estimators, random_state=42)
    clf.fit(X_adv_train, y_adv_train)
    
    # Evaluate
    y_proba = clf.predict_proba(X_adv_val)[:, 1]
    auc = roc_auc_score(y_adv_val, y_proba)
    
    # Feature importances (which features differ most)
    importances = clf.feature_importances_
    
    # Interpretation
    if auc < 0.55:
        interpretation = "GOOD: Distributions appear similar"
    elif auc < 0.65:
        interpretation = "WARNING: Some distribution difference detected"
    else:
        interpretation = "ALERT: Significant distribution mismatch!"
    
    return {
        'auc': auc,
        'interpretation': interpretation,
        'feature_importances': importances
    }
 
 
def feature_distribution_comparison(X_train, X_test, feature_names=None):
    """
    Compare feature distributions between train and test.
    Returns p-values for each feature (low p-value = different distributions).
    """
    n_features = X_train.shape[1]
    if feature_names is None:
        feature_names = [f'Feature_{i}' for i in range(n_features)]
    
    results = []
    
    for i in range(n_features):
        train_values = X_train[:, i]
        test_values = X_test[:, i]
        
        # Kolmogorov-Smirnov test
        ks_stat, ks_pval = stats.ks_2samp(train_values, test_values)
        
        # Basic statistics
        results.append({
            'feature': feature_names[i],
            'train_mean': np.mean(train_values),
            'test_mean': np.mean(test_values),
            'mean_diff': np.mean(test_values) - np.mean(train_values),
            'ks_statistic': ks_stat,
            'ks_pvalue': ks_pval,
            'significant': ks_pval < 0.01
        })
    
    # Sort by significance
    results.sort(key=lambda x: x['ks_pvalue'])
    
    return results
 
 
# Usage
result = adversarial_validation(X_train, X_test)
print(f"Adversarial AUC: {result['auc']:.4f}")
print(f"Interpretation: {result['interpretation']}")

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
115
116

import numpy as np
from statsmodels.stats.proportion import proportion_confint
from scipy import stats
 
class HoldoutWithUncertainty:
    """
    Holdout evaluation with proper uncertainty quantification.
    """
    
    def __init__(self, y_true, y_pred, y_proba=None):
        self.y_true = np.array(y_true)
        self.y_pred = np.array(y_pred)
        self.y_proba = np.array(y_proba) if y_proba is not None else None
        self.n_samples = len(y_true)
    
    def accuracy_with_ci(self, confidence=0.95, method='wilson'):
        """
        Compute accuracy with confidence interval.
        
        Methods:
        - 'wilson': Wilson score interval (recommended)
        - 'agresti_coull': Agresti-Coull interval
        - 'normal': Normal approximation (use only for large n)
        """
        correct = np.sum(self.y_pred == self.y_true)
        accuracy = correct / self.n_samples
        
        alpha = 1 - confidence
        ci_low, ci_high = proportion_confint(correct, self.n_samples, 
                                              alpha=alpha, method=method)
        
        return {
            'accuracy': accuracy,
            'ci_low': ci_low,
            'ci_high': ci_high,
            'ci_width': ci_high - ci_low,
            'n_samples': self.n_samples,
            'method': method
        }
    
    def bootstrap_metrics(self, n_bootstrap=1000, metrics=['accuracy']):
        """
        Bootstrap estimation of multiple metrics.
        """
        results = {metric: [] for metric in metrics}
        
        for _ in range(n_bootstrap):
            indices = np.random.choice(self.n_samples, self.n_samples, replace=True)
            y_true_boot = self.y_true[indices]
            y_pred_boot = self.y_pred[indices]
            
            if 'accuracy' in metrics:
                results['accuracy'].append(np.mean(y_pred_boot == y_true_boot))
        
        summary = {}
        for metric in metrics:
            values = np.array(results[metric])
            summary[metric] = {
                'mean': np.mean(values),
                'std': np.std(values),
                'ci_low': np.percentile(values, 2.5),
                'ci_high': np.percentile(values, 97.5)
            }
        
        return summary
    
    def minimum_detectable_difference(self, power=0.8, alpha=0.05):
        """
        What performance difference can we reliably detect with this test size?
        """
        # Two-proportion z-test power calculation
        # Approximate: MDD ≈ 2.8 * sqrt(p*(1-p)/n) for 80% power, alpha=0.05
        p = np.mean(self.y_pred == self.y_true)
        mdd = 2.8 * np.sqrt(p * (1 - p) / self.n_samples)
        
        return {
            'minimum_detectable_difference': mdd,
            'interpretation': f"Need >{mdd:.1%} accuracy difference to detect with {power:.0%} power"
        }
    
    def generate_report(self) -> str:
        """Generate a complete uncertainty report."""
        acc = self.accuracy_with_ci()
        boot = self.bootstrap_metrics()
        mdd = self.minimum_detectable_difference()
        
        lines = [
            "=" * 60,
            "HOLDOUT EVALUATION REPORT WITH UNCERTAINTY",
            "=" * 60,
            f"",
            f"Sample size: {self.n_samples}",
            f"",
            f"ACCURACY:",
            f"  Point estimate: {acc['accuracy']:.4f}",
            f"  95% CI (Wilson): [{acc['ci_low']:.4f}, {acc['ci_high']:.4f}]",
            f"  CI width: {acc['ci_width']:.4f}",
            f"",
            f"BOOTSTRAP (1000 resamples):",
            f"  Mean: {boot['accuracy']['mean']:.4f}",
            f"  Std:  {boot['accuracy']['std']:.4f}",
            f"  95% CI: [{boot['accuracy']['ci_low']:.4f}, {boot['accuracy']['ci_high']:.4f}]",
            f"",
            f"STATISTICAL POWER:",
            f"  {mdd['interpretation']}",
            f"",
            "=" * 60
        ]
        
        return "
".join(lines)
 
 
# Usage
evaluator = HoldoutWithUncertainty(y_test, y_pred)
print(evaluator.generate_report())