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

import numpy as np
from sklearn.ensemble import RandomForestClassifier, BaggingClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.datasets import make_classification
from itertools import combinations
 
def measure_tree_correlation(estimators, X, proba_class=1):
    """
    Measure average pairwise correlation between tree predictions.
    
    Parameters:
    -----------
    estimators : list of fitted trees
    X : test data to make predictions on
    proba_class : which class probability to use (for classification)
    
    Returns:
    --------
    float : average pairwise Pearson correlation
    """
    # Get predictions from each tree
    predictions = []
    for tree in estimators:
        if hasattr(tree, 'predict_proba'):
            pred = tree.predict_proba(X)[:, proba_class]
        else:
            pred = tree.predict(X).astype(float)
        predictions.append(pred)
    
    predictions = np.array(predictions)  # (n_trees, n_samples)
    
    # Calculate all pairwise correlations
    correlations = []
    for i, j in combinations(range(len(estimators)), 2):
        # Pearson correlation between tree i and tree j predictions
        corr = np.corrcoef(predictions[i], predictions[j])[0, 1]
        if not np.isnan(corr):
            correlations.append(corr)
    
    return np.mean(correlations), np.std(correlations)
 
 
def compare_correlation_by_method(X, y, n_estimators=50):
    """
    Compare tree correlation across different ensemble methods.
    """
    # Method 1: Bagged trees (no feature randomization)
    bagging = BaggingClassifier(
        estimator=DecisionTreeClassifier(),
        n_estimators=n_estimators,
        bootstrap=True,
        max_features=1.0,  # Using all features
        random_state=42
    )
    bagging.fit(X, y)
    
    # Method 2: Random Forest (sqrt features)
    rf_sqrt = RandomForestClassifier(
        n_estimators=n_estimators,
        max_features='sqrt',
        random_state=42
    )
    rf_sqrt.fit(X, y)
    
    # Method 3: Random Forest (log2 features)
    rf_log2 = RandomForestClassifier(
        n_estimators=n_estimators,
        max_features='log2',
        random_state=42
    )
    rf_log2.fit(X, y)
    
    # Method 4: Extremely aggressive randomization
    rf_extreme = RandomForestClassifier(
        n_estimators=n_estimators,
        max_features=2,  # Just 2 features per split
        random_state=42
    )
    rf_extreme.fit(X, y)
    
    results = {}
    
    print("Tree Correlation Comparison")
    print("=" * 60)
    
    for name, model in [
        ("Bagged Trees (all features)", bagging),
        ("Random Forest (sqrt)", rf_sqrt),
        ("Random Forest (log2)", rf_log2),
        ("Extreme Random (m=2)", rf_extreme),
    ]:
        mean_corr, std_corr = measure_tree_correlation(model.estimators_, X)
        accuracy = model.score(X, y)
        results[name] = {'correlation': mean_corr, 'accuracy': accuracy}
        print(f"{name:30} | ρ = {mean_corr:.4f} (±{std_corr:.4f}) | Acc = {accuracy:.4f}")
    
    return results
 
 
# Example with synthetic data
X, y = make_classification(
    n_samples=1000,
    n_features=50,
    n_informative=10,
    n_redundant=10,
    n_clusters_per_class=3,
    random_state=42
)
 
compare_correlation_by_method(X, y)

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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_iris
from scipy.stats import pearsonr
from collections import defaultdict
 
def comprehensive_correlation_analysis(rf, X, y):
    """
    Perform comprehensive correlation analysis on a fitted Random Forest.
    """
    n_samples = X.shape[0]
    n_trees = len(rf.estimators_)
    
    # ============================================
    # Method 1: Prediction Correlation
    # ============================================
    # Get probability predictions from each tree
    tree_probas = np.array([
        tree.predict_proba(X)[:, 1] 
        for tree in rf.estimators_
    ])  # Shape: (n_trees, n_samples)
    
    # Compute pairwise correlations
    correlations = []
    for i in range(n_trees):
        for j in range(i+1, n_trees):
            r = np.corrcoef(tree_probas[i], tree_probas[j])[0, 1]
            if not np.isnan(r):
                correlations.append(r)
    
    print("=" * 60)
    print("Random Forest Correlation Diagnostics")
    print("=" * 60)
    print(f"
1. PREDICTION CORRELATION")
    print(f"   Average pairwise ρ: {np.mean(correlations):.4f}")
    print(f"   Std pairwise ρ: {np.std(correlations):.4f}")
    print(f"   Min/Max ρ: [{np.min(correlations):.4f}, {np.max(correlations):.4f}]")
    
    # ============================================
    # Method 2: Vote Disagreement (Classification)
    # ============================================
    tree_predictions = np.array([
        tree.predict(X) for tree in rf.estimators_
    ])  # Shape: (n_trees, n_samples)
    
    # For each sample, what fraction of trees agree with majority?
    majority_agreement = []
    for i in range(n_samples):
        votes = tree_predictions[:, i]
        majority_class = np.bincount(votes.astype(int)).argmax()
        agreement = (votes == majority_class).mean()
        majority_agreement.append(agreement)
    
    avg_agreement = np.mean(majority_agreement)
    diversity_score = 1 - avg_agreement  # Higher = more diverse
    
    print(f"
2. VOTE DIVERSITY")
    print(f"   Average majority agreement: {avg_agreement:.4f}")
    print(f"   Diversity score: {diversity_score:.4f}")
    print(f"   (Higher diversity = lower correlation)")
    
    # ============================================
    # Method 3: Structural Analysis
    # ============================================
    # Count feature usage at each depth
    depth_feature_usage = defaultdict(lambda: defaultdict(int))
    
    for tree in rf.estimators_:
        tree_struct = tree.tree_
        n_nodes = tree_struct.node_count
        
        # DFS to track depth
        stack = [(0, 0)]  # (node_id, depth)
        while stack:
            node_id, depth = stack.pop()
            
            if tree_struct.feature[node_id] >= 0:  # Not a leaf
                feature = tree_struct.feature[node_id]
                depth_feature_usage[depth][feature] += 1
                
                stack.append((tree_struct.children_left[node_id], depth + 1))
                stack.append((tree_struct.children_right[node_id], depth + 1))
    
    print(f"
3. STRUCTURAL DIVERSITY (Feature usage at depth)")
    for depth in sorted(depth_feature_usage.keys())[:5]:  # First 5 depths
        features_used = depth_feature_usage[depth]
        unique_features = len(features_used)
        total_splits = sum(features_used.values())
        top_feature = max(features_used.items(), key=lambda x: x[1])
        top_feature_pct = top_feature[1] / total_splits * 100
        
        print(f"   Depth {depth}: {unique_features} unique features, "
              f"top feature {top_feature[0]} used {top_feature_pct:.1f}% of splits")
    
    # ============================================
    # Method 4: Implied Variance Floor
    # ============================================
    ensemble_var = np.var(rf.predict_proba(X)[:, 1])
    individual_vars = np.var(tree_probas, axis=1)
    avg_individual_var = np.mean(individual_vars)
    
    # Implied correlation from variance floor
    # Var_ensemble = ρ * σ² as T → ∞
    implied_rho = ensemble_var / avg_individual_var if avg_individual_var > 0 else 0
    
    print(f"
4. VARIANCE ANALYSIS")
    print(f"   Avg individual tree variance: {avg_individual_var:.6f}")
    print(f"   Ensemble variance: {ensemble_var:.6f}")
    print(f"   Variance reduction ratio: {ensemble_var/avg_individual_var:.4f}")
    print(f"   Implied minimum correlation: {implied_rho:.4f}")
    
    return {
        'avg_correlation': np.mean(correlations),
        'diversity_score': diversity_score,
        'variance_reduction': ensemble_var/avg_individual_var
    }
 
 
# Example usage
from sklearn.datasets import make_classification
 
X, y = make_classification(
    n_samples=500, n_features=30,
    n_informative=10, n_redundant=5,
    random_state=42
)
 
rf = RandomForestClassifier(n_estimators=100, max_features='sqrt', random_state=42)
rf.fit(X, y)
 
results = comprehensive_correlation_analysis(rf, X, y)

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

from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_score
import numpy as np
 
def compare_rf_et(X, y, n_estimators=100):
    """
    Compare Random Forest vs Extra-Trees on correlation and accuracy.
    """
    rf = RandomForestClassifier(n_estimators=n_estimators, random_state=42)
    et = ExtraTreesClassifier(n_estimators=n_estimators, random_state=42)
    
    rf.fit(X, y)
    et.fit(X, y)
    
    def measure_avg_correlation(estimators, X):
        probas = np.array([t.predict_proba(X)[:, 1] for t in estimators])
        corrs = []
        for i in range(len(estimators)):
            for j in range(i+1, len(estimators)):
                r = np.corrcoef(probas[i], probas[j])[0, 1]
                if not np.isnan(r):
                    corrs.append(r)
        return np.mean(corrs)
    
    rf_corr = measure_avg_correlation(rf.estimators_, X)
    et_corr = measure_avg_correlation(et.estimators_, X)
    
    rf_cv = cross_val_score(rf, X, y, cv=5).mean()
    et_cv = cross_val_score(et, X, y, cv=5).mean()
    
    print("Random Forest vs Extra-Trees")
    print("=" * 50)
    print(f"Random Forest:  ρ = {rf_corr:.4f}, CV Accuracy = {rf_cv:.4f}")
    print(f"Extra-Trees:    ρ = {et_corr:.4f}, CV Accuracy = {et_cv:.4f}")
    print(f"Correlation reduction: {(rf_corr - et_corr)/rf_corr*100:.1f}%")
 
 
X, y = make_classification(
    n_samples=1000, n_features=50,
    n_informative=15, random_state=42
)
 
compare_rf_et(X, y)

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

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
import numpy as np
 
def diagnose_and_tune_correlation(X, y):
    """
    Diagnose correlation issues and suggest remediation.
    """
    n_features = X.shape[1]
    
    print("Random Forest Correlation Diagnostic")
    print("=" * 60)
    
    # Test different configurations
    configs = [
        ('Default (sqrt)', {'max_features': 'sqrt'}),
        ('p/3', {'max_features': max(1, n_features // 3)}),
        ('p/4', {'max_features': max(1, n_features // 4)}),
        ('log2', {'max_features': 'log2'}),
        ('3 features', {'max_features': 3}),
    ]
    
    results = []
    for name, params in configs:
        rf = RandomForestClassifier(n_estimators=100, random_state=42, **params)
        rf.fit(X, y)
        
        # Measure correlation
        probas = np.array([t.predict_proba(X)[:, 1] for t in rf.estimators_])
        corrs = []
        for i in range(100):
            for j in range(i+1, 100):
                r = np.corrcoef(probas[i], probas[j])[0, 1]
                if not np.isnan(r):
                    corrs.append(r)
        avg_corr = np.mean(corrs)
        
        # Cross-validation score
        cv_score = cross_val_score(rf, X, y, cv=5).mean()
        
        results.append((name, avg_corr, cv_score))
        print(f"{name:20} | ρ = {avg_corr:.4f} | CV = {cv_score:.4f}")
    
    # Find best configuration
    best_idx = np.argmax([r[2] for r in results])
    print(f"
Recommendation: {results[best_idx][0]}")
    print(f"  Achieved lowest correlation with highest CV score")
    
    # Check if correlation is problematic
    default_corr = results[0][1]
    if default_corr > 0.6:
        print("
⚠️  WARNING: High tree correlation detected!")
        print("   Consider reducing max_features or using Extra-Trees")
    elif default_corr > 0.4:
        print("
📊 MODERATE correlation - some room for improvement")
    else:
        print("
✅ LOW correlation - ensemble is well-diversified")
    
    return results
 
 
# Example
from sklearn.datasets import make_classification
X, y = make_classification(
    n_samples=1000, n_features=50,
    n_informative=10, n_redundant=10,
    random_state=42
)
 
diagnose_and_tune_correlation(X, y)