Phương pháp Bootstrap cho Độ lệch chuẩn

Bootstrap: Cuộc cách mạng thống kê thời đại máy tính

Tái lấy mẫu bootstrap là kỹ thuật thống kê mạnh mẽ ước lượng phân phối lấy mẫu của bất kỳ thống kê nào bằng cách lấy mẫu lặp đi lặp lại từ dữ liệu quan sát. Được Bradley Efron giới thiệu năm 1979, nó đã cách mạng hóa suy luận thống kê bằng cách cho phép phân tích các thống kê phức tạp mà không cần dựa vào công thức toán học hay giả định phân phối.

Ý tưởng cốt lõi của bootstrap cực kỳ đơn giản: mẫu của bạn là ước lượng tốt nhất cho tổng thể. Bằng cách tái lấy mẫu từ mẫu của bạn (có hoàn lại), bạn mô phỏng điều sẽ xảy ra nếu bạn có thể lấy mẫu lặp đi lặp lại từ tổng thể. Phương pháp này đặc biệt có giá trị cho độ lệch chuẩn, nơi công thức khoảng tin cậy truyền thống giả định tính chuẩn—một giả định thường không đúng trong thực tế.

Bootstrap đã trở nên thiết yếu trong khoa học dữ liệu hiện đại vì nó hoạt động với bất kỳ thống kê nào (trung vị, tương quan, hệ số hồi quy, trọng số mạng nơ-ron) và không đưa ra giả định về phân phối cơ bản của dữ liệu.

Tại sao Bootstrap cho Độ lệch chuẩn?

Khoảng tin cậy truyền thống cho độ lệch chuẩn giả định dữ liệu đến từ phân phối chuẩn. Khi giả định này thất bại (điều phổ biến), các khoảng này có thể sai lệch nghiêm trọng. Bootstrap cung cấp phương án thay thế không phụ thuộc phân phối.

Ưu điểm chính của bootstrap cho độ lệch chuẩn:

• Không giả định phân phối: Hoạt động tốt như nhau với dữ liệu chuẩn, lệch hoặc đuôi nặng
• Hiệu suất mẫu nhỏ: Thường chính xác hơn phương pháp tham số với n < 30
• Xử lý thống kê phức tạp: Cùng phương pháp hoạt động cho SD cắt ngọn, MAD hoặc đại lượng biến thiên tùy chỉnh
• Trực quan: Phân phối bootstrap cho bạn thấy điều đang xảy ra, không chỉ con số cuối cùng

Quy trình Bootstrap

Thuật toán bootstrap cực kỳ đơn giản. Từ mẫu gốc n quan sát:

Bao nhiêu mẫu bootstrap? B = 1.000 thường đủ cho ước lượng sơ bộ và kiểm định giả thuyết. Cho khoảng tin cậy, B = 10.000 cho phân vị ổn định. Cho khoảng BCa chất lượng công bố, B = 15.000+ được khuyến nghị.

Các phương pháp khoảng tin cậy Bootstrap

Có nhiều phương pháp xây dựng khoảng tin cậy từ mẫu bootstrap, mỗi phương pháp có ưu nhược điểm:

1. Phương pháp phân vị (Đơn giản nhất)

Cách tiếp cận trực quan nhất: lấy phân vị của phân phối bootstrap trực tiếp.

Với 10.000 mẫu bootstrap, đây là giá trị thứ 250 và 9.750 theo thứ tự. Đơn giản nhưng có thể chệch khi phân phối bootstrap lệch.

2. Bootstrap cơ bản (Trục)

Sử dụng mối quan hệ giữa thống kê mẫu và thống kê bootstrap:

Trong đó θ̂ là SD mẫu gốc. Phương pháp này “phản chiếu” khoảng phân vị quanh ước lượng mẫu.

3. BCa (Hiệu chỉnh chệch và Gia tốc)

Tiêu chuẩn vàng về độ chính xác. BCa điều chỉnh cho cả chệch trong phân phối bootstrap và gia tốc (cách sai số chuẩn thay đổi theo giá trị tham số). Phức tạp hơn để tính nhưng cho khoảng chính xác bậc hai.

Ví dụ có lời giải: Dữ liệu không chuẩn

Xem xét 15 phép đo thời gian phản ứng (ms): 245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312. Dữ liệu này lệch phải (một số phản ứng rất chậm).

CI bootstrap bất đối xứng (rộng hơn phía cao), phản ánh tính lệch phải của dữ liệu. CI chi bình phương không nắm bắt tính bất đối xứng này.

Triển khai Python

Triển khai bootstrap hoàn chỉnh với nhiều phương pháp CI:

import numpy as np
from scipy import stats

def bootstrap_sd_ci(data, n_bootstrap=10000, ci=0.95, method='percentile'):
    """
    Bootstrap confidence interval for standard deviation.

    Parameters:
    -----------
    data : array-like - Original sample
    n_bootstrap : int - Number of bootstrap samples
    ci : float - Confidence level (e.g., 0.95)
    method : str - 'percentile', 'basic', or 'bca'

    Returns:
    --------
    tuple : (lower_bound, upper_bound, bootstrap_sds)
    """
    data = np.array(data)
    n = len(data)
    original_sd = np.std(data, ddof=1)

    # Generate bootstrap samples and calculate SDs
    bootstrap_sds = np.array([
        np.std(np.random.choice(data, size=n, replace=True), ddof=1)
        for _ in range(n_bootstrap)
    ])

    alpha = 1 - ci

    if method == 'percentile':
        lower = np.percentile(bootstrap_sds, 100 * alpha/2)
        upper = np.percentile(bootstrap_sds, 100 * (1 - alpha/2))

    elif method == 'basic':
        lower = 2*original_sd - np.percentile(bootstrap_sds, 100*(1-alpha/2))
        upper = 2*original_sd - np.percentile(bootstrap_sds, 100*alpha/2)

    elif method == 'bca':
        # Bias correction
        prop_less = np.mean(bootstrap_sds < original_sd)
        z0 = stats.norm.ppf(prop_less)

        # Acceleration (jackknife estimate)
        jackknife_sds = np.array([
            np.std(np.delete(data, i), ddof=1) for i in range(n)
        ])
        jack_mean = jackknife_sds.mean()
        a = np.sum((jack_mean - jackknife_sds)**3) / \
            (6 * np.sum((jack_mean - jackknife_sds)**2)**1.5)

        # Adjusted percentiles
        z_alpha = stats.norm.ppf([alpha/2, 1-alpha/2])
        adj_percentiles = stats.norm.cdf(
            z0 + (z0 + z_alpha) / (1 - a*(z0 + z_alpha))
        ) * 100
        lower = np.percentile(bootstrap_sds, adj_percentiles[0])
        upper = np.percentile(bootstrap_sds, adj_percentiles[1])

    return lower, upper, bootstrap_sds

# Example usage
response_times = [245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312]

for method in ['percentile', 'basic', 'bca']:
    lower, upper, _ = bootstrap_sd_ci(response_times, method=method)
    print(f"{method.upper():12s} 95% CI: [{lower:.1f}, {upper:.1f}]")