Bootstrap Resampling

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

import numpy as np
from typing import Callable, List, Tuple
 
def bootstrap_resample(data: np.ndarray, 
                       statistic_fn: Callable[[np.ndarray], float],
                       n_bootstrap: int = 10000,
                       random_state: int = None) -> np.ndarray:
    """
    Perform non-parametric bootstrap resampling.
    
    Parameters
    ----------
    data : np.ndarray
        Original sample data of shape (n_samples,) or (n_samples, n_features)
    statistic_fn : Callable
        Function that computes the statistic of interest from a sample
    n_bootstrap : int
        Number of bootstrap iterations (B)
    random_state : int, optional
        Random seed for reproducibility
    
    Returns
    -------
    np.ndarray
        Array of B bootstrap statistics
    
    The algorithm:
    1. For each bootstrap iteration b = 1, ..., B:
       a. Draw n samples WITH REPLACEMENT from data
       b. Compute statistic on resampled data
    2. Return all B statistics
    """
    rng = np.random.RandomState(random_state)
    n_samples = len(data)
    
    bootstrap_statistics = np.zeros(n_bootstrap)
    
    for b in range(n_bootstrap):
        # Step 1: Sample WITH REPLACEMENT
        # Each observation has probability 1/n of being selected each time
        indices = rng.choice(n_samples, size=n_samples, replace=True)
        bootstrap_sample = data[indices]
        
        # Step 2: Compute statistic on bootstrap sample
        bootstrap_statistics[b] = statistic_fn(bootstrap_sample)
    
    return bootstrap_statistics
 
 
def compute_bootstrap_standard_error(data: np.ndarray,
                                     statistic_fn: Callable[[np.ndarray], float],
                                     n_bootstrap: int = 10000) -> Tuple[float, float]:
    """
    Compute bootstrap estimate of standard error.
    
    This is the most basic bootstrap application: estimating the
    standard deviation of the sampling distribution.
    
    Returns
    -------
    Tuple[float, float]
        (original_statistic, bootstrap_standard_error)
    """
    # Original statistic
    theta_hat = statistic_fn(data)
    
    # Bootstrap distribution
    theta_stars = bootstrap_resample(data, statistic_fn, n_bootstrap)
    
    # Bootstrap SE is just the standard deviation of bootstrap statistics
    se_boot = np.std(theta_stars, ddof=1)
    
    return theta_hat, se_boot
 
 
# Example: Bootstrap SE of the sample mean
if __name__ == "__main__":
    # Generate sample data
    np.random.seed(42)
    sample = np.random.exponential(scale=2.0, size=50)
    
    # Statistic: sample mean
    theta_hat, se_boot = compute_bootstrap_standard_error(
        sample, 
        statistic_fn=np.mean,
        n_bootstrap=10000
    )
    
    # Compare with analytical SE (known for the mean)
    se_analytical = np.std(sample, ddof=1) / np.sqrt(len(sample))
    
    print(f"Sample mean: {theta_hat:.4f}")
    print(f"Bootstrap SE: {se_boot:.4f}")
    print(f"Analytical SE: {se_analytical:.4f}")
    print(f"Ratio (Boot/Analytical): {se_boot/se_analytical:.4f}")

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

import numpy as np
from collections import Counter
 
def analyze_bootstrap_inclusion(n_samples: int, 
                                 n_bootstrap: int = 10000,
                                 random_state: int = 42) -> dict:
    """
    Empirically verify bootstrap inclusion probabilities.
    
    Parameters
    ----------
    n_samples : int
        Size of the dataset
    n_bootstrap : int
        Number of bootstrap samples to generate
    random_state : int
        Random seed
        
    Returns
    -------
    dict
        Statistics about inclusion patterns
    """
    rng = np.random.RandomState(random_state)
    
    # Track inclusion across bootstrap samples
    inclusion_counts = np.zeros((n_bootstrap, n_samples))
    unique_counts = np.zeros(n_bootstrap)
    
    for b in range(n_bootstrap):
        # Generate bootstrap indices
        indices = rng.choice(n_samples, size=n_samples, replace=True)
        
        # Count unique observations
        unique_counts[b] = len(np.unique(indices))
        
        # Track which observations were included
        for idx in indices:
            inclusion_counts[b, idx] += 1
    
    # Theoretical values
    theoretical_exclusion = (1 - 1/n_samples)**n_samples
    theoretical_inclusion = 1 - theoretical_exclusion
    
    # Empirical values
    empirical_unique_mean = np.mean(unique_counts)
    empirical_unique_std = np.std(unique_counts)
    
    # For each observation, count how often it was excluded
    never_included_per_bootstrap = np.sum(inclusion_counts == 0, axis=1)
    empirical_exclusion_rate = np.mean(never_included_per_bootstrap) / n_samples
    
    return {
        'n_samples': n_samples,
        'theoretical_exclusion': theoretical_exclusion,
        'theoretical_inclusion': theoretical_inclusion,
        'empirical_exclusion_rate': empirical_exclusion_rate,
        'mean_unique_obs': empirical_unique_mean,
        'std_unique_obs': empirical_unique_std,
        'theoretical_unique': n_samples * theoretical_inclusion,
    }
 
 
# Demonstrate the 1/e phenomenon
if __name__ == "__main__":
    import math
    
    print("Bootstrap Inclusion Probability Analysis")
    print("=" * 60)
    print(f"Theoretical limit (n→∞): P(excluded) = 1/e ≈ {1/math.e:.6f}")
    print()
    
    for n in [10, 50, 100, 500, 1000]:
        results = analyze_bootstrap_inclusion(n, n_bootstrap=5000)
        print(f"n = {n:4d}: Theoretical = {results['theoretical_exclusion']:.4f}, "
              f"Empirical = {results['empirical_exclusion_rate']:.4f}, "
              f"Unique obs = {results['mean_unique_obs']:.1f} ± {results['std_unique_obs']:.1f}")

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

import numpy as np
from sklearn.base import clone
from typing import List, Tuple, Callable
 
def out_of_bootstrap_error(X: np.ndarray,
                            y: np.ndarray,
                            model,
                            n_bootstrap: int = 200,
                            scoring_fn: Callable = None,
                            random_state: int = 42) -> dict:
    """
    Compute the Out-of-Bootstrap (OOB) error estimate.
    
    For each bootstrap iteration:
    1. Create bootstrap sample (with replacement)
    2. Train model on bootstrap sample
    3. Evaluate on observations NOT in bootstrap sample
    
    Parameters
    ----------
    X : np.ndarray
        Feature matrix of shape (n_samples, n_features)
    y : np.ndarray
        Target vector of shape (n_samples,)
    model : sklearn estimator
        Model with fit() and predict() methods
    n_bootstrap : int
        Number of bootstrap iterations
    scoring_fn : Callable, optional
        Custom scoring function(y_true, y_pred) -> float
        Default: mean squared error for regression
    random_state : int
        Random seed
    
    Returns
    -------
    dict
        Dictionary containing OOB error estimate and diagnostics
    """
    if scoring_fn is None:
        # Default: MSE for regression
        scoring_fn = lambda y_true, y_pred: np.mean((y_true - y_pred)**2)
    
    rng = np.random.RandomState(random_state)
    n_samples = len(y)
    
    # Track predictions for each observation when it's OOB
    oob_predictions = {i: [] for i in range(n_samples)}
    bootstrap_errors = []
    oob_sizes = []
    
    for b in range(n_bootstrap):
        # Generate bootstrap indices
        bootstrap_indices = rng.choice(n_samples, size=n_samples, replace=True)
        
        # Identify in-bag and out-of-bag observations
        in_bag = set(bootstrap_indices)
        out_of_bag = [i for i in range(n_samples) if i not in in_bag]
        
        if len(out_of_bag) == 0:
            continue  # Rare case: all observations selected
        
        oob_sizes.append(len(out_of_bag))
        
        # Train on bootstrap sample
        X_train = X[bootstrap_indices]
        y_train = y[bootstrap_indices]
        
        model_clone = clone(model)
        model_clone.fit(X_train, y_train)
        
        # Predict on OOB observations
        X_oob = X[out_of_bag]
        y_oob_true = y[out_of_bag]
        y_oob_pred = model_clone.predict(X_oob)
        
        # Store predictions for aggregation
        for i, obs_idx in enumerate(out_of_bag):
            oob_predictions[obs_idx].append(y_oob_pred[i])
        
        # Compute error for this bootstrap iteration
        bootstrap_errors.append(scoring_fn(y_oob_true, y_oob_pred))
    
    # Aggregate OOB predictions per observation
    aggregated_predictions = []
    aggregated_true = []
    prediction_counts = []
    
    for i in range(n_samples):
        if len(oob_predictions[i]) > 0:
            aggregated_predictions.append(np.mean(oob_predictions[i]))
            aggregated_true.append(y[i])
            prediction_counts.append(len(oob_predictions[i]))
    
    # Final OOB error estimate (using aggregated predictions)
    oob_error = scoring_fn(
        np.array(aggregated_true), 
        np.array(aggregated_predictions)
    )
    
    return {
        'oob_error': oob_error,
        'oob_error_std': np.std(bootstrap_errors),
        'mean_bootstrap_error': np.mean(bootstrap_errors),
        'mean_oob_size': np.mean(oob_sizes),
        'mean_oob_fraction': np.mean(oob_sizes) / n_samples,
        'mean_predictions_per_obs': np.mean(prediction_counts),
        'observations_with_predictions': len(aggregated_true),
    }
 
 
# Example usage
if __name__ == "__main__":
    from sklearn.linear_model import Ridge
    from sklearn.datasets import make_regression
    
    # Generate sample data
    X, y = make_regression(n_samples=200, n_features=10, noise=10, random_state=42)
    
    # Compute OOB error
    model = Ridge(alpha=1.0)
    results = out_of_bootstrap_error(X, y, model, n_bootstrap=500)
    
    print("Out-of-Bootstrap Error Estimation")
    print("=" * 50)
    print(f"OOB Error (MSE): {results['oob_error']:.4f}")
    print(f"OOB Error Std: {results['oob_error_std']:.4f}")
    print(f"Mean OOB fraction: {results['mean_oob_fraction']:.4f}")
    print(f"Mean predictions per observation: {results['mean_predictions_per_obs']:.1f}")
    print(f"Observations covered: {results['observations_with_predictions']}/{len(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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174

import numpy as np
from typing import List, Tuple
 
def stratified_bootstrap(data: np.ndarray, 
                          strata: np.ndarray,
                          statistic_fn,
                          n_bootstrap: int = 1000,
                          random_state: int = 42) -> np.ndarray:
    """
    Perform stratified bootstrap resampling.
    
    Resamples WITHIN each stratum separately, preserving
    the stratification structure.
    
    Parameters
    ----------
    data : np.ndarray
        Data array of shape (n_samples,) or (n_samples, n_features)
    strata : np.ndarray
        Stratum labels for each observation, shape (n_samples,)
    statistic_fn : Callable
        Function computing statistic from resampled data
    n_bootstrap : int
        Number of bootstrap iterations
    random_state : int
        Random seed
    
    Returns
    -------
    np.ndarray
        Array of B bootstrap statistics
    """
    rng = np.random.RandomState(random_state)
    unique_strata = np.unique(strata)
    
    # Indices for each stratum
    stratum_indices = {s: np.where(strata == s)[0] for s in unique_strata}
    
    bootstrap_stats = np.zeros(n_bootstrap)
    
    for b in range(n_bootstrap):
        resampled_indices = []
        
        # Resample within each stratum
        for s in unique_strata:
            indices = stratum_indices[s]
            n_stratum = len(indices)
            
            # Sample with replacement within stratum
            selected = rng.choice(indices, size=n_stratum, replace=True)
            resampled_indices.extend(selected)
        
        # Extract resampled data
        resampled_data = data[resampled_indices]
        
        # Compute statistic
        bootstrap_stats[b] = statistic_fn(resampled_data)
    
    return bootstrap_stats
 
 
def cluster_bootstrap(data: np.ndarray,
                       clusters: np.ndarray,
                       statistic_fn,
                       n_bootstrap: int = 1000,
                       random_state: int = 42) -> np.ndarray:
    """
    Perform cluster bootstrap resampling.
    
    Resamples CLUSTERS with replacement, then includes all
    observations within selected clusters.
    
    Parameters
    ----------
    data : np.ndarray
        Data array of shape (n_samples,) or (n_samples, n_features)
    clusters : np.ndarray
        Cluster labels for each observation, shape (n_samples,)
    statistic_fn : Callable
        Function computing statistic from resampled data
    n_bootstrap : int
        Number of bootstrap iterations
    random_state : int
        Random seed
    
    Returns
    -------
    np.ndarray
        Array of B bootstrap statistics
    """
    rng = np.random.RandomState(random_state)
    unique_clusters = np.unique(clusters)
    n_clusters = len(unique_clusters)
    
    # Indices for each cluster
    cluster_indices = {c: np.where(clusters == c)[0] for c in unique_clusters}
    
    bootstrap_stats = np.zeros(n_bootstrap)
    
    for b in range(n_bootstrap):
        # Sample clusters with replacement
        selected_clusters = rng.choice(unique_clusters, size=n_clusters, replace=True)
        
        # Collect all observations from selected clusters
        resampled_indices = []
        for c in selected_clusters:
            resampled_indices.extend(cluster_indices[c])
        
        # Extract resampled data
        resampled_data = data[resampled_indices]
        
        # Compute statistic
        bootstrap_stats[b] = statistic_fn(resampled_data)
    
    return bootstrap_stats
 
 
# Example: Time series block bootstrap
def block_bootstrap(data: np.ndarray,
                     block_size: int,
                     statistic_fn,
                     n_bootstrap: int = 1000,
                     random_state: int = 42) -> np.ndarray:
    """
    Perform block bootstrap for time series data.
    
    Resamples contiguous blocks to preserve temporal dependence.
    Uses the moving block bootstrap (MBB) approach.
    
    Parameters
    ----------
    data : np.ndarray
        Time series of shape (n_timepoints,)
    block_size : int
        Size of each block (should reflect autocorrelation range)
    statistic_fn : Callable
        Function computing statistic from resampled series
    n_bootstrap : int
        Number of bootstrap iterations
    random_state : int
        Random seed
    
    Returns
    -------
    np.ndarray
        Array of B bootstrap statistics
    """
    rng = np.random.RandomState(random_state)
    n = len(data)
    
    # Number of blocks needed to cover at least n observations
    n_blocks_needed = int(np.ceil(n / block_size))
    
    # Possible starting positions for blocks
    max_start = n - block_size
    
    bootstrap_stats = np.zeros(n_bootstrap)
    
    for b in range(n_bootstrap):
        # Sample block starting positions
        starts = rng.randint(0, max_start + 1, size=n_blocks_needed)
        
        # Construct resampled series from blocks
        resampled = []
        for start in starts:
            resampled.extend(data[start:start + block_size])
        
        # Trim to exactly n observations
        resampled = np.array(resampled[:n])
        
        # Compute statistic
        bootstrap_stats[b] = statistic_fn(resampled)
    
    return bootstrap_stats