Hay veces, más de las que se cree, en que nos interesa estimar un cuantil en vez de la media. Si tenemos una variable dependinte \(y\) y una o varias independientes \(X\), lo que se suele hacer es una regresión cuantil.
Si visitamos uno de los papers originales de dicha técnica Computing Regression Quantiles vemos que trata de minimizar la siguiente expresión.
\[
\arg\min_b R_{\theta}(b) = \sum_{i = 1}^{n}\rho_\theta \left( y_i - x_ib\right)
\] Con \(\theta \in (0,1)\) y
\[
\begin{equation}
\rho_\theta(u) =
\begin{cases}
\theta u & u \geq 0\\
(\theta -1) & u < 0 \\
\end{cases}
\end{equation}
\]
Lo cual es simplemente “penalizar” por \(\theta\) cuando el residuo sea mayor o igual que 0, es decir, cuando nos equivocamos por arriba y por \((\theta -1)\) si nos equivocamos por abajo.
Ejemplo, si \(y_i = 40\) y \(f(x) = 50\) y queremos estimar el cuantil 0.95. Entonces como el residuo es menor que 0, se pondera por 0.05
\[\rho_(40 - 50) = (0.95 -1) (40 - 50) = 0.5 \] Si en cambio \(f(x) = 30\), es decir, nos equivocamos por abajo, pero a la misma distancia del valor real entonces
\[\rho(40-30) = 0.95 (40-30) = 9.5 \]
Y por tanto la función a minimizar \(\arg\min_b R_{\theta}(b)\) cuando \(\theta > 0.5\) va a tener un valor mucho mayor cuando nos “equivocamos” por abajo que por arriba. Y debido a cómo está definido \(\rho_\theta(u)\) se consigue la regresión cuantil con cuantil igual a \(\theta\). En el paper (de 1987) viene mejor explicado y el algoritmo para resolverlo en el caso de que \(f(x)\) sea lineal.
Warning
Pero, ¿qué me estás contando??
¡¡Estamos en la segunda década del siglo XXI y ahora todo es IA y Machín Lenin!!
Fuera coñas, el caso es que la gente de yandex en su librería catboost han utilizado esto para hacer la regresión cuantil, simplemente utilizando la expresión anterior como función de pérdida. Aquí se puede ver las diferentes funciones de pérdida que usan según el caso.
Como vemos, es lo mismo que se cuenta en el paper de 1987. Pero al meterlo como función de pérdida en el algoritmo sirve para el algoritmo de boosting que se utiliza en la librería.
Important
Ojo, que los de yandex han ido un poquito más allá
La gente de catboost, atinadamente ha dicho, y ¿por qué no construimos un función de pérdida que minimice globalmente varios cuantiles? Lo cual es algo así como “encuéntrame la distribución de los parámetros que mejor se ajusta a estos datos en vez de un sólo parámetro”.
Pero esto son arbolitos y boosting, no hay lo que se dice un parámetro de la función propiamente dicho, por lo que al final lo que se “aprende” debe ser la configuración de árboles que minimiza globalmente los cuantiles indicados.
Bueno, la función de pérdida “multi-quantile” es una modificación simple de la anterior.
El ejemplo no es mío, lo he visto por algún sitio que no me acuerdo.
Tip
catboost se puede utilizar en R y python.
Mostrar / ocultar código
import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as snsfrom catboost import CatBoostRegressorsns.set()n =800# X aleatoriasx_train = np.random.rand(n)x_test = np.random.rand(n)# un poquito de ruido gaussianonoise_train = np.random.normal(0, 0.3, n)noise_test = np.random.normal(0, 0.3, n)# Simulamos y_train e y _x como y = 2 + 3 * x + ruidoa, b =2, 3# al lioy_train = a * x_train + b + noise_trainy_test = a * x_test + b + noise_test
Pintamos
Mostrar / ocultar código
sns.scatterplot(x = x_train, y = y_train).set(title ="Ejemplillo")
Vaos a predecir 10 cuantiles
Mostrar / ocultar código
quantiles = [q/10for q inrange(1, 10)]# se ponen en string separados por commasquantile_str =str(quantiles).replace('[','').replace(']','')print(quantile_str)#> 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9
fig, ax = plt.subplots(figsize=(10, 6))ax.scatter(x_test, y_test)for col in ['pred_0.1', 'pred_0.5', 'pred_0.9']: ax.scatter(x_test.reshape(-1,1), preds[col], alpha=0.50, label=col)ax.legend()
Y ya estaría, no parece mala alternativa si uno tiene que hacer este tipo de cosas.
Tip
Ojalá le sirva a mi amigo Kenet para una cosa que estaba bicheando.
Pues poco más. Feliz domingo
Con R también se puede, como no.
Mostrar / ocultar código
library(reticulate) # para comunicar R y python y poder convertir datos y funciones de uno a otro bidireccionalmentelibrary(catboost)X_train <-as.matrix(py$x_train) # catboost en R espera una matrizY_train <-as.matrix(py$y_train)X_test <-as.matrix(py$x_test) Y_test <-as.matrix(py$y_test)head(X_train) ; head(Y_train)#> [,1]#> [1,] 0.6499924#> [2,] 0.9125034#> [3,] 0.4820761#> [4,] 0.2060426#> [5,] 0.5411845#> [6,] 0.1779657#> [,1]#> [1,] 4.165367#> [2,] 5.018479#> [3,] 4.146041#> [4,] 2.942043#> [5,] 4.287877#> [6,] 3.589826(quantiles_str <- py$quantile_str)#> [1] "0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9"
