La Dra. Grace Augustine fue una de las primeras personas en llegar a Pandora. Mientras elaboraba su libro sobre la flora de dicho planeta, conoció a los Na'vi, unos habitantes de Pandora que se caracterizaban por tener aspecto humanoide y medir hasta 3.9 mts. Ella creo una escuela para enseñarles sobre nuestra cultura y para construirla tuvo que entender mejor como se distribuyen las estaturas de los Na'vi. Felizmente ella, siendo xenobotánica, también dominaba la Ciencia de Datos, así que para lograrlo empezó con un censo de estaturas similar a este:
import random
_example_size = 500
_interval_size = 5
_xlim_min = 250
_xlim_max = 400
def get_random_serie(count, mu, sigma):
# Genera una serie de numeros aleatorios con distribución normal
for i in range(count):
yield round(random.gauss(mu,sigma),2)
naviHeights = list(get_random_serie(_example_size,300,24))
print(naviHeights)
[317.21, 293.32, 290.84, 305.5, 310.75, 245.35, 280.41, 280.03, 309.73, 338.67, 314.55, 289.66, 316.0, 288.67, 307.38, 263.22, 296.74, 281.21, 270.06, 362.5, 254.13, 296.14, 275.77, 303.41, 336.97, 268.0, 293.18, 265.16, 324.78, 327.91, 285.1, 309.59, 284.15, 290.76, 323.94, 313.13, 287.95, 305.12, 298.96, 306.19, 277.28, 310.06, 273.73, 275.59, 293.57, 353.44, 267.63, 301.13, 340.15, 272.62, 308.79, 268.84, 299.84, 306.83, 306.64, 344.65, 301.47, 307.3, 275.09, 325.78, 251.11, 255.86, 260.53, 294.76, 258.01, 259.5, 308.23, 312.96, 335.87, 306.72, 305.73, 330.2, 287.49, 330.86, 289.95, 310.6, 286.26, 325.15, 343.57, 248.48, 347.22, 277.41, 306.15, 333.74, 288.46, 283.55, 287.15, 291.3, 328.36, 283.58, 279.62, 297.41, 284.41, 318.76, 313.11, 267.13, 259.68, 284.37, 290.46, 291.32, 331.67, 285.27, 285.97, 286.68, 298.85, 316.03, 297.75, 311.56, 300.05, 309.17, 295.77, 341.28, 315.21, 325.3, 296.68, 351.12, 341.31, 299.75, 305.33, 304.11, 318.81, 300.8, 301.86, 279.54, 304.48, 316.38, 310.3, 353.04, 325.46, 316.42, 316.79, 287.61, 289.68, 266.27, 254.86, 275.24, 314.27, 331.26, 309.22, 303.24, 304.62, 306.54, 310.08, 299.73, 259.83, 294.84, 298.78, 322.32, 300.03, 332.8, 300.79, 314.62, 323.0, 330.61, 329.87, 302.1, 273.43, 313.06, 283.7, 300.35, 287.48, 318.76, 230.91, 311.25, 276.98, 360.4, 305.9, 323.84, 286.92, 330.34, 276.31, 331.73, 299.03, 306.26, 302.18, 311.52, 312.74, 266.8, 290.62, 306.87, 308.19, 247.96, 304.48, 301.51, 291.21, 317.72, 297.68, 299.67, 315.99, 337.94, 259.23, 292.71, 346.95, 285.9, 298.37, 265.36, 304.73, 295.7, 296.67, 330.58, 282.68, 334.5, 313.27, 309.56, 265.67, 295.64, 290.58, 321.26, 337.43, 316.36, 256.98, 280.89, 272.58, 281.69, 285.29, 346.57, 311.77, 279.2, 286.6, 271.57, 309.5, 300.41, 291.04, 319.87, 282.88, 295.76, 300.16, 318.55, 315.8, 303.29, 264.03, 313.73, 290.52, 251.59, 289.22, 314.77, 326.74, 308.62, 223.54, 321.06, 291.02, 267.94, 269.4, 276.72, 310.65, 295.72, 332.49, 280.11, 303.07, 268.56, 319.74, 291.34, 299.87, 299.59, 267.33, 335.72, 318.12, 318.23, 312.41, 289.66, 313.02, 302.81, 311.56, 289.35, 303.8, 282.63, 322.49, 256.23, 299.49, 330.87, 277.52, 280.82, 285.15, 312.41, 339.36, 306.71, 316.59, 281.93, 273.72, 342.22, 303.64, 340.75, 290.9, 285.09, 301.84, 312.74, 271.6, 266.79, 300.64, 292.23, 310.68, 329.21, 274.33, 296.75, 277.56, 307.22, 320.87, 286.27, 268.6, 318.61, 312.74, 307.82, 335.41, 294.24, 317.39, 314.25, 315.51, 276.3, 265.18, 249.31, 326.07, 307.21, 327.33, 285.01, 266.97, 342.32, 273.93, 299.42, 323.78, 290.46, 294.01, 268.41, 336.95, 311.46, 305.06, 327.49, 340.53, 313.34, 298.1, 295.71, 294.81, 346.41, 291.77, 288.35, 299.21, 303.93, 294.44, 296.33, 312.48, 330.48, 320.88, 327.56, 322.16, 327.64, 310.98, 262.61, 247.59, 286.96, 319.7, 284.71, 264.06, 305.15, 303.66, 297.2, 307.99, 339.89, 277.28, 282.74, 265.48, 299.64, 263.35, 323.58, 245.84, 344.34, 304.11, 288.09, 326.78, 286.0, 269.17, 275.41, 308.22, 304.47, 299.93, 301.26, 301.23, 287.82, 325.29, 261.72, 302.21, 294.09, 301.83, 305.92, 291.77, 319.44, 296.29, 298.15, 295.89, 310.33, 288.69, 295.04, 261.15, 285.92, 294.74, 270.31, 327.35, 278.52, 295.95, 283.06, 293.09, 325.34, 306.75, 291.77, 323.5, 308.77, 301.6, 320.16, 297.04, 291.6, 333.57, 301.52, 292.39, 258.96, 249.53, 313.76, 280.42, 289.33, 320.75, 278.92, 292.33, 302.64, 278.58, 307.32, 311.07, 311.26, 340.43, 322.85, 307.83, 334.54, 304.26, 299.77, 309.29, 270.42, 309.87, 299.91, 304.43, 292.88, 335.83, 298.22, 273.89, 311.03, 309.07, 310.6, 301.17, 279.0, 299.94, 295.22, 319.07, 290.73, 309.07, 304.07, 271.1, 292.56, 230.97, 246.65, 270.31, 279.28, 284.21, 302.96, 265.59, 293.31, 306.46, 297.94, 302.35, 338.31, 305.67, 294.9, 288.52, 284.06, 300.09, 328.33, 313.81, 328.78, 274.62, 326.15, 294.02, 287.35, 283.85, 266.49, 283.72, 322.68, 293.52, 283.65, 312.41, 282.32, 287.0, 280.15, 297.28, 276.53, 316.96, 243.65, 286.0, 314.78, 298.72, 266.07, 303.32, 300.35, 352.17, 278.63, 296.52, 327.05]
Si graficamos los datos en un histograma con intervalos de 5.0 cm, obtenemos lo siguiente:
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
colors = ["#0086CD", "#000"]
sns.set_context("talk")
sns.set_palette(sns.color_palette(colors))
sns.set_style("dark")
plt.subplots(figsize=(18,7))
plt.ylabel("Conteo")
plt.xlabel("Estatura (cm)")
sns.histplot(naviHeights, kde=False, binwidth=_interval_size, binrange=(_xlim_min,_xlim_max), label = "Estaturas")
plt.legend()
<matplotlib.legend.Legend at 0x15f5d668190>
Una de las tareas que se realizan en Machine Learning y Ciencia de Datos es ajustar un modelo a los datos para representar una distribución subyacente. Por ejemplo, podemos usar como modelo la distribución Normal (o Gaussiana) ya que esta tiene forma de campana y podría asemejarse a nuestro histograma. La ecuación de la distribución normal es $g(x) = \frac {1}{\sigma \sqrt{2\pi}} e^{ - \frac {(x - \mu)^2}{2\sigma^2} }$ y podemos apreciarla a continuación:
import numpy as np
import scipy.stats as stats
def get_random_model(count, mu, sigma):
serie = pd.Series(get_random_serie(count, mu, sigma))
sigma = serie.std()
mu = serie.mean()
x = np.linspace(mu - 3*sigma, mu + 3*sigma, serie.count())
return (x, stats.norm.pdf(x, mu, sigma))
normModel1 = get_random_model(_example_size,275,10)
plt.subplots(figsize=(18,7))
plt.xlabel("Estatura (cm)")
plt.ylabel("Densidad de probabilidad")
sns.lineplot(normModel1[0], normModel1[1], color="#000", label = "Modelo 1")
plt.show()
Una forma de cuantificar qué tan bueno es el ajuste de un modelo frente a otros modelos es calculando los "residuos", que no son más que la diferencia entre nuestros datos medidos y la predicción modelada para cada intervalo de histograma.
Hagámonos una primera idea con la siguiente imagen, todo lo que no es parte de la intersección de ambos gráficos serían como los residuos para nuestro modelo.
plt.subplots(figsize=(18,7))
plt.xlabel("Estatura (cm)")
plt.ylabel("Densidad de probabilidad")
sns.lineplot(normModel1[0], normModel1[1], color="#000", label = "Modelo 1")
plt.fill_between(normModel1[0], normModel1[1], alpha=0.3, color="#000")
sns.distplot(naviHeights, label = "Estaturas")
plt.legend()
<matplotlib.legend.Legend at 0x15f5d88eeb0>
Para cuantificar el rendimiento de un modelo con un solo número, podemos hacer lo siguiente: tomamos todos los residuos (la diferencia entre los datos medidos y los pronosticados), los elevamos al cuadrado y los sumamos. Si hacemos esto con todos los modelos, el modelo con menor $SSR$ sería el que tiene mejor rendimiento.
En nuestro problema, la ecuación de la distribución normal tiene 2 parámetros: $\mu$ (mu) y $\sigma$ (sigma). Esto quiere decir que, para cada par de valores de $\mu$ y $\sigma$, podríamos tener un $SSR$ diferente.
Para hallar el $SSR$ de nuestros modelos tendremos que construir dos vectores:
El primer vector tiene que ver con nuestro histograma. Debemos elegir el tamaño del intervalo de agrupación de estaturas y calcular la frecuencia relativa para cada uno de los intervalos, nuestro experimento podría perjudicarse si el intervalo es muy grande o muy pequeño para la cantidad de datos que tenemos, así que elijamos uno razonable como 5.0 cm. Ahora supongamos que $f(a,b)$ es una función que me devuelve la frecuencia relativa de un intervalo donde $a$ es el límite inferior y $b$ es el límite superior. Con esto ya podemos construir nuestro primer vector:
$$F = \begin{bmatrix} f(250.0,255.0) \\ f(255.0,260.0) \\ f(260.0,265.0) \\ . \\ . \\ . \end{bmatrix}$$El segundo vector lo construimos empleando la función de distribución normal $g(x) = \frac {1}{\sigma \sqrt{2\pi}} e^{ - \frac {(x - \mu)^2}{2\sigma^2} }$, los valores de $\mu$ y $\sigma$ dependerán del modelo en evaluación.
$$G = \begin{bmatrix} g(x_1) \\ g(x_2) \\ g(x_3) \\ . \\ . \\ . \end{bmatrix}$$Antes de empezar a calcular el $SSR$ debemos normalizar los vectores ya que las frecuencias relativas y densidades de probabilidad se encuentran en diferentes escalas.
$$|(norm(F) - norm(G))^2| = SSR$$import math
intervals = []
for i in range(round((_xlim_max-_xlim_min)/_interval_size)):
intervals.append([_xlim_min+i*_interval_size, _xlim_min+(i+1)*_interval_size])
def f_x(x):
# Calcula el valor de la frecuencia relativa correspondientes a una estatura x según los intervalos definidos
a = 0
b = 0
for item in intervals:
if((item[0] <= x) & (x <= item[1])):
a = item[0]
b = item[1]
array = np.array(naviHeights)
return ((a <= array) & (array < b)).sum() / array.size
def y_x(x, mu, sigma):
# Calcula el valor de y usando la función de la distribución normal
return 1/(sigma*math.sqrt(2*math.pi)) * (math.e**( -1 * (x-mu)**2 / (2 * sigma**2) ))
def get_ssr(mu, sigma):
y0 = np.array([f_x(x) for x in naviHeights])
y0 = y0 / np.linalg.norm(y0, keepdims=True)
y1 = np.array([y_x(x, mu, sigma) for x in naviHeights])
y1 = y1 / np.linalg.norm(y1, keepdims=True)
return np.sum(np.square(y0-y1))
rss1 = get_ssr(275,10)
normModel2 = get_random_model(_example_size,280,19)
rss2 = get_ssr(280,19)
normModel3 = get_random_model(_example_size,329,17)
rss3 = get_ssr(329,17)
normModel4 = get_random_model(_example_size,307,15)
rss4 = get_ssr(307,15)
plt.subplots(figsize=(18,7))
plt.xlabel("Estatura (cm)")
plt.ylabel("Densidad de probabilidad")
plt.xlim(_xlim_min,_xlim_max)
sns.lineplot(normModel1[0], normModel1[1], label = f'Modelo 1: μ = 275, σ = 10, SSR = {round(rss1,7)}', color="#000")
sns.lineplot(normModel2[0], normModel2[1], label = f'Modelo 2: μ = 280, σ = 19, SSR = {round(rss2,7)}', color="#f00")
sns.lineplot(normModel3[0], normModel3[1], label = f'Modelo 3: μ = 329, σ = 17, SSR = {round(rss3,7)}', color="#0f0")
sns.lineplot(normModel4[0], normModel4[1], label = f'Modelo 4: μ = 307, σ = 15, SSR = {round(rss4,7)}', color="#f0f")
plt.legend()
<matplotlib.legend.Legend at 0x15f5dda0790>
Con estos experimentos podríamos decir que el Modelo 4 es el mejor de los modelos propuestos, pero aún es muy temprano para asegurar que es el modelo definitivo o ideal. Veamos que pasa cuando calculamos el $SSR$ para cada valor de $\mu$ en el intervalo $[275,315]$ y $\sigma$ en el intervalo $[1,40]$:
mu_array = np.linspace(275, 315, 30)
sigma_array = np.linspace(1, 40, 30)
X, Y = np.meshgrid(mu_array, sigma_array)
Z = []
for i in range(30):
z = []
for j in range(30):
z.append(get_ssr(X[i][j], Y[i][j]))
Z.append(z)
plt.subplots(figsize=(18,7))
ax = plt.axes(projection='3d')
ax.contour3D(X, Y, Z, 50, cmap='hot', edgecolor='none')
ax.set_xlabel('μ')
ax.set_ylabel('σ')
ax.set_zlabel('SSR')
Text(0.5, 0, 'SSR')
Cada punto de esta superficie representa el $SSR$ para cada selección de parámetros $[\mu$, $\sigma$]. Notemos que el $SSR$ llega a su punto más bajo cuando $\mu$ se encuentra en el intervalo $[290,310]$ y $\sigma$ en el intervalo $[20,30]$ y, si lo vemos desde arriba, que cada anillo si bien tiene valores diferentes para $\mu$ y $\sigma$, tiene un valor constante para el $SSR$.
plt.subplots(figsize=(18,7))
ax = plt.axes(projection='3d')
ax.view_init(90, 90)
ax.contour3D(X, Y, Z, 50, cmap='hot', edgecolor='none')
ax.set_xlabel('μ')
ax.set_ylabel('σ')
ax.set_zlabel('SSR')
Text(0.5, 0, 'SSR')
Si empleamos los valores $[\mu$, $\sigma$] que estan a lo largo de una misma curva de nivel (las líneas rojas) no vamos a mejorar al $SSR$. Sin embargo, saltar a otro nivel o moverse perpendicularmente a ellos sí puede mejorar o reducir significativamente el rendimiento del modelo. Por lo que vemos aquí, podemos decir que $[297,22]$ tendría un mejor rendimiento que los anteriores modelos propuestos ya que se encontraría en un nivel más bajo.
normModel5 = get_random_model(_example_size,297,22)
rss5 = get_ssr(297,22)
plt.subplots(figsize=(18,7))
plt.xlabel("Estatura (cm)")
plt.ylabel("Densidad de probabilidad")
plt.xlim(_xlim_min,_xlim_max)
sns.lineplot(normModel4[0], normModel4[1], label = f'Modelo 4: μ = 307, σ = 14, SSR = {round(rss4,7)}', color="#f0f")
sns.lineplot(normModel5[0], normModel5[1], label = f'Modelo 5: μ = 297, σ = 22, SSR = {round(rss5,7)}', color="#00f")
plt.legend()
plt.subplots(figsize=(18,7))
plt.xlabel("Estatura (cm)")
plt.ylabel("Conteo")
sns.histplot(naviHeights, kde=False, binwidth=_interval_size, binrange=(_xlim_min,_xlim_max), label = "Estaturas")
plt.legend()
<matplotlib.legend.Legend at 0x15f5e135c70>