Un intervalo de confianza para cualquier cosa

En el capítulo 10 aprendimos a ponerle un intervalo al promedio. Y quedó algo colgando que igual ni notaste 🥾

Ese intervalo salió de una fórmula, y esa fórmula solo existe para el promedio. Para la mediana no la hay. Ni para el percentil 90, ni para la razón entre dos cosas, ni para casi nada de lo que te van a preguntar en una reunión de verdad.

import numpy as np
import pandas as pd
from scipy import stats

URL = 'https://missyera.com/static/datasets/ventas-miss-yera.csv'


def carga_limpia(url):
    """La misma del capítulo 2."""
    v = pd.read_csv(url).drop_duplicates()
    v['ciudad'] = (v['ciudad'].str.strip().str.lower()
                   .str.normalize('NFKD')
                   .str.encode('ascii', 'ignore').str.decode('utf-8'))
    v['monto'] = pd.to_numeric(v['monto'].str.replace(',', '.'))
    for col in ['fecha', 'fecha_ultima_compra']:
        f = pd.to_datetime(v[col], format='%Y-%m-%d', errors='coerce')
        falta = f.isna() & v[col].notna()
        f[falta] = pd.to_datetime(v.loc[falta, col], format='%d/%m/%Y',
                                  errors='coerce')
        v[col] = f
    return v


v = carga_limpia(URL)
m = v['monto'].dropna()

ee = m.std() / np.sqrt(len(m))
bajo, alto = stats.t.interval(0.95, len(m) - 1, m.mean(), ee)
print('promedio %.2f, y su intervalo va de %.2f a %.2f' % (m.mean(), bajo, alto))
print('mediana  %.2f, y su intervalo va de ...' % m.median())

promedio 803.76, y su intervalo va de 776.84 a 830.68
mediana  533.63, y su intervalo va de ...

Ahí está el problema, escrito con puntos suspensivos 😅

Y fíjate que la mediana es la que más falta hace aquí: 533,63 contra 803,76 de promedio, porque los montos tienen cola y el promedio se va para arriba. Eso ya lo vimos en el capítulo 3. O sea que el número que mejor describe estas ventas es justo el que no sabemos acompañar.

Antes de seguir, piénsalo tú: ¿de dónde saldría el intervalo si pudieras repetir el experimento mil veces? Guárdate esa respuesta, que es todo el capítulo 🐣

La idea, que es de las más bonitas que hay

Un intervalo de confianza contesta "cuánto se movería este número si volviera a tomar la muestra". El problema es que no puedes volver a tomarla: los 3.000 pedidos son los que hay.

El truco del bootstrap es descarado y funciona: si no puedes repetir el experimento, repite la muestra. Se saca de tus propios datos una muestra del mismo tamaño, con reemplazo, o sea permitiendo que una fila salga dos veces y otra ninguna. Eso da una muestra parecida pero no idéntica. Y se hace cinco mil veces.

Lo de "con reemplazo" no es un detalle técnico, es la idea entera. Sin reemplazo estarías barajando las mismas 3.000 filas y saldría siempre el mismo número.

def bootstrap(datos, resumen, veces=5000, semilla=7):
    """El intervalo de lo que sea, sin una sola formula."""
    generador = np.random.default_rng(semilla)
    datos = np.asarray(datos)
    repes = [resumen(generador.choice(datos, len(datos), replace=True))
             for _ in range(veces)]
    return np.percentile(repes, [2.5, 97.5])


print('la formula del promedio dice: de %.2f a %.2f' % (bajo, alto))
print('y el bootstrap dice         : de %.2f a %.2f'
      % tuple(bootstrap(m.values, np.mean)))

la formula del promedio dice: de 776.84 a 830.68
y el bootstrap dice         : de 776.94 a 830.85

Iguales hasta el primer decimal 🎯

Y ese es el orden que hay que respetar siempre: primero se comprueba el método donde ya sabes la respuesta, y después se usa donde no. Si el bootstrap no hubiera reproducido la fórmula, no habría por qué creerle en lo que viene ahora.

Fíjate también en lo que la función NO tiene: ninguna distribución, ninguna t de Student, ningún supuesto de normalidad. Solo remuestrear y mirar los percentiles 2,5 y 97,5, que dejan el 95% en medio.

Y ahora la mediana, que era lo que no se podía

print('mediana: %.2f' % m.median())
print('su intervalo del 95%%: de %.2f a %.2f'
      % tuple(bootstrap(m.values, np.median)))

mediana: 533.63
su intervalo del 95%: de 510.11 a 552.60

Ya está. Cambiando np.mean por np.median 💜

Eso es lo que hace al bootstrap tan útil: la función no sabe ni le importa qué le pasas. Cualquier cosa que sepas calcular sobre una muestra, sabes acompañarla de su intervalo.

Algo que ninguna fórmula te da

Ahora una de verdad, de las que te piden. Cuánto vale un punto de satisfacción, medido en soles de mediana:

sat = v[['monto', 'satisfaccion']].dropna()
observada = sat['monto'].median() / sat['satisfaccion'].median()

generador = np.random.default_rng(7)
filas = np.arange(len(sat))
repes = []
for _ in range(2000):
    t = sat.iloc[generador.choice(filas, len(filas), replace=True)]
    repes.append(t['monto'].median() / t['satisfaccion'].median())

print('filas con las dos columnas: %d' % len(sat))
print('soles de mediana por punto de satisfaccion: %.4f' % observada)
print('intervalo del 95%%: de %.4f a %.4f' % tuple(np.percentile(repes, [2.5, 97.5])))

filas con las dos columnas: 2769
soles de mediana por punto de satisfaccion: 177.1533
intervalo del 95%: de 169.0233 a 184.0771

Una razón entre dos medianas, con su intervalo. Búscale la fórmula a eso 🙃

Y mira una cosa importante del código: se remuestrean las filas enteras, no cada columna por su lado. Si remuestrearas monto y satisfacción por separado romperías la relación entre las dos, que es justo lo que estás midiendo. Cuando el número junta dos columnas, la unidad que se remuestrea es la fila.

Comparar dos grupos sin ninguna prueba

Y esto es lo que más te va a servir. La pregunta de siempre: ¿vende más caro Web o WhatsApp?

canales = v.dropna(subset=['monto', 'canal'])
web = canales[canales['canal'] == 'Web']['monto'].values
wsp = canales[canales['canal'] == 'WhatsApp']['monto'].values

generador = np.random.default_rng(7)
difs = [np.median(generador.choice(web, len(web), replace=True))
        - np.median(generador.choice(wsp, len(wsp), replace=True))
        for _ in range(3000)]

print('Web %d pedidos, WhatsApp %d' % (len(web), len(wsp)))
print('la mediana de Web supera a la de WhatsApp en %.2f soles'
      % (np.median(web) - np.median(wsp)))
print('intervalo del 95%%: de %.2f a %.2f' % tuple(np.percentile(difs, [2.5, 97.5])))

Web 792 pedidos, WhatsApp 719
la mediana de Web supera a la de WhatsApp en 12.17 soles
intervalo del 95%: de -45.08 a 64.80

Doce soles de diferencia, y el intervalo va de -45,08 a 64,80 🫠

Cruza el cero. Con estos datos, la mediana de Web podría estar 45 soles por debajo de la de WhatsApp igual de bien que 65 por encima. La diferencia que viste no aguanta.

Y esto es lo que hace el bootstrap tan honesto: no te da un sí o un no como una prueba de hipótesis, te da el rango de lo que es compatible con tus datos, y ahí se ve solo que doce soles no significan nada. Es la misma lección del capítulo 14, contada por otro camino.

Dónde se rompe, que es lo que casi nadie cuenta

Te lo enseño porque un método que parece servir para todo da miedo, y con razón. Pídele al bootstrap un intervalo para el pedido más caro:

extremos = [np.max(np.random.default_rng(s).choice(m.values, len(m), replace=True))
            for s in range(3000)]
print('el pedido mas caro que existe: %.2f' % m.max())
print('intervalo bootstrap del maximo: de %.2f a %.2f'
      % tuple(np.percentile(extremos, [2.5, 97.5])))
print('repeticiones que devuelven ese mismo maximo: %d de 3000'
      % int(np.sum(np.array(extremos) == m.max())))

el pedido mas caro que existe: 4236.71
intervalo bootstrap del maximo: de 3751.02 a 4236.71
repeticiones que devuelven ese mismo maximo: 1895 de 3000

Mira el techo del intervalo 😳

Es exactamente 4.236,71, que es el máximo observado. Y no es casualidad: el bootstrap remuestrea de lo que hay, así que nunca puede devolver un número mayor que el más grande que viste. De hecho 1.895 de las 3.000 repeticiones devuelven ese mismísimo valor, porque en una muestra de 3.000 con reemplazo es rarísimo que el más caro no salga.

Y ese intervalo es falso. El pedido más caro que tu negocio puede recibir el mes que viene no está acotado por el más caro que recibió hasta hoy. La pregunta es legítima; el bootstrap no la sabe contestar.

La regla, y sirve para más cosas que esta: el bootstrap funciona para números que dependen de toda la muestra, y falla para los que dependen de un solo dato del borde. El promedio y la mediana usan todo. El máximo, el mínimo y el percentil 99 con pocos datos dependen de una fila, y ahí no hay remuestreo que valga.

El error que sale al escribirlo mal

Y este error es bonito porque el que revienta es justo el concepto:

generador = np.random.default_rng(7)
generador.choice(wsp, len(canales), replace=False)

ValueError: Cannot take a larger sample than population when replace is False

No puedes tomar una muestra más grande que la población sin reemplazo. Pues eso es el bootstrap dicho al revés: solo se puede remuestrear porque se permite repetir filas. Quítale el reemplazo y el método deja de existir 🔑

Y en tu trabajo, ¿cuándo lo sacas?

La última fila es la que hay que saber decir en voz alta en una reunión, y es justo la que este capítulo mide 🔍

Si quieres el original, el bootstrap lo publicó Bradley Efron en Bootstrap Methods: Another Look at the Jackknife (1979), y scipy trae una versión con más opciones en scipy.stats.bootstrap, que hace lo mismo que nuestras tres líneas y además corrige el sesgo 📄

repes = []
for _ in range(5000):
    muestra = generador.choice(montos, len(montos), replace=False)
    repes.append(np.median(muestra))

print(np.percentile(repes, [2.5, 97.5]))   # [533.63, 533.63]

Ejercicios

Seis. El 1 es el que quiero que hagas sí o sí, porque es la trampa de arriba ejecutada 💛

generador = np.random.default_rng(7)
mal = [np.median(generador.choice(m.values, len(m), replace=False))
       for _ in range(200)]
print('valores distintos que salieron:', len(set(mal)))
print('intervalo:', np.percentile(mal, [2.5, 97.5]))

p90 = lambda x: np.percentile(x, 90)
print('mediana   :', np.round(bootstrap(m.values, np.median), 2))
print('percentil 90:', np.round(bootstrap(m.values, p90), 2))

from scipy.stats import mannwhitneyu
print('valor p:', round(mannwhitneyu(web, wsp).pvalue, 4))

for veces in [100, 1000, 10000]:
    izq, der = np.percentile(bootstrap(v['monto'], np.mean, veces=veces),
                             [2.5, 97.5])
    print('con %6d remuestreos: de %7.2f a %7.2f, ancho %6.2f'
          % (veces, izq, der, der - izq))

con    100 remuestreos: de  779.02 a  823.08, ancho  44.07
con   1000 remuestreos: de  780.37 a  828.75, ancho  48.38
con  10000 remuestreos: de  778.36 a  829.61, ancho  51.25

for ciudad, t in v.groupby('ciudad'):
    izq, der = np.percentile(bootstrap(t['monto'], np.mean), [2.5, 97.5])
    print('%-10s n=%4d  promedio %7.2f  de %7.2f a %7.2f'
          % (ciudad, len(t), t['monto'].mean(), izq, der))

arequipa   n= 546  promedio  801.92  de  744.09 a  867.10
chiclayo   n= 510  promedio  813.91  de  750.86 a  876.31
cusco      n= 496  promedio  780.70  de  722.27 a  842.35
lima       n= 471  promedio  792.74  de  730.13 a  855.81
piura      n= 506  promedio  792.16  de  734.38 a  854.71
trujillo   n= 471  promedio  842.69  de  771.11 a  916.17

Lo que te llevas

• 🥾 Bootstrap es remuestrear tus datos con reemplazo miles de veces y mirar los percentiles 2,5 y 97,5 de lo que salga.
• ✅ Primero se comprueba donde ya hay fórmula: 776,94 a 830,85 contra 776,84 a 830,68. Después se usa donde no la hay.
• 🔁 Sirve para la mediana, para una razón, para una diferencia y para cualquier cosa que sepas calcular sobre una muestra.
• 🧩 Cuando el número junta dos columnas, se remuestrean filas enteras.
• 🚫 Y falla para lo que depende de un dato del borde: el intervalo del máximo tiene por techo el máximo que ya viste, y eso es mentira.

En el capítulo 12 el valor p se calculó barajando etiquetas, que es este mismo truco con otra cara: cuando no hay fórmula, se simula. Y el vocabulario suelto está en el glosario de IA 📖

Que tengas lindo día! 🌸