วิธีบูตสแตรปสำหรับส่วนเบี่ยงเบนมาตรฐาน

บูตสแตรป: การปฏิวัติทางสถิติยุคคอมพิวเตอร์

การสุ่มตัวอย่างซ้ำแบบบูตสแตรป เป็นเทคนิคทางสถิติที่ทรงพลังในการประมาณการแจกแจงการสุ่มตัวอย่างของสถิติใดๆ โดยการสุ่มตัวอย่างซ้ำจากข้อมูลที่สังเกตได้ เปิดตัวโดย Bradley Efron ในปี 1979 มันปฏิวัติการอนุมานทางสถิติโดยทำให้การวิเคราะห์สถิติซับซ้อนเป็นไปได้โดยไม่ต้องอาศัยสูตรทางคณิตศาสตร์หรือข้อสมมติเกี่ยวกับการแจกแจง

ข้อมูลเชิงลึกสำคัญเบื้องหลังบูตสแตรปนั้นเรียบง่ายอย่างหรูหรา: ตัวอย่างของคุณคือค่าประมาณที่ดีที่สุดของประชากร โดยการสุ่มตัวอย่างซ้ำจากตัวอย่างของคุณ (แบบใส่คืน) คุณจำลองว่าจะเกิดอะไรขึ้นถ้าคุณสามารถสุ่มตัวอย่างจากประชากรซ้ำๆ ได้ วิธีนี้มีคุณค่าเป็นพิเศษสำหรับส่วนเบี่ยงเบนมาตรฐาน ซึ่งสูตรช่วงความเชื่อมั่นแบบดั้งเดิมสมมติความเป็นปกติ ข้อสมมติที่มักล้มเหลวในทางปฏิบัติ

บูตสแตรปกลายเป็นสิ่งจำเป็นในวิทยาศาสตร์ข้อมูลสมัยใหม่เพราะทำงานได้กับสถิติใดก็ได้ (มัธยฐาน สหสัมพันธ์ สัมประสิทธิ์การถดถอย น้ำหนักเครือข่ายประสาท) และไม่มีข้อสมมติเกี่ยวกับการแจกแจงพื้นฐานของข้อมูล

ทำไมต้องบูตสแตรปสำหรับส่วนเบี่ยงเบนมาตรฐาน?

ช่วงความเชื่อมั่นแบบดั้งเดิมสำหรับส่วนเบี่ยงเบนมาตรฐานสมมติว่าข้อมูลมาจากการแจกแจงปกติ เมื่อข้อสมมตินี้ล้มเหลว (ซึ่งเกิดขึ้นบ่อย) ช่วงเหล่านี้อาจไม่แม่นยำอย่างมาก บูตสแตรปให้ทางเลือกที่ไม่ขึ้นกับการแจกแจง

ข้อได้เปรียบสำคัญของบูตสแตรปสำหรับส่วนเบี่ยงเบนมาตรฐาน:

• ไม่มีข้อสมมติเกี่ยวกับการแจกแจง: ทำงานได้ดีเท่ากันกับข้อมูลปกติ เบ้ หรือหางหนัก
• ประสิทธิภาพกับตัวอย่างเล็ก: มักแม่นยำกว่าวิธีแบบพารามิเตอร์เมื่อ n < 30
• จัดการสถิติซับซ้อน: วิธีเดียวกันใช้ได้กับ SD ที่ตัดแต่ง MAD หรือตัววัดความแปรผันที่กำหนดเอง
• ข้อมูลเชิงลึกเชิงภาพ: การแจกแจงบูตสแตรปแสดงให้คุณเห็นว่าเกิดอะไรขึ้น ไม่ใช่แค่ตัวเลขสุดท้าย

ขั้นตอนบูตสแตรป

อัลกอริทึมบูตสแตรปตรงไปตรงมาอย่างน่าทึ่ง จากตัวอย่างเดิมที่มี n ข้อสังเกต:

ต้องการตัวอย่างบูตสแตรปกี่ตัว? B = 1,000 มักเพียงพอสำหรับค่าประมาณคร่าวๆ และการทดสอบสมมติฐาน สำหรับช่วงความเชื่อมั่น B = 10,000 ให้เปอร์เซ็นไทล์ที่เสถียร สำหรับช่วง BCa คุณภาพสิ่งพิมพ์ แนะนำ B = 15,000+

วิธีช่วงความเชื่อมั่นแบบบูตสแตรป

มีหลายวิธีในการสร้างช่วงความเชื่อมั่นจากตัวอย่างบูตสแตรป แต่ละวิธีมีข้อดีข้อเสีย:

1. วิธีเปอร์เซ็นไทล์ (ง่ายที่สุด)

วิธีที่เข้าใจง่ายที่สุด: ใช้เปอร์เซ็นไทล์ของการแจกแจงบูตสแตรปโดยตรง

สำหรับ 10,000 ตัวอย่างบูตสแตรป นี่คือค่าที่เรียงลำดับตำแหน่งที่ 250 และ 9,750 ง่ายแต่อาจมีความเอนเอียงเมื่อการแจกแจงบูตสแตรปเบ้

2. บูตสแตรปพื้นฐาน (แบบจุดหมุน)

ใช้ความสัมพันธ์ระหว่างสถิติตัวอย่างและสถิติบูตสแตรป:

โดยที่ θ̂ คือ SD ตัวอย่างเดิม สิ่งนี้ “สะท้อน” ช่วงเปอร์เซ็นไทล์รอบค่าประมาณตัวอย่าง

3. BCa (แก้ไขความเอนเอียงและเร่ง)

มาตรฐานทองคำสำหรับความแม่นยำ BCa ปรับแก้ทั้งความเอนเอียงในการแจกแจงบูตสแตรปและการเร่ง (ความคลาดเคลื่อนมาตรฐานเปลี่ยนไปอย่างไรตามค่าพารามิเตอร์) ซับซ้อนกว่าในการคำนวณแต่ให้ช่วงที่แม่นยำในระดับที่สอง

ตัวอย่างพร้อมวิธีทำ: ข้อมูลที่ไม่เป็นปกติ

พิจารณา 15 การวัดเวลาตอบสนอง (เป็น ms): 245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312 ข้อมูลนี้เบ้ขวา (บางค่าตอบสนองช้ามาก)

CI บูตสแตรปไม่สมมาตร (กว้างกว่าด้านสูง) สะท้อนลักษณะเบ้ขวาของข้อมูล CI ไคสแควร์ไม่จับความไม่สมมาตรนี้

การนำไปใช้ด้วย Python

การนำบูตสแตรปไปใช้อย่างสมบูรณ์พร้อมหลายวิธี 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}]")