Métodos bootstrap para la desviación estándar

Bootstrap: la revolución estadística de la era informática

El remuestreo bootstrap es una técnica estadística poderosa que estima la distribución muestral de cualquier estadístico mediante el remuestreo repetido de los datos observados. Introducido por Bradley Efron en 1979, revolucionó la inferencia estadística al permitir el análisis de estadísticos complejos sin depender de fórmulas matemáticas o supuestos distribucionales.

La idea clave detrás del bootstrap es elegantemente simple: su muestra es la mejor estimación de la población. Al remuestrear de su muestra (con reemplazo), se simula lo que ocurriría si pudiera muestrear repetidamente de la población. Este enfoque es particularmente valioso para la desviación estándar, donde las fórmulas tradicionales de intervalos de confianza asumen normalidad, un supuesto que frecuentemente no se cumple en la práctica.

El bootstrap se ha vuelto esencial en la ciencia de datos moderna porque funciona con cualquier estadístico (mediana, correlación, coeficientes de regresión, pesos de redes neuronales) y no hace supuestos sobre la distribución subyacente de sus datos.

¿Por qué bootstrap para la desviación estándar?

Los intervalos de confianza tradicionales para la desviación estándar asumen que sus datos provienen de una distribución normal. Cuando este supuesto falla (lo cual es habitual), estos intervalos pueden ser muy imprecisos. El bootstrap proporciona una alternativa libre de distribución.

Ventajas clave del bootstrap para la desviación estándar:

• Sin supuestos distribucionales: Funciona igualmente bien con datos normales, sesgados o de colas pesadas
• Rendimiento con muestras pequeñas: Frecuentemente más preciso que los métodos paramétricos con n < 30
• Manejo de estadísticos complejos: El mismo enfoque funciona para DE recortada, MAD o medidas de variabilidad personalizadas
• Visión visual: La distribución bootstrap muestra lo que está ocurriendo, no solo los números finales

El procedimiento bootstrap

El algoritmo bootstrap es notablemente sencillo. A partir de su muestra original de n observaciones:

¿Cuántas muestras bootstrap? B = 1,000 suele ser suficiente para estimaciones aproximadas y pruebas de hipótesis. Para intervalos de confianza, B = 10,000 proporciona percentiles estables. Para intervalos BCa de calidad de publicación, se recomiendan B = 15,000 o más.

Métodos de intervalos de confianza bootstrap

Existen varios métodos para construir intervalos de confianza a partir de muestras bootstrap, cada uno con sus ventajas y desventajas:

1. Método de percentiles (el más simple)

El enfoque más intuitivo: tome los percentiles de la distribución bootstrap directamente.

Para 10,000 muestras bootstrap, estos son los valores ordenados en las posiciones 250 y 9,750. Simple pero puede ser sesgado cuando la distribución bootstrap es asimétrica.

2. Bootstrap básico (pivotal)

Utiliza la relación entre el estadístico muestral y los estadísticos bootstrap:

Donde θ̂ es la DE muestral original. Esto "refleja" el intervalo de percentiles alrededor de la estimación muestral.

3. BCa (corregido por sesgo y acelerado)

El estándar de referencia en precisión. El BCa ajusta tanto el sesgo en la distribución bootstrap como la aceleración (cómo cambia el error estándar con el valor del parámetro). Más complejo de calcular pero proporciona intervalos con precisión de segundo orden.

Ejemplo resuelto: datos no normales

Considere 15 mediciones de tiempos de respuesta (en ms): 245, 312, 287, 456, 234, 298, 267, 523, 289, 301, 278, 645, 256, 289, 312. Estos datos tienen sesgo a la derecha (algunas respuestas muy lentas).

El IC bootstrap es asimétrico (más amplio en el lado alto), reflejando la naturaleza sesgada a la derecha de los datos. El IC chi-cuadrado no captura esta asimetría.

Implementación en Python

Implementación completa de bootstrap con múltiples métodos de IC:

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