Estadística y probabilidad básicas

El problema

Una cooperativa registra diez manzanas sintéticas de dos huertas. Cada manzana tiene dos medidas en milímetros, un peso en gramos y una marca de daño:

Llega una manzana nueva de 74 mm de diámetro y 65 mm de altura. Te piden encontrar la manzana registrada más parecida. También recibes una afirmación: «la huerta A da manzanas más pesadas y se daña menos».

Vectores y producto punto

Un vector es una lista ordenada de números. Cada posición conserva el mismo significado. Aquí cada manzana es un vector de dos coordenadas: diámetro y altura, ambas en milímetros.

import numpy as np

ids = np.array(
    ["M-01", "M-02", "M-03", "M-04", "M-05", "M-06", "M-07", "M-08", "M-09", "M-10"]
)
huerta = np.array(["A", "A", "A", "A", "A", "B", "B", "B", "B", "B"])
medidas = np.array(
    [
        [72, 65],
        [70, 63],
        [75, 68],
        [68, 62],
        [84, 74],
        [73, 66],
        [71, 64],
        [74, 67],
        [70, 63],
        [72, 65],
    ],
    dtype=np.float64,
)
peso = np.array(
    [150, 146, 158, 141, 205, 152, 149, 155, 147, 151], dtype=np.float64
)
danada = np.array(
    [False, False, True, False, False, True, False, False, True, False]
)

print(medidas[:3])
print(medidas.shape)

[[72. 65.]
 [70. 63.]
 [75. 68.]]
(10, 2)

La celda plegada copia la tabla en arreglos de NumPy. medidas tiene una fila por manzana y dos columnas, diámetro y altura. ids, huerta, peso y danada siguen el mismo orden, de M-01 a M-10.

El producto punto de dos vectores multiplica posición por posición y suma los productos: u @ v = u[0] * v[0] + u[1] * v[1]. Con el vector de coeficientes w = [0.5, 0.5], el producto punto combina las dos medidas de una manzana en un solo número. np.dot(u, w) es la forma de función del operador @.

Pregunta. Calcula a mano u @ w para M-01. Después predice la forma de medidas @ w.

u = medidas[0]
w = np.array([0.5, 0.5])

print(u @ w)
print(np.dot(u, w))
tamano = medidas @ w
print(tamano)
print(tamano.shape)

68.5
68.5
[68.5 66.5 71.5 65.  79.  69.5 67.5 70.5 66.5 68.5]
(10,)

Norma y distancia

La norma de un vector mide su longitud: ‖u‖ = √(u @ u). La distancia euclidiana entre dos vectores es la norma de su diferencia: ‖u − v‖. En dos coordenadas es el teorema de Pitágoras.

np.linalg.norm calcula la norma de un vector, o de cada fila si recibe axis=1. np.argmin devuelve la posición del valor más pequeño.

Pregunta. Calcula a mano la distancia entre M-01 (72, 65) y la manzana nueva (74, 65). Después elige, mirando la tabla, qué manzana registrada quedará más cerca.

nueva = np.array([74.0, 65.0])

print(round(float(np.sqrt(u @ u)), 3))
print(round(float(np.linalg.norm(u)), 3))

distancias = np.linalg.norm(medidas - nueva, axis=1)
print(np.round(distancias, 3))

cercana = np.argmin(distancias)
print(ids[cercana], round(float(distancias[cercana]), 3))

97.0
97.0
[ 2.     4.472  3.162  6.708 13.454  1.414  3.162  2.     4.472  2.   ]
M-06 1.414

Centro y dispersión

La media suma los valores y divide entre n. La mediana es el valor central después de ordenar. La varianza promedia los cuadrados de las desviaciones respecto de la media. La desviación estándar es su raíz cuadrada y vuelve a la unidad original: gramos.

peso[huerta == "A"] usa una máscara booleana para quedarse con los pesos de A. np.mean, np.median, np.var y np.std calculan los cuatro resúmenes.

Pregunta. Para la huerta A, ¿qué será mayor: la media o la mediana? ¿Qué huerta tendrá mayor desviación estándar?

peso_a = peso[huerta == "A"]
peso_b = peso[huerta == "B"]

for nombre, valores in [("A", peso_a), ("B", peso_b)]:
    print(
        nombre,
        np.mean(valores),
        np.median(valores),
        round(float(np.var(valores)), 2),
        round(float(np.std(valores)), 2),
    )

A 160.0 150.0 537.2 23.18
B 150.8 151.0 7.36 2.71

M-05 pesa 205 g, muy por encima del resto. Falta medir cuánto depende cada resumen de esa sola manzana.

Pregunta. Si quitas M-05 de la huerta A, ¿qué cambiará más: la media o la mediana?

Probabilidad condicional

La probabilidad empírica de un evento es su frecuencia en los datos: casos favorables entre casos totales. Con un arreglo booleano, sum() cuenta los True y mean() da la proporción.

La probabilidad condicional P(dañada | B) cuenta solo dentro de la huerta B: manzanas dañadas de B entre manzanas de B. P(B | dañada) cuenta dentro de las manzanas dañadas.

En el código, ~es_b invierte la máscara (huerta A) y danada & es_b exige las dos condiciones. danada.size es el total de elementos.

Pregunta. Escribe como fracción P(dañada), P(dañada | A), P(dañada | B) y P(B | dañada).

es_b = huerta == "B"

print(danada.sum(), danada.size)
print(danada.mean())
print(danada[~es_b].mean())
print((danada & es_b).sum(), es_b.sum())
print(danada[es_b].mean())
print(round(float(es_b[danada].mean()), 3))  # B entre las dañadas

3 10
0.3
0.2
2 5
0.4
0.667

Variación entre muestras

Diez manzanas son una muestra de lo que producen las huertas. Otra muestra daría otra media. Para ver cuánto cambia, simula una huerta con peso medio de 150 g y desviación estándar de 8 g. Después toma 2000 muestras de 5 manzanas y 2000 muestras de 50.

np.random.default_rng(0) crea un generador aleatorio con semilla 0, así que los números se repiten en cada ejecución. rng.normal(150, 8, size=(2000, 5)) devuelve una matriz de 2000 filas con 5 pesos cada una, y .mean(axis=1) deja la media de cada fila.

Pregunta. ¿Qué conjunto de medias variará más: el de muestras de 5 o el de muestras de 50?

rng = np.random.default_rng(0)

medias_5 = rng.normal(150, 8, size=(2000, 5)).mean(axis=1)
medias_50 = rng.normal(150, 8, size=(2000, 50)).mean(axis=1)

print(round(float(medias_5.std()), 2), round(float(medias_50.std()), 2))
print(round(float(medias_5.min()), 1), round(float(medias_5.max()), 1))
print(round(float(medias_50.min()), 1), round(float(medias_50.max()), 1))

3.57 1.14
135.8 162.5
146.3 153.5

Ahora simula dos huertas idénticas. Si una diferencia aparece a menudo entre huertas iguales, observarla no basta para afirmar que las huertas son distintas.

np.abs toma el valor absoluto de cada diferencia.

Pregunta. Con dos huertas idénticas y 5 manzanas por huerta, ¿con qué frecuencia aproximada las medias diferirán en 2 g o más? ¿Con qué frecuencia ocurrirá con 50 manzanas por huerta?

media_x = rng.normal(150, 8, size=(2000, 5)).mean(axis=1)
media_y = rng.normal(150, 8, size=(2000, 5)).mean(axis=1)
print((np.abs(media_x - media_y) >= 2).mean())

media_x_50 = rng.normal(150, 8, size=(2000, 50)).mean(axis=1)
media_y_50 = rng.normal(150, 8, size=(2000, 50)).mean(axis=1)
print((np.abs(media_x_50 - media_y_50) >= 2).mean())

0.6955
0.208

Para el daño, otra simulación marca cada manzana como dañada con probabilidad 0.3, de forma independiente, y cuenta las dañadas de cada huerta.

Pregunta. Si ambas huertas tienen probabilidad de daño 0.3, ¿con qué frecuencia diferirá en al menos una manzana el número de dañadas entre dos huertas de 5?

Práctica

Traza guiada

Para estos dos vectores, escribe antes de procesar: a @ b, np.linalg.norm(a) y la distancia entre a y b.

a = np.array([3.0, 4.0])
b = np.array([1.0, 2.0])

Comprueba con a @ b, np.linalg.norm(a) y np.linalg.norm(a - b). Si una predicción falla, escribe la suma de productos o la suma de cuadrados término por término.

Otro tamaño de muestra

Repite la simulación de medias con muestras de 20 manzanas. Antes de procesar, predice la desviación estándar de las medias con 8 / √20. Después calcula la frecuencia con la que dos huertas idénticas de 20 manzanas difieren en 2 g o más, y ubícala entre los resultados de 5 y 50.

Comparación propia

Construye dos arreglos de seis pesos inventados y un arreglo booleano de daño para cada uno. Calcula:

1. media, mediana y desviación estándar de cada grupo
2. P(dañada | grupo) para ambos grupos
3. la diferencia entre medias
4. la frecuencia de una diferencia igual o mayor en una simulación de dos grupos idénticos de seis valores; usa la media y la desviación estándar de los doce pesos juntos

Escribe una frase que diga si la diferencia observada es frecuente entre grupos iguales.

Reto opcional

Busca el menor tamaño de muestra n para el que dos huertas idénticas difieren en 2 g o más en menos del 5 % de 2000 simulaciones. Recorre n de 10 en 10 con la misma desviación estándar de 8 g.

Errores que revisar

• Si np.linalg.norm(medidas - nueva) devuelve un solo número, falta axis=1: sin eje, NumPy calcula la norma de toda la matriz.
• Si tu desviación estándar no coincide con la de pandas, revisa ddof: np.std usa ddof=0 y Series.std() usa ddof=1.
• Si la media y la mediana de un grupo están lejos, busca valores extremos antes de resumir con la media.
• Si confundes P(dañada | B) con P(B | dañada), escribe el denominador de cada una antes de calcular.
• Si la simulación da otros números al repetir una celda, vuelve a ejecutar desde la celda que crea rng = np.random.default_rng(0) y sigue el orden de la lección: cada llamada consume el generador.

Cierre

Responde a la afirmación «la huerta A da manzanas más pesadas y se daña menos» con dos frases: una sobre el peso y otra sobre el daño, cada una con un número de esta lección.

Lista de comprobación

• ¿Puedes calcular a mano el producto punto y la distancia entre dos vectores de dos coordenadas?
• ¿Puedes explicar por qué quitar M-05 cambia más la media que la mediana?
• ¿Distingues el denominador de P(dañada | B) y el de P(B | dañada)?
• ¿Puedes explicar por qué la media de 50 manzanas varía menos entre muestras que la de 5?

Documentación oficial

• numpy.dot: producto punto entre vectores y su relación con @.
• numpy.linalg.norm: normas de vectores y matrices, y el parámetro axis.
• Estadísticas en NumPy: mean, median, var, std y el parámetro ddof.
• Generador aleatorio de NumPy: default_rng, semillas y normal.