שיטות בוטסטרפ לסטיית תקן

בוטסטרפ: המהפכה הסטטיסטית של עידן המחשב

דגימת בוטסטרפ היא טכניקה סטטיסטית עוצמתית שמאמידה את התפלגות הדגימה של כל סטטיסטי באמצעות דגימה חוזרת מהנתונים הנצפים. השיטה הוצגה על ידי בראדלי אפרון ב-1979 וחוללה מהפכה בהסקה סטטיסטית, בכך שאיפשרה ניתוח של סטטיסטים מורכבים ללא הסתמכות על נוסחאות מתמטיות או הנחות על התפלגות הנתונים.

התובנה המרכזית מאחורי בוטסטרפ פשוטה ואלגנטית: המדגם שלכם הוא האומדן הטוב ביותר של האוכלוסייה. על ידי דגימה חוזרת מהמדגם (עם החזרה), אתם מדמים את מה שהיה קורה אילו יכולתם לדגום שוב ושוב מהאוכלוסייה. גישה זו חשובה במיוחד עבור סטיית תקן, כאשר נוסחאות מסורתיות לרווחי סמך מניחות התפלגות נורמלית — הנחה שלעתים קרובות אינה מתקיימת בפועל.

בוטסטרפ הפך לכלי חיוני במדעי הנתונים המודרניים מכיוון שהוא עובד עם כל סטטיסטי (חציון, מתאם, מקדמי רגרסיה, משקולות רשתות עצביות) ואינו מניח דבר על התפלגות הנתונים.

למה להשתמש בבוטסטרפ לסטיית תקן?

רווחי סמך מסורתיים לסטיית תקן מניחים שהנתונים מגיעים מהתפלגות נורמלית. כאשר הנחה זו אינה מתקיימת (מצב שכיח), רווחים אלה עלולים להיות לא מדויקים באופן קיצוני. בוטסטרפ מספק חלופה שאינה תלויה בהתפלגות.

יתרונות מרכזיים של בוטסטרפ לסטיית תקן:

• ללא הנחות על התפלגות: עובד באותה מידה עם נתונים נורמליים, א-סימטריים או בעלי זנבות כבדים
• ביצועים טובים במדגמים קטנים: לעתים קרובות מדויק יותר משיטות פרמטריות כאשר n < 30
• מתאים לסטטיסטים מורכבים: אותה גישה עובדת עבור סטיית תקן חתוכה, MAD או מדדי פיזור מותאמים אישית
• תובנה חזותית: התפלגות הבוטסטרפ מראה לכם מה קורה, לא רק מספרים סופיים

תהליך הבוטסטרפ

אלגוריתם הבוטסטרפ פשוט להפליא. מתוך המדגם המקורי של n תצפיות:

כמה מדגמי בוטסטרפ? B = 1,000 מספיק בדרך כלל לאומדנים גסים ובדיקות השערות. לרווחי סמך, B = 10,000 מספק אחוזונים יציבים. לרווחי BCa באיכות פרסום, מומלץ B = 15,000 ומעלה.

שיטות לרווחי סמך מבוססי בוטסטרפ

קיימות מספר שיטות לבניית רווחי סמך ממדגמי בוטסטרפ, כל אחת עם יתרונות וחסרונות:

1. שיטת האחוזון (הפשוטה ביותר)

הגישה האינטואיטיבית ביותר: קחו את האחוזונים של התפלגות הבוטסטרפ ישירות.

עבור 10,000 מדגמי בוטסטרפ, אלה הערכים ה-250 וה-9,750 בסדר עולה. פשוט אך עלול להיות מוטה כאשר התפלגות הבוטסטרפ א-סימטרית.

2. בוטסטרפ בסיסי (ציר)

משתמש ביחס בין הסטטיסטי של המדגם לסטטיסטים של הבוטסטרפ:

כאשר θ̂ היא סטיית התקן של המדגם המקורי. שיטה זו “משקפת” את רווח האחוזון סביב אומדן המדגם.

3. BCa (מתוקנת הטיה ומואצת)

תקן הזהב לדיוק. BCa מתקנת הן הטיה בהתפלגות הבוטסטרפ והן האצה (כיצד שגיאת התקן משתנה עם ערך הפרמטר). מורכבת יותר לחישוב אך מספקת רווחים מדויקים מסדר שני.

דוגמה מפורטת: נתונים לא נורמליים

נניח 15 מדידות של זמני תגובה (באלפיות שנייה): 245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312. נתונים אלה מוטים ימינה (חלק מזמני התגובה ארוכים במיוחד).

רווח הסמך של הבוטסטרפ א-סימטרי (רחב יותר בצד הגבוה), מה שמשקף את האופי הא-סימטרי ימינה של הנתונים. רווח הסמך של חי-בריבוע אינו לוכד א-סימטריה זו.

מימוש ב-Python

מימוש בוטסטרפ מלא עם מספר שיטות לרווחי סמך:

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}]")