In [2]:
"""
IPython Notebook v4.0 para python 2.7
Librerías adicionales: numpy, matplotlib
Contenido bajo licencia CC-BY 4.0. Código bajo licencia MIT. (c) Sebastian Flores.
"""
# Configuracion para recargar módulos y librerías
%reload_ext autoreload
%autoreload 2
%matplotlib inline
from IPython.core.display import HTML
HTML(open("style/mat281.css", "r").read())
Out[2]:
Consideremos los siguientes datos: link
In [5]:
%%bash
head data/x01.txt -n 60
In [4]:
import numpy as np
from matplotlib import pyplot as plt
# Plot of data
data = np.loadtxt("data/x01.txt", skiprows=33)
x = data[:,1]
y = data[:,2]
plt.figure(figsize=(16,8))
plt.plot(x, y, 'rs')
plt.xlabel("brain weight")
plt.ylabel("body weight")
plt.show()
In [6]:
import numpy as np
from matplotlib import pyplot as plt
# Plot of data
data = np.loadtxt("data/x01.txt", skiprows=33)
x = np.log(data[:,1])
y = np.log(data[:,2])
plt.figure(figsize=(16,8))
plt.plot(x, y, 'rs')
plt.xlabel("log brain weight")
plt.ylabel("log body weight")
plt.show()
In [7]:
%%bash
head data/x06.txt -n 40
In [9]:
%matplotlib gtk
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Plot of data
data = np.loadtxt("data/x06.txt", skiprows=37)
x = data[:,1]
y = data[:,2]
z = data[:,3]
fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, 'rs')
plt.xlabel("age [days]")
plt.ylabel("Water Temperature [C]")
plt.title("Length")
plt.show()
In [10]:
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Plot of data
data = np.loadtxt("data/x06.txt", skiprows=37)
x = np.log(data[:,1])
y = np.log(data[:,2])
z = np.log(data[:,3])
fig = plt.figure(figsize=(16,8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, 'rs')
plt.xlabel("age [days]")
plt.ylabel("Water Temperature [C]")
plt.title("Length")
plt.show()
Nuestra hipótesis de modelo lineal puede escribirse como
$$\begin{aligned} h_{\theta}(x) &= \theta_0 + \theta_1 x_1 + \theta_2 x_2 + ... + \theta_n x_n \\ &= \begin{bmatrix}\theta_0 & \theta_1 & \theta_2 & \dots & \theta_n\end{bmatrix} \begin{bmatrix}1 \\ x_1 \\x_2 \\ \vdots \\ x_n\end{bmatrix} \\ &= \theta^T \begin{bmatrix}1\\x\end{bmatrix} = \begin{bmatrix}1 & x^T\end{bmatrix} \theta \end{aligned}$$Definamos las matrices
$$\begin{aligned} Y &= \begin{bmatrix}y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(m)}\end{bmatrix}\end{aligned}$$y
$$\begin{aligned} X = \begin{bmatrix} 1 & x^{(1)}_1 & \dots & x^{(1)}_n \\ 1 & x^{(2)}_1 & \dots & x^{(2)}_n \\ \vdots & \vdots & & \vdots \\ 1 & x^{(m)}_1 & \dots & x^{(m)}_n \\ \end{bmatrix} = \begin{bmatrix} - (x^{(1)})^T - \\ - (x^{(2)})^T - \\ \vdots \\ - (x^{(m)})^T - \\ \end{bmatrix}\end{aligned}$$Luego la evaluación de todos los datos puede escribirse matricialmente como
$$\begin{aligned} X \theta &= \begin{bmatrix} 1 & x_1^{(1)} & ... & x_n^{(1)} \\ \vdots & \vdots & & \vdots \\ 1 & x_1^{(m)} & ... & x_n^{(m)} \\ \end{bmatrix} \begin{bmatrix}\theta_0 \\ \theta_1 \\ \vdots \\ \theta_n\end{bmatrix} \\ & = \begin{bmatrix} 1 \theta_0 + x^{(1)}_1 \theta_1 + ... + x^{(1)}_n \theta_n \\ \vdots \\ 1 \theta_0 + x^{(m)}_1 \theta_1 + ... + x^{(m)}_n \theta_n \\ \end{bmatrix} \\ & = \begin{bmatrix} h(x^{(1)}) \\ \vdots \\ h(x^{(m)}) \end{bmatrix}\end{aligned}$$Nuestro problema es encontrar un “buen” conjunto de valores $\theta$ de modo que
$$\begin{aligned} \begin{bmatrix} h(x^{(1)}) \\ h(x^{(2)}) \\ \vdots \\ h(x^{(m)}) \end{bmatrix} \approx \begin{bmatrix}y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(m)}\end{bmatrix}\end{aligned}$$es decir, que $$X \theta \approx Y$$
Para encontrar el mejor vector $\theta$ podríamos definir una función de costo $J(\theta)$ de la siguiente manera:
$$J(\theta) = \frac{1}{2} \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right)^2$$El mejor vector $\theta$ sería aquel que permite minimizar la norma 2 entre la predicción y el valor real.
Bueno, si $A \in \mathbb{R}^{m \times n}$, entonces $A^T \in \mathbb{R}^{n \times m}$ y la multiplicación está bien definida y obtengo un sistema lineal $n \times n$. $$(A^T A) \ \theta = A^T b$$ Si la matriz $A^T A$ es invertible, el sistema se puede solucionar “sin mayor reparo”. $$\theta = (A^T A)^{-1} A^T b$$
En nuestro caso, obtendríamos $$\theta = (X^T X)^{-1} X^T Y$$ Esta respuesta, aunque correcta, no admite interpretaciones y no permite generalizar a otros casos más generales.
En particular...
¿Qué relación tiene con la función de costo (no) utilizada?
¿Qué pasa si $A^T A$ no es invertible?
¿Cómo podemos obtener una buena aproximación para $\theta$?
Queremos encontrar $\theta^*$ que minimice $J(\theta)$.
Basta con utilizar una buena rutina de optimización para cumplir con dicho objetivo.
En particular, una elección natural es tomar la dirección de mayor descenso, es decir, el método del máximo descenso (gradient descent). $$\theta^{(n+1)} = \theta^{(n)} - \alpha \nabla_{\theta} J(\theta^{(n)})$$ donde $\alpha >0$ es la tasa de aprendizaje.
En nuestro caso, puesto que tenemos $$J(\theta) = \frac{1}{2} \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right)^2$$ se tiene que
$$\begin{aligned} \frac{\partial J(\theta)}{\partial \theta_k} &= \frac{\partial }{\partial \theta_k} \frac{1}{2} \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right)^2 \\ &= \frac{1}{2} \sum_{i=1}^{m} 2 \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) \frac{\partial h_{\theta}(x^{(i)})}{\partial \theta_k} \\ &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) x^{(i)}_k\end{aligned}$$Este algoritmo se llama Least Mean Squares
$$\begin{aligned} \theta^{(n+1)} & = \theta^{(n)} - \alpha \nabla_{\theta} J(\theta^{(n)}) \\ \frac{\partial J(\theta)}{\partial \theta_k} &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) x^{(i)}_k\end{aligned}$$OBS: La elección de $\alpha$ es crucial para la convergencia. En particular, $0.01/m$ funciona bien.
In [24]:
import numpy as np
from numpy.linalg import norm
def lms_regression_slow(X, Y, theta, tol=1E-6):
converged = False
alpha = 0.01/len(Y)
while not converged:
gradient = 0.
for xiT, yi in zip(X,Y):
hi = np.dot(theta, xiT)
gradient += (hi - yi)*xiT.T
new_theta = theta - alpha * gradient
converged = norm(theta-new_theta) < tol * norm(theta)
theta = new_theta
return theta
def lms_regression_fast(X, Y, theta, tol=1E-6):
converged = False
alpha = 0.01/len(Y)
theta = theta.reshape(X.shape[1], 1)
A = np.dot(X.T,X)
b = np.dot(X.T, Y)
while not converged:
gradient = np.dot(A, theta) - b
new_theta = theta - alpha * gradient
converged = norm(theta-new_theta) < tol * norm(theta)
theta = new_theta
return theta
m = 1000
t = np.linspace(0,1,m)
x = 2 + 2*t
y = 300 + 100*t
X = np.array([np.ones(m), x]).T
Y = y.reshape(m,1)
theta_0 = np.array([[0.0,0.0]])
theta = lms_regression_fast(X, Y, theta_0)
print theta
Se tiene
$$\begin{aligned} J(\theta) &= \frac{1}{2} \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right)^2 \\ &= \frac{1}{2} \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) \\ &= \frac{1}{2} \left( X \theta - Y \right)^T \left( X \theta - Y \right) \\ &= \frac{1}{2} \left( \theta^T X^T - Y^T \right) \left( X \theta - Y \right) \\ &= \frac{1}{2} \left( \theta^T X^T X \theta - \theta^T X^T Y - Y^T X \theta + Y^T Y \right) \\ &= \frac{1}{2} \left( \theta^T X^T X \theta - 2 (Y^T X) \theta + Y^T Y \right)\end{aligned}$$Aplicando a cada uno de los términos, obtenemos:
$$\begin{aligned} \nabla_\theta ( \theta^T X^T X \theta ) &= X^T X \theta + (X^T X)^T \theta \\ & = 2 X^T X \theta\end{aligned}$$también se tiene
$$\begin{aligned} \nabla_\theta ( Y^T X \theta ) &= (Y^T X) ^T\\ &= X^T Y\end{aligned}$$y por último
$$\begin{aligned} \nabla_\theta ( Y^T Y ) = 0\end{aligned}$$Esto significa que el problema $$\min_\theta J(\theta)$$ se resuelve al hacer todas las derivadas parciales iguales a cero (ie, gradiente igual a cero) $$\nabla_\theta J(\theta) = 0$$ lo cual en nuestro caso se convierte convenientemente a la ecuación normal $$X^T X \theta = X^T Y$$ y se tiene $$\theta = (X^T X)^{-1} X^T Y$$
In [21]:
import numpy as np
from numpy.linalg import norm
def matrix_regression(X, Y, theta, tol=1E-6):
A = np.dot(X.T,X)
b = np.dot(X.T,Y)
sol = np.linalg.solve(A,b)
return sol.flatten()
m = 100
t = np.linspace(0,1,m)
x = 2 + 2*t
y = 300 + 100*t
X = np.array([np.ones(m), x]).T
Y = y.reshape(m,1)
theta_0 = np.array([[0.0,0.0]])
theta = matrix_regression(X, Y, theta_0)
print theta
Asumamos que outputs e inputs están relacionados mediante $$y^{(i)}= \theta^T x^{(i)}+ \varepsilon^{(i)}$$ donde $\varepsilon^{(i)}$ es un error que captura efectos sin modelar o ruido de medición.
Supongamos que los $\varepsilon^{(i)}$ se distribuyen de manera idéntica e independientemente de acuerdo a una distribución gausiana de media $0$ y varianza $\sigma^2$. $$\varepsilon^{(i)}\sim \mathcal{N}(0, \sigma^2)$$
Cabe destacar que:
$\theta$ no es una variable aleatoria, es un parámetro (desconocido).
$\varepsilon^{(i)}$, $x^{(i)}$ y $y^{(i)}$ son variables aleatorias.
$\varepsilon^{(i)}\sim \mathcal{N}(0, \sigma^2)$
$y^{(i)} \ | \ x^{(i)}; \theta \sim \mathcal{N}(\theta^T x^{(i)}, \sigma^2)$ pues $y^{(i)}= \theta^T x^{(i)}+ \varepsilon^{(i)}$
Tenemos entonces que $$\mathbb{P}[\varepsilon^{(i)}] = \frac{1}{\sqrt{2\pi}\sigma} \exp\left( -\frac{(\varepsilon^{(i)})^2}{2\sigma^2} \right)$$ y por tanto $$\mathbb{P}[y^{(i)}\ | \ x^{(i)}; \theta ] = \frac{1}{\sqrt{2\pi}\sigma} \exp\left( -\frac{(y^{(i)}- \theta^T x^{(i)})^2}{2\sigma^2} \right)$$
La función de verosimilitud $L(\theta)$ nos permite entender que tan probable es encontrar los datos observados, para una elección del parámetro $\theta$.
$$\begin{aligned} L(\theta) &= \prod_{i=1}^{m} \mathbb{P}[y^{(i)}| x^{(i)}; \theta ] \\ &= \prod_{i=1}^{m} \frac{1}{\sqrt{2\pi}\sigma} \exp\left( -\frac{(y^{(i)}- \theta^T x^{(i)})^2}{2\sigma^2} \right) \\ &= \frac{1}{(\sqrt{2\pi}\sigma)^{m}} \exp\left( - \sum_{i=1}^{m} \frac{(y^{(i)}- \theta^T x^{(i)})^2}{2\sigma^2} \right)\end{aligned}$$Nos gustaría encontrar el parámetro $\theta$ que más probablemente haya generado los datos observados, es decir, el parámetro $\theta$ que maximiza la función de verosimilitud.
Es decir, la función costo $J(\theta)$ cuadrática puede interpretarse como un intento por encontrar el parámetro $\theta$ más probable bajo la hipótesis que el error en $y$ es gaussiano.
Al realizar regresión, conviene normalizar/estandarizar los datos, es decir transformarlos para que tengan una escala común:
Utilizando la media y la desviación estándar $$\frac{x_i-\overline{x_i}}{\sigma_{x_i}}$$
Utilizando mínimos y máximos $$\frac{x_i-\min{x_i}}{\max{x_i} - \min{x_i}}$$
In [27]:
from sklearn import datasets
# The data
boston = datasets.load_boston()
X = boston.data
Y = boston.target
#print X
#print Y
# What’s the data
print boston.DESCR
In [29]:
from sklearn import datasets
import numpy as np
#from mat281_code.lms import lms_regression as lms
from sklearn import linear_model
# The data
boston = datasets.load_boston()
X = boston.data
Y = boston.target
m = X.shape[0]
# Normalization of data
#X_train_aux = (X-X.min(axis=0))/(X.max(axis=0)-X.min(axis=0))
#Y_train_aux = (Y-Y.min())/(Y.max()-Y.min())
X_train_aux = (X-X.mean(axis=0))/X.std(axis=0)
Y_train_aux = (Y-Y.mean())/Y.std()
# Put in shape for normal equations
X_train = np.hstack([np.ones([m,1]), X_train_aux])
Y_train = Y_train_aux.reshape(m,1)
# Direct Solution
theta = np.linalg.solve(np.dot(X_train.T, X_train),
np.dot(X_train.T, Y_train))
print theta.flatten()
# sklearn solution
regr = linear_model.LinearRegression()
regr.fit(X_train_aux, Y_train_aux)
theta = regr.intercept_, regr.coef_
print theta
# LMS Solution - Tarda mucho
#theta0 = Y_train.mean()/X_train.mean(axis=0)/X.shape[1]
#theta = lms(X_train, Y_train, theta0)
#print theta
Las ecuaciones normales y sklearn entregan el mismo resultado:
$$\begin{aligned} \theta = (& 0.00, -0.10, 0.12, 0.02, 0.07, -0.22, \\ & 0.29, 0.00, -0.34, 0.29, -0.23, -0.22, 0.09, -0.41 )\end{aligned}$$Mientras que el algoritmo lms entrega
$$\begin{aligned} \theta = (&0.00, -0.10, 0.12, 0.02, 0.07, -0.21, \\ &0.29, 0.00, -0.34, 0.29, -0.23, -0.21, 0.09, -0.41 )\end{aligned}$$Si las variables son
CRIM, ZN, INDUS, CHAS, NOX, RM, AGE, DIS, RAD, TAX, PTRATIO, B, LSTAT, MEDV
¿Cuáles variables tienen más impacto en el precio de la vivienda?
$\theta_0=+0.00$.
$\theta_1 = -0.10$: CRIM, per capita crime rate by town.
$\theta_2 = +0.12$: ZN, proportion of residential land zoned for lots over 25,000 sq.ft.
$\theta_3 = +0.02$: INDUS, proportion of non-retail business acres per town
$\theta_4 = +0.07$: CHAS, Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
$\theta_5 = -0.22$: NOX, nitric oxides concentration (parts per 10 million)
$\theta_6 = +0.29$: RM, average number of rooms per dwelling
$\theta_7 = +0.00$: AGE, proportion of owner-occupied units built prior to 1940
$\theta_8 = -0.34$: DIS, weighted distances to five Boston employment centres
$\theta_9 = +0.29$: RAD, index of accessibility to radial highways
$\theta_{10} = -0.23$: TAX, full-value property-tax rate per \$10,000
$\theta_{11} = -0.22$: PTRATIO pupil-teacher ratio by town
$\theta_{12} = +0.09$: B, $1000(Bk - 0.63)^2$ where Bk is the proportion of blacks by town
$\theta_{13} = -0.41$: LSTAT, % lower status of the population
In [30]:
%matplotlib inline
from sklearn import datasets
import matplotlib.pyplot as plt
iris = datasets.load_iris()
def plot(dataset, ax, i, j):
ax.scatter(dataset.data[:,i], dataset.data[:,j], c=dataset.target, s=50)
ax.set_xlabel(dataset.feature_names[i], fontsize=20)
ax.set_ylabel(dataset.feature_names[j], fontsize=20)
# row and column sharing
f, ((ax1, ax2), (ax3, ax4), (ax5,ax6)) = plt.subplots(3, 2, figsize=(16,16))
plot(iris, ax1, 0, 1)
plot(iris, ax2, 0, 2)
plot(iris, ax3, 1, 2)
plot(iris, ax4, 0, 3)
plot(iris, ax5, 1, 3)
plot(iris, ax6, 2, 3)
f.tight_layout()
plt.show()
In [31]:
# REVISAR
import numpy as np
from sklearn import datasets
from sklearn import linear_model
# Loading the data
iris = datasets.load_iris()
X = iris.data
iris_label = iris.target
# Apply linear regression to each model
predictions = {}
regr = linear_model.LinearRegression(fit_intercept=True, normalize=False)
for label in range(0,3):
X_train = X[iris_label==label][:,:-1]
Y_train = X[iris_label==label][:,-1]
regr.fit(X_train, Y_train) # Still must add the column of 1
theta = regr.intercept_, regr.coef_
print theta
Y_pred = regr.predict(X_train)
predictions[label] = Y_pred
print "Error", np.linalg.norm(Y_train-Y_pred,2)/len(Y_pred)
In [32]:
%matplotlib inline
from sklearn import datasets
import matplotlib.pyplot as plt
iris = datasets.load_iris()
def plot(dataset, ax, i, j):
colors = {0:"r", 1:"b", 2:"g"}
markers = {0:"s", 1:"o", 2:"<"}
for label in range(3):
x = dataset.data[:,i][dataset.target==label]
y = dataset.data[:,j][dataset.target==label]
ax.scatter(x, y, c=colors[label], marker=markers[label], s=50)
if j==3:
ax.scatter(x, predictions[label], c="w", marker=markers[label],
s=50)
ax.set_xlabel(dataset.feature_names[i], fontsize=20)
ax.set_ylabel(dataset.feature_names[j], fontsize=20)
# row and column sharing
f, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(16,16))
plot(iris, ax1, 0, 3)
plot(iris, ax2, 1, 3)
plot(iris, ax3, 2, 3)
f.tight_layout()
plt.show()