# Regresión Lineal Desde una Perspectiva Probabilística



> Al día de hoy todos debemos haber ajustado una línea recta (plano, o hiperplano) $f(x) = x^T w$ a un conjunto de puntos $\{(x_i, y_i)\}_{i=1}^{N}$, con $x_i\in\mathbb{R}^{d}$ y $y_i\in\mathbb{R}$, y parámetros $w\in\mathbb{R}^{d}$. Esto normalmente se hace bajo el enfoque de mínimos cuadrados, es decir:
>
> $$
  \hat{w} = \arg \min_{w} ||y - Xw||^2,
  $$
  
> donde 
> 
> $$
  X = \left[\begin{array}{ccc}
  - & x_1^T  & - \\
  - & x_2^T  & - \\
    & \vdots &   \\
  - & x_N^T  & - \\
  \end{array}\right] \in \mathbb{R}^{N \times d}, \qquad y = \left[\begin{array}{ccc}
  y_1 \\
  y_2 \\
  \vdots \\
  y_N \\
  \end{array}\right] \in \mathbb{R}^{N}.
  $$

> Adicionalmente, para evitar el fenómeno de overfitting, una práctica común es incluir un término de regularización penalizando el tamaño de los parámetros $w$:
>
> $$
  \hat{w} = \arg \min_{w} ||y - Xw||^2 + \lambda ||w||^2
  $$

> Pues bien, hoy veremos cómo se explica todo esto desde una perspectiva probabilística.

> **Objetivos:**
> - Explicar la regresión lineal por mínimos cuadrados regularizada desde una perspectiva probabilística.

> **Referencias:**
> 
> - Bayesian Methods for Machine Learning course, HSE University, Coursera.


```python
# Importamos librerías
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
from matplotlib import pyplot as plt
import numpy as np
```


```python
# Generamos datos ficticios
N = 100
x = np.linspace(0, 1, N)
y = np.sin(2 * np.pi * x) + 0.3 * np.random.normal(size=(N,))
```


```python
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)
```


```python
# Gráfica de los datos
plt.figure(figsize=(6, 4))
plt.plot(x_train, y_train, 'or', ms=5, label="Datos de entrenamiento")
plt.plot(x_test, y_test, 'ob', ms=5, label="Datos de prueba")
plt.legend(loc="upper left", bbox_to_anchor=(1.05, 1))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.grid()
```


```python
# Definición del modelo
model_pipeline = Pipeline([
    ("poly", PolynomialFeatures(degree=3)),
    ("reg", LinearRegression())
])
```


```python
# Ajuste del modelo
model_pipeline.fit(x_train.reshape((-1, 1)), y_train)
```




    Pipeline(memory=None,
             steps=[('poly',
                     PolynomialFeatures(degree=3, include_bias=True,
                                        interaction_only=False, order='C')),
                    ('reg',
                     LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None,
                                      normalize=False))],
             verbose=False)




```python
# Visualización del modelo
plt.figure(figsize=(6, 4))
plt.plot(x_train, y_train, 'or', ms=5, label="Datos de entrenamiento")
plt.plot(x_train, model_pipeline.predict(x_train.reshape((-1, 1))), '*g', lw=2, label="Polinomio de grado 3 ajustado")
plt.legend(loc="upper left", bbox_to_anchor=(1.05, 1))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.grid()
```


```python
# Score sobre datos de entrenamiento
model_pipeline.score(x_train.reshape((-1, 1)), y_train)
```




    0.8299202521309362




```python
plt.figure(figsize=(6, 4))
plt.plot(x_test, y_test, 'or', ms=5, label="Datos de prueba")
plt.plot(x_test, model_pipeline.predict(x_test.reshape((-1, 1))), '*g', lw=2, label="Polinomio de grado 3 ajustado")
plt.legend(loc="upper left", bbox_to_anchor=(1.05, 1))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.grid()
```


```python
# Score sobre datos de prueba
model_pipeline.score(x_test.reshape((-1, 1)), y_test)
```




    0.8027169204325145




```python
# Definición del modelo
model_pipeline_ = Pipeline([
    ("poly", PolynomialFeatures(degree=8)),
    ("reg", LinearRegression())
])
```


```python
model_pipeline_.fit(x_train.reshape((-1, 1)), y_train)
```




    Pipeline(memory=None,
             steps=[('poly',
                     PolynomialFeatures(degree=8, include_bias=True,
                                        interaction_only=False, order='C')),
                    ('reg',
                     LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None,
                                      normalize=False))],
             verbose=False)




```python
# Visualización del modelo
plt.figure(figsize=(6, 4))
plt.plot(x_train, y_train, 'or', ms=5, label="Datos de entrenamiento")
plt.plot(x_train, model_pipeline_.predict(x_train.reshape((-1, 1))), '*g', lw=2, label="Polinomio de grado 8 ajustado")
plt.legend(loc="upper left", bbox_to_anchor=(1.05, 1))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.grid()
```


```python
# Score sobre datos de entrenaimiento
model_pipeline_.score(x_train.reshape((-1, 1)), y_train)
```




    0.8358048987005741




```python
plt.figure(figsize=(6, 4))
plt.plot(x_test, y_test, 'or', ms=5, label="Datos de prueba")
plt.plot(x_test, model_pipeline_.predict(x_test.reshape((-1, 1))), '*g', lw=2, label="Polinomio de grado 8 ajustado")
plt.legend(loc="upper left", bbox_to_anchor=(1.05, 1))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.grid()
```


```python
# Score sobre datos de prueba
model_pipeline_.score(x_test.reshape((-1, 1)), y_test)
```




    0.8137585467475051



### Overfitting


```python
model_pipeline.named_steps["reg"].coef_
```




    array([  0.        ,  13.71177813, -39.06580598,  25.87020616])




```python
model_pipeline_.named_steps["reg"].coef_
```




    array([    0.        ,     5.60686017,    68.77894649,  -523.40718161,
            1347.72005326, -1681.42816826,   958.15242705,  -100.12905775,
             -75.10849758])



### Regularización


```python
model_pipeline_reg = Pipeline([
    ("poly", PolynomialFeatures(degree=8)),
    ("reg", Ridge(alpha=0.01))
])
```


```python
model_pipeline_reg.fit(x_train.reshape((-1, 1)), y_train)
```




    Pipeline(memory=None,
             steps=[('poly',
                     PolynomialFeatures(degree=8, include_bias=True,
                                        interaction_only=False, order='C')),
                    ('reg',
                     Ridge(alpha=0.01, copy_X=True, fit_intercept=True,
                           max_iter=None, normalize=False, random_state=None,
                           solver='auto', tol=0.001))],
             verbose=False)




```python
# Visualización del modelo
plt.figure(figsize=(6, 4))
plt.plot(x_train, y_train, 'or', ms=5, label="Datos de entrenamiento")
plt.plot(x_train, model_pipeline_reg.predict(x_train.reshape((-1, 1))), '*g', lw=2, label="Polinomio de grado 8 ajustado")
plt.legend(loc="upper left", bbox_to_anchor=(1.05, 1))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.grid()
```


```python
model_pipeline_reg.score(x_train.reshape((-1, 1)), y_train)
```




    0.7863977898634911




```python
plt.figure(figsize=(6, 4))
plt.plot(x_test, y_test, 'or', ms=5, label="Datos de prueba")
plt.plot(x_test, model_pipeline_reg.predict(x_test.reshape((-1, 1))), '*g', lw=2, label="Polinomio de grado 8 ajustado")
plt.legend(loc="upper left", bbox_to_anchor=(1.05, 1))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.grid()
```


```python
model_pipeline_reg.score(x_test.reshape((-1, 1)), y_test)
```




    0.7731855600449242




```python
# Sin regularizacion
model_pipeline_.named_steps["reg"].coef_
```




    array([    0.        ,     5.60686017,    68.77894649,  -523.40718161,
            1347.72005326, -1681.42816826,   958.15242705,  -100.12905775,
             -75.10849758])




```python
# Regularizado
model_pipeline_reg.named_steps["reg"].coef_
```




    array([  0.        ,   5.72317488, -11.24124265,  -3.93324143,
             2.39265945,   4.53547432,   3.66532282,   1.12546011,
            -2.18374113])



## 1. Distribución normal

Antes de entrar en detalles repasemos un poco la distribución normal.

### 1.1. Distribución normal univariada

Se dice que una VA distribuye normal si su función de densidad de probabilidad es:

$$
\mathcal{N}(x|\mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left\{-\frac{(x - \mu)^2}{2 \sigma^2}\right\}
$$

con parámetros $\mu \in \mathbb{R}$: media de la VA X, y $\sigma^2 \in \mathbb{R}_{\geq0}$: varianza de la VA X.

**Ejercicio.** Demostrar que:

1. 
$$
\int_{-\infty}^{\infty} \mathcal{N}(x|\mu, \sigma^2) = 1.
$$

2. 
$$
\mu = \arg \max_{x} \mathcal{N}(x|\mu, \sigma^2)
$$

¿Cómo luce esta densidad?


```python
# Importar scipy.stats
from scipy import stats
# Importar numpy
import numpy as np
# Importar matplitlib.pyplot
from matplotlib import pyplot as plt
```


```python
# Tres VA normales
X = stats.norm(loc=0, scale=1)
Y = stats.norm(loc=1, scale=1)
Z = stats.norm(loc=0, scale=2)
```


```python
# Vector x para graficar
x = np.linspace(-6, 6, 100)
```


```python
# Graficar
fig, ax = plt.subplots(1, 1)
ax.plot(x, X.pdf(x), label=r"$\mu=0, \sigma=1$")
ax.plot(x, Y.pdf(x), label=r"$\mu=1, \sigma=1$")
ax.axvline(x=0, color='k')
ax.axvline(x=1, color='k')
ax.legend()
plt.show()
```


```python
fig, ax = plt.subplots(1, 1)
ax.plot(x, X.pdf(x), label=r"$\mu=0, \sigma=1$")
ax.plot(x, Z.pdf(x), label=r"$\mu=0, \sigma=2$")
ax.axvline(x=0, color='k')
ax.legend()
plt.show()
```

### 1.2. Distribución normal multivariada

Haciendo una extensión natural de lo anterior a vectores, obtenemos la función de densidad de probabilidad para distribución normal multivariada:

$$
\mathcal{N}(x|\mu, \Sigma) = \frac{1}{\sqrt{\det(2 \pi \Sigma)}} \exp \left\{-\frac{1}{2}(x - \mu)^T \Sigma^{-1} (x - \mu)\right\},
$$

con parámetros $\mu \in \mathbb{R}^n$: vector de medias de la VA X y $\Sigma \in \mathbb{R}^{n \times n}$: matriz de covarianzas de la VA X (simétrica y definida positiva).

¿Cómo luce esta densidad?


```python
# VA normales
X = stats.multivariate_normal(mean=np.zeros(2),
                              cov=np.eye(2))
Y = stats.multivariate_normal(mean=np.array([2, 3]),
                              cov=np.array([[1, 0.8],
                                            [0.8, 2]]))
```


```python
from mpl_toolkits import mplot3d
```


```python
# Gráfico
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
x, y = np.meshgrid(x, y)
z = X.pdf(np.dstack([x, y]))
fig = plt.figure(figsize=(8, 4))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(x, y, z)
ax2 = fig.add_subplot(122)
ax2.contour(x, y, z)
plt.show()
```


```python
# Gráfico
x = np.linspace(-2, 6, 100)
y = np.linspace(-3, 9, 100)
x, y = np.meshgrid(x, y)
z = Y.pdf(np.dstack([x, y]))
fig = plt.figure(figsize=(8, 4))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(x, y, z)
ax2 = fig.add_subplot(122)
ax2.contour(x, y, z)
plt.show()
```

## 2. Regresión lineal desde una perspectiva probabilística

Para ver algunos fundamentos de regresión lineal ver la introducción de este cuaderno.

Desde una perspectiva Bayesiana, podemos representar las VA objetivo $y$ y parámetros $w$, datos $X$ de un modelo lineal $y=Xw$ como:



¿Cuál es la factorización que induce este modelo?

$$
P(y, w, X) = P(w) P(X) P(y | w, X)
$$

Como no nos interesa modelar los datos, sino los parámetros, estudiamos la distribución posterior $P(w | X, y)$, la cual se puede reescribir como:

$$
P(w | X, y) = \frac{P(y, w | X)}{P(y | X)}.
$$

Por otra parte, usando la definición de probabilidad condicional y la factorización inducida, tenemos que:

$$
P(y, w | X) = P(y | X, w) P(w),
$$

de donde

$$
P(w | X, y) = \frac{P(y, w | X)}{P(y | X)} = \frac{P(y | X, w) P(w)}{P(y | X)}.
$$

Las distribuciones posteriores son en general difíciles de encontrar exactamente, dada la complejidad de modelar la distribución de evidencia $P(y | X)$. Durante el curso estaremos viendo técnicas para encontrarlas tanto de manera exacta (cuando sea posible analíticamente) como de manera aproximada.

Para este ejercicio, una suposición plausible es que la función de verosimilitud $P(y | X, w)$ y la distribución previa $P(w)$ son distribuciones normales de la forma:

$$
P(y | X, w) = \mathcal{N}(y | X w, \sigma^2 I),
$$

$$
P(w) = \mathcal{N}(w | 0, \gamma^2 I),
$$

mientras que la evidencia la podemos pensar como un factor de normalización (para que la distribución posterior integre a 1). En general nos interesa encontrar la **distribución posterior**, sin embargo, para este ejemplo, encontraremos sólo **la moda** de la misma.

De modo que para encontrar el valor más probable de los parámetros según la evidencia de los datos y el objetivo, hacemos:

$$
\arg \max_{w} P(w | X, y) = \arg \max_{w} \frac{P(y, w | X)}{P(y | X)} = \arg \max_{w} P(y, w | X) = \arg \max_{w} P(y | X, w) P(w)
$$

___
**Propiedad de los logaritmos:** dado que el logaritmo es una función creciente,

$$
\arg \max_x f(x) = \arg \max_x \log f(x) = \arg \min_x - \log f(x)
$$
___

Así que el anterior problema lo podemos reescribir como:

$$
\arg \min_w -\log P(y | X, w) - \log P(w)
$$

Desarrollando los logaritmos <font color=green>en el pizarrón</font>:

\begin{align}
-\log P(y | X, w) - \log P(w) & = \dots \\
                              & = \dots \\
                              & = \dots \\
                              & = \log C + \frac{1}{2\sigma^2} (y - X w)^T (y - X w) + \frac{1}{2\gamma^2} w^T w \\
                              & = \log C + \frac{1}{2\sigma^2} ||y - X w||^2 + \frac{1}{2\gamma^2} ||w||^2
\end{align}

Finalmente, podemos resolver el problema 

- Ridge:

$$
\arg \min_w ||y - X w||^2 + \lambda ||w||_2^2,
$$

- Lasso (previa de Laplace para los parámetros):

$$
\arg \min_w ||y - X w||^2 + \lambda ||w||_1,
$$

con $\lambda = \frac{\sigma^2}{\gamma^2}$. <font color=green> Interpretación de $\lambda$</font>.

Se puede demostrar que la solución es:

$$
w^{*}_{reg} = (X^T X + \lambda I)^{-1} X^T y
$$

**Ejemplo:** Usar mínimos cuadrados regularizados para ajustar una recta (plano) a los datos de precios de casas en Portland, OR.


```python
# Importar pandas
import pandas as pd
```


```python
# Leer datos
data = pd.read_csv("housing_prices.csv")
X = data.drop(columns=["price"])
y = data["price"]
X, y
```




    (    size  n_bedrooms
     0   2104           3
     1   1600           3
     2   2400           3
     3   1416           2
     4   3000           4
     5   1985           4
     6   1534           3
     7   1427           3
     8   1380           3
     9   1494           3
     10  1940           4
     11  2000           3
     12  1890           3
     13  4478           5
     14  1268           3
     15  2300           4
     16  1320           2
     17  1236           3
     18  2609           4
     19  3031           4
     20  1767           3
     21  1888           2
     22  1604           3
     23  1962           4
     24  3890           3
     25  1100           3
     26  1458           3
     27  2526           3
     28  2200           3
     29  2637           3
     30  1839           2
     31  1000           1
     32  2040           4
     33  3137           3
     34  1811           4
     35  1437           3
     36  1239           3
     37  2132           4
     38  4215           4
     39  2162           4
     40  1664           2
     41  2238           3
     42  2567           4
     43  1200           3
     44   852           2
     45  1852           4
     46  1203           3,
     0     399900
     1     329900
     2     369000
     3     232000
     4     539900
     5     299900
     6     314900
     7     198999
     8     212000
     9     242500
     10    239999
     11    347000
     12    329999
     13    699900
     14    259900
     15    449900
     16    299900
     17    199900
     18    499998
     19    599000
     20    252900
     21    255000
     22    242900
     23    259900
     24    573900
     25    249900
     26    464500
     27    469000
     28    475000
     29    299900
     30    349900
     31    169900
     32    314900
     33    579900
     34    285900
     35    249900
     36    229900
     37    345000
     38    549000
     39    287000
     40    368500
     41    329900
     42    314000
     43    299000
     44    179900
     45    299900
     46    239500
     Name: price, dtype: int64)




```python
# Importar sklearn.linear_model.Ridge
from sklearn.linear_model import Ridge, Lasso
# Importar sklearn.pipeline.Pipeline
from sklearn.pipeline import Pipeline
# Importart sklearn.preprocessing.StandardScaler
from sklearn.preprocessing import StandardScaler
```


```python
# Scaler-regressor pipeline
ridge_pipeline = Pipeline([
    # ("scaler", StandardScaler()),
    ("reg", Ridge(alpha=1))
])
lasso_pipeline = Pipeline([
    # ("scaler", StandardScaler()),
    ("reg", Lasso(alpha=1))
])
```


```python
# Aprendizaje
ridge_pipeline.fit(X, y)
```




    Pipeline(memory=None,
             steps=[('reg',
                     Ridge(alpha=1, copy_X=True, fit_intercept=True, max_iter=None,
                           normalize=False, random_state=None, solver='auto',
                           tol=0.001))],
             verbose=False)




```python
ridge_pipeline.named_steps["reg"].coef_, ridge_pipeline.named_steps["reg"].intercept_
```




    (array([  138.96771858, -8284.92788471]), 88647.59023379174)




```python
lasso_pipeline.fit(X, y)
```




    Pipeline(memory=None,
             steps=[('reg',
                     Lasso(alpha=1, copy_X=True, fit_intercept=True, max_iter=1000,
                           normalize=False, positive=False, precompute=False,
                           random_state=None, selection='cyclic', tol=0.0001,
                           warm_start=False))],
             verbose=False)




```python
lasso_pipeline.named_steps["reg"].coef_, lasso_pipeline.named_steps["reg"].intercept_
```




    (array([  139.20929338, -8735.4473684 ]), 89592.51877296084)




```python
# Gráfico
plt3d = plt.figure().gca(projection='3d')

xx1, xx2 = np.meshgrid(np.linspace(500, 4500), np.linspace(1, 5, 5))
yp_ridge = ridge_pipeline.predict(np.array([xx1.ravel(), xx2.ravel()]).T).reshape(xx1.shape)
yp_lasso = lasso_pipeline.predict(np.array([xx1.ravel(), xx2.ravel()]).T).reshape(xx1.shape)
plt3d.plot_surface(xx1, xx2, yp_ridge, alpha=0.2)
plt3d.plot_surface(xx1, xx2, yp_lasso, alpha=0.2)
plt3d.scatter(X['size'], X['n_bedrooms'], y);
plt3d.view_init(elev=10, azim=110)
```

### Distribución de Laplace


```python
X = stats.laplace(loc=0, scale=1)
Y = stats.laplace(loc=1, scale=1)
Z = stats.laplace(loc=0, scale=2)
```


```python
x = np.linspace(-5, 5, 100)
```


```python
# Graficar
fig, ax = plt.subplots(1, 1)
ax.plot(x, X.pdf(x), label=r"$\mu=0, b=1$")
ax.plot(x, Y.pdf(x), label=r"$\mu=1, b=1$")
ax.axvline(x=0, color='k')
ax.axvline(x=1, color='k')
ax.legend()
plt.show()
```


```python
# Graficar
fig, ax = plt.subplots(1, 1)
ax.plot(x, X.pdf(x), label=r"$\mu=0, b=1$")
ax.plot(x, Z.pdf(x), label=r"$\mu=0, b=2$")
ax.axvline(x=0, color='k')
ax.legend()
plt.show()
```

# Anuncios parroquiales:

## 1. Quiz de tema 3 para el miércoles 22 de Septiembre.



<footer id="attribution" style="float:right; color:#808080; background:#fff;">
Created with Jupyter by Esteban Jiménez Rodríguez.
</footer>

f08709a3a68845d20048a0b7386132a34ddf8f85
ipynb
Jupyter Notebook
modulo1/tema3/4_regresion_lineal.ipynb
G4ll4rd0/mebo2021
11096e6bb9c897c38ae02f9c30f90e1c878ee619
MIT
modulo1/tema3/4_regresion_lineal.ipynb
G4ll4rd0/mebo2021
11096e6bb9c897c38ae02f9c30f90e1c878ee619
MIT
modulo1/tema3/4_regresion_lineal.ipynb
G4ll4rd0/mebo2021
11096e6bb9c897c38ae02f9c30f90e1c878ee619
MIT
2022-02-03T22:50:53.000Z
2022-02-03T22:50:53.000Z
Qwen/Qwen-72B
1. YES
2. YES