Cuándo la campana 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.

Y contéstate esto antes de seguir: ¿cuántas veces has usado la regla del 68% sin mirar si tus datos eran una campana? Yo, muchas, y por eso este capítulo existe 🔔

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 14 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ó.

Y dos más, de intentar meter en la campana lo que no cabe

stats.norm.fit(v['ciudad'])

TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Ajustar una normal es buscar la media y la desviación que mejor encajan, y las dos son cuentas. Sobre nombres de ciudad no hay nada que ajustar, y scipy lo dice con un mensaje bastante peor de lo que merece la pregunta 😅

stats.zscore(v['ciudad'])

TypeError: unsupported operand type(s) for /: 'str' and 'int'

La puntuación z es restar la media y dividir entre la desviación, y ahí se ve la división imposible en el mensaje. Todo este capítulo vive sobre variables numéricas, y las categóricas tienen su propio capítulo en el 13.

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

for col in ['monto', 'unidades', 'satisfaccion']:
    s = v[col].dropna()
    dentro = ((s > s.mean() - s.std()) & (s < s.mean() + s.std())).mean()
    print('%-13s dentro de una desviacion: %.1f%%   (la campana dice 68,3%%)'
          % (col, dentro * 100))

monto         dentro de una desviacion: 81.8%   (la campana dice 68,3%)
unidades      dentro de una desviacion: 68.7%   (la campana dice 68,3%)
satisfaccion  dentro de una desviacion: 59.4%   (la campana dice 68,3%)

import numpy as np
log = np.log(v['monto'][v['monto'] > 0])
print('asimetria del monto:      %.3f' % v['monto'].skew())
print('asimetria del logaritmo:  %.3f' % log.skew())
print('la prueba de normalidad sobre el logaritmo: p = %.6f'
      % stats.shapiro(log.sample(500, random_state=0)).pvalue)

asimetria del monto:      1.340
asimetria del logaritmo:  -0.158
la prueba de normalidad sobre el logaritmo: p = 0.000011

El 68% que salió 81,6%

Antes de cerrar, la trampa de este capítulo. La regla es correcta, la cuenta es correcta, y el resultado no cuadra 🔔

mu, sd = v['monto'].mean(), v['monto'].std()
dentro = ((v['monto'] > mu - sd) & (v['monto'] < mu + sd)).mean()

print(round(dentro * 100, 1))
# 81.8

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.

Y si de todo el capítulo te llevas una sola frase, que sea esta:

Y para mirar la forma de tus datos sin escribir código, un histograma en un tablero hace exactamente el mismo trabajo: guía de Power BI 📈

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 🎯

Preguntas frecuentes