La campana, y cuándo no aplica

La campana de Gauss es la figura más famosa de la estadística, y también la más aplicada donde no toca 🔔

Importa por una razón muy concreta: media docena de herramientas del libro la dan por supuesta. La regla del 68%, los intervalos de confianza clásicos, la prueba t, la detección de atípicos por desviaciones. Si tus datos no son campana, esas herramientas no se rompen con un error rojo: te devuelven un número equivocado con toda la calma 😶

Así que vamos a aprender a comprobarlo. Y de paso vamos a ver que la forma "oficial" de comprobarlo, la prueba de normalidad, es bastante inútil.

Primero, mirarlo

import pandas as pd
import numpy as np
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)


def histograma(serie, tramos=10):
    """Un histograma de texto, que se ve igual de bien y no necesita gráficos."""
    cuenta, bordes = np.histogram(serie.dropna(), bins=tramos)
    for i in range(tramos):
        print('%8.0f a %8.0f | %4d %s'
              % (bordes[i], bordes[i + 1], cuenta[i], '#' * int(cuenta[i] / 20)))


histograma(v['monto'])

   -2498 a    -1824 |    2
   -1824 a    -1151 |    0
   -1151 a     -477 |    3
    -477 a      196 |  474 #######################
     196 a      869 | 1571 ##############################################################################
     869 a     1543 |  440 ######################
    1543 a     2216 |  288 ##############
    2216 a     2890 |  173 ########
    2890 a     3563 |   43 ##
    3563 a     4237 |    6

Con eso ya sabes que monto no es una campana, y no hizo falta ninguna prueba 👀

Una campana sube y baja simétrica. Esto sube de golpe, tiene su pico entre 196 y 869, y luego baja con una cola larga por la derecha que llega hasta 4.237. Y por la izquierda hay cinco filas sueltas en territorio negativo.

Los dos números que resumen la forma

for col in ['monto', 'unidades', 'satisfaccion', 'descuento']:
    s = v[col].dropna()
    print('%-13s asimetria %7.4f | curtosis %7.4f' % (col, s.skew(), s.kurtosis()))

monto         asimetria  1.3405 | curtosis  1.5458
unidades      asimetria  0.1843 | curtosis -0.3351
satisfaccion  asimetria -0.0149 | curtosis -1.3089
descuento     asimetria -0.0157 | curtosis -1.1604

La asimetría dice hacia dónde cae la cola. Cero es simétrico, positivo es cola a la derecha, negativo a la izquierda. Como regla de campo, entre -0,5 y 0,5 es bastante simétrico y por encima de 1 la cola ya es evidente.

La curtosis (la que devuelve pandas es la de exceso) dice si la campana es puntiaguda o achatada. Cero es como la normal, positivo es más picuda con colas más pesadas, negativo es más plana.

Ahora léelo como si fueran cuatro personas 🙂

• 💰 monto: asimetría 1,3405. Cola a la derecha, confirmado.
• 📦 unidades: 0,1843 y -0,3351. Esto sí parece una campana.
• ⭐ satisfaccion: simétrica pero con curtosis -1,3089, o sea muy plana. Eso no es campana, es una meseta: los cinco valores salen casi igual de veces.
• 🏷️ descuento: lo mismo, plana.

Ahora la prueba oficial, y la decepción

La prueba de Shapiro-Wilk contesta "¿podrían estos datos venir de una normal?". Si el valor p es chico, la respuesta es no.

Se la aplicamos a unidades, la que mejor pinta tenía:

print(stats.shapiro(v['unidades']))

ShapiroResult(statistic=np.float64(0.9869264852139725), pvalue=np.float64(5.180807008091816e-16))

p = 0,00000000000000052 😬

O sea que la prueba dice, con toda la contundencia del mundo, que unidades no es normal. La misma columna que tenía asimetría 0,18 y que en el capítulo 4 cumplía la regla del 68% clavada (68,70%).

¿Quién miente? Ninguno. Mira lo que pasa cuando le doy los mismos datos en trozos de distinto tamaño:

for n in [20, 50, 100, 300, 1000, 3000]:
    muestra = v['unidades'].sample(n, random_state=7)
    print('con %5d filas -> p = %.6f' % (n, stats.shapiro(muestra).pvalue))

con    20 filas -> p = 0.442212
con    50 filas -> p = 0.235741
con   100 filas -> p = 0.045462
con   300 filas -> p = 0.002306
con  1000 filas -> p = 0.000000
con  3000 filas -> p = 0.000000

Ahí está el truco entero 🎩

Los mismos datos. Con 50 filas la prueba dice "normal, adelante". Con 1.000 dice "de ninguna manera". La columna no cambió: cambió cuánta evidencia tienes.

Y no es que Shapiro sea mala. Mira qué pasa con datos que sí vienen de una normal de verdad:

generador = np.random.default_rng(7)

for n in [100, 1000, 3000]:
    print('normal de verdad con %5d datos -> p = %.4f'
          % (n, stats.shapiro(generador.normal(size=n)).pvalue))

normal de verdad con   100 datos -> p = 0.8599
normal de verdad con  1000 datos -> p = 0.4335
normal de verdad con  3000 datos -> p = 0.3883

Con datos normales de verdad, el p se queda alto por muchos datos que le des ✅

O sea que la prueba funciona perfectamente. Lo que pasa es que responde a una pregunta que no es la tuya:

Ningún dato real es exactamente normal. Con suficientes filas siempre se nota. Por eso, en la práctica: mira el histograma y la asimetría, y usa la prueba como acompañamiento, nunca como sentencia 🧭

Guárdate esta idea, porque es la misma que va a salir en el capítulo 12 con todas las pruebas: significativo no quiere decir importante. Con muchos datos, todo sale significativo.

Lo que sí decide: comparar cuantiles

Mi comprobación favorita, porque contesta la pregunta práctica de verdad: "¿me equivoco mucho si trato esto como normal?"

u = v['unidades']
m, s = u.mean(), u.std()

for q in [0.05, 0.25, 0.50, 0.75, 0.95]:
    print('cuantil %.2f -> observado %6.2f | si fuera normal %6.2f'
          % (q, u.quantile(q), stats.norm.ppf(q, m, s)))

cuantil 0.05 -> observado   2.00 | si fuera normal   2.10
cuantil 0.25 -> observado   7.00 | si fuera normal   7.71
cuantil 0.50 -> observado  12.00 | si fuera normal  11.60
cuantil 0.75 -> observado  15.00 | si fuera normal  15.50
cuantil 0.95 -> observado  21.00 | si fuera normal  21.10

Mira qué cerca 😍 En los extremos, que es donde más duele equivocarse, la diferencia es de una décima de unidad.

Entonces, ¿es unidades normal? Formalmente no, y Shapiro tiene razón. ¿Puedo tratarla como normal para decidir cosas de negocio? Perfectamente.

Esa es la respuesta madura, y es la que quiero que te lleves: la normalidad no es un sí o un no, es cuánto me cuesta suponerla 💡

Cuando no hay campana: transformar

Volvamos a monto, que sí estaba torcido de verdad. Hay un truco clásico para colas por la derecha: aplicar el logaritmo.

positivos = v.loc[v['monto'] > 0, 'monto']

print('filas positivas: %d de %d' % (len(positivos), len(v)))
print('asimetria antes:   %.4f' % positivos.skew())
print('asimetria despues: %.4f' % np.log(positivos).skew())

filas positivas: 2979 de 3000
asimetria antes:   1.4278
asimetria despues: -0.1585

De 1,4278 a -0,1585. El logaritmo aplasta la cola y deja algo casi simétrico 🪄

Funciona porque el logaritmo comprime los números grandes mucho más que los chicos: la distancia entre 100 y 1.000 pasa a ser la misma que entre 1.000 y 10.000.

Fíjate en el precio de la operación: hay que tirar 21 filas, las que tienen monto negativo o cero, porque el logaritmo de un número negativo no existe. Nunca transformes sin contar cuántas filas te dejas 🧾

Y hay un regalo escondido en esto:

print('mediana de las positivas:  %.2f' % positivos.median())
print('exp(media del logaritmo):  %.2f' % np.exp(np.log(positivos).mean()))

mediana de las positivas:  536.65
exp(media del logaritmo):  538.39

Casi el mismo número. Lo que hiciste al promediar en escala logarítmica y volver fue calcular la media geométrica, la del capítulo 3. Y en una distribución con cola por la derecha, la media geométrica cae encima de la mediana 🎯

Lo que cuesta transformar es la interpretación. Después de la transformación tus resultados están en "log de soles", que no significa nada para nadie. Hay que deshacerla para contar el resultado, y ahí es donde se cometen errores. Yo transformo solo cuando la herramienta lo necesita de verdad.

El error del capítulo

Vamos a intentar partir la satisfacción en diez tramos iguales, que es lo que harías para ver su forma con más detalle:

pd.qcut(v['satisfaccion'], 10)

ValueError: Bin edges must be unique: Index([1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 5.0], dtype='float64', name='satisfaccion').
You can drop duplicate edges by setting the 'duplicates' kwarg

Este error es de los buenos, porque te está enseñando algo del dato y no del código 🎓

qcut parte en tramos con la misma cantidad de filas cada uno. Para hacer diez tramos necesita diez cortes distintos, y satisfaccion solo tiene cinco valores posibles. Al buscar el borde del decil 10 y del decil 20 encuentra el mismo 1.0 en los dos, y se planta.

La tentación es hacerle caso al mensaje y poner duplicates='drop'. No lo hagas sin pensar: eso taparía el problema y te dejaría tramos de tamaños muy distintos disfrazados de deciles.

Lo correcto es aceptar lo que el error te está diciendo: esa columna no tiene diez tramos, tiene cinco categorías. Se cuenta con value_counts y se acabó.

Practica 💪

print(v['satisfaccion'].value_counts().sort_index())

satisfaccion
1.0    560
2.0    537
3.0    551
4.0    558
5.0    563
Name: count, dtype: int64

generador = np.random.default_rng(7)
uniforme = pd.Series(generador.integers(1, 6, size=3000))

print('satisfaccion -> curtosis %.4f | asimetria %.4f'
      % (v['satisfaccion'].kurtosis(), v['satisfaccion'].skew()))
print('uniforme 1-5 -> curtosis %.4f | asimetria %.4f'
      % (uniforme.kurtosis(), uniforme.skew()))

satisfaccion -> curtosis -1.3089 | asimetria -0.0149
uniforme 1-5 -> curtosis -1.2944 | asimetria -0.0219

for seg in ['Bodega', 'Mayorista']:
    s = v.loc[v['segmento'] == seg, 'monto']
    print('%-10s asimetria %7.4f | curtosis %7.4f | minimo %8.2f'
          % (seg, s.skew(), s.kurtosis(), s.min()))

Bodega     asimetria -0.7752 | curtosis  3.4590 | minimo  -224.61
Mayorista  asimetria -0.4560 | curtosis  3.0510 | minimo -2497.72

m, s = v['monto'].mean(), v['monto'].std()

for q in [0.05, 0.25, 0.50, 0.75, 0.95, 0.99]:
    obs = v['monto'].quantile(q)
    teo = stats.norm.ppf(q, m, s)
    print('cuantil %.2f -> observado %8.2f | normal %8.2f | error %+8.2f'
          % (q, obs, teo, obs - teo))

cuantil 0.05 -> observado   119.03 | normal  -433.15 | error  +552.18
cuantil 0.25 -> observado   252.69 | normal   296.55 | error   -43.86
cuantil 0.50 -> observado   533.63 | normal   803.76 | error  -270.13
cuantil 0.75 -> observado  1068.68 | normal  1310.97 | error  -242.29
cuantil 0.95 -> observado  2431.76 | normal  2040.67 | error  +391.08
cuantil 0.99 -> observado  3058.59 | normal  2553.15 | error  +505.44

z = (v['monto'] - v['monto'].mean()) / v['monto'].std()

print('media de z:      %.6f' % z.mean())
print('desviacion de z: %.6f' % z.std())
print('asimetria de z:  %.4f' % z.skew())
print()
print('la venta mas grande esta a %.2f desviaciones' % z.max())
print('la mas chica, a %.2f' % z.min())

media de z:      -0.000000
desviacion de z: 1.000000
asimetria de z:  1.3405

la venta mas grande esta a 4.57 desviaciones
la mas chica, a -4.39

generador = np.random.default_rng(11)

for n in [100, 500, 2000, 5000]:
    casi = np.concatenate([
        generador.normal(size=int(n * 0.99)),
        generador.normal(loc=6, size=int(n * 0.01)),
    ])
    print('n=%6d con 1%% de intrusos -> p = %.6f'
          % (n, stats.shapiro(casi).pvalue))

n=   100 con 1% de intrusos -> p = 0.000005
n=   500 con 1% de intrusos -> p = 0.000000
n=  2000 con 1% de intrusos -> p = 0.000000
n=  5000 con 1% de intrusos -> p = 0.000000

def que_forma_tiene(serie, nombre):
    s = serie.dropna()
    a, k = s.skew(), s.kurtosis()

    if abs(a) > 1:
        forma = 'cola larga hacia la ' + ('derecha' if a > 0 else 'izquierda')
    elif k < -1:
        forma = 'plana, parece uniforme'
    elif abs(a) < 0.5 and abs(k) < 1:
        forma = 'razonablemente acampanada'
    else:
        forma = 'algo torcida, mirala'

    return ('%-13s asimetria %+7.4f | curtosis %+7.4f -> %s'
            % (nombre, a, k, forma))


for col in ['monto', 'unidades', 'satisfaccion', 'descuento']:
    print(que_forma_tiene(v[col], col))

monto         asimetria +1.3405 | curtosis +1.5458 -> cola larga hacia la derecha
unidades      asimetria +0.1843 | curtosis -0.3351 -> razonablemente acampanada
satisfaccion  asimetria -0.0149 | curtosis -1.3089 -> plana, parece uniforme
descuento     asimetria -0.0157 | curtosis -1.1604 -> plana, parece uniforme

Lo que te llevas

• 👀 El histograma decide más que cualquier prueba. Diez líneas de texto y ya sabías que monto no era una campana.
• 📐 Asimetría dice hacia dónde cae la cola; curtosis, si es picuda o plana. Una curtosis de -1,3 en cinco valores es la firma de una uniforme.
• 🎭 Las pruebas de normalidad no contestan lo que crees. Los mismos datos dan p = 0,44 con 50 filas y p = 0,000000 con 1.000.
• 🧭 La pregunta útil no es "¿es normal?" sino "¿cuánto me cuesta suponer que lo es?". En unidades el error de cuantiles era de una décima; en monto, la normal se inventaba ventas de -433 soles.
• 🪄 El logaritmo endereza colas por la derecha: la asimetría del monto pasó de 1,4278 a -0,1585. Pero cuesta 21 filas y complica la interpretación.
• 📏 Estandarizar no normaliza. La z deja media 0 y desviación 1, y la asimetría intacta.

Qué viene ahora

Han salido tres veces ya: las ventas negativas y las de cuatro mil soles. En el capítulo 6 vamos a por ellas, con las tres formas de detectarlas y la pregunta que casi nadie hace, que es si hay que quitarlas o no 🎯