In [7]:
"""
IPython Notebook v4.0 para python 2.7
Librerías adicionales: numpy, matplotlib, sklearn
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[7]:
Los datos corresponden a 3 cultivos diferentes de vinos de la misma región de Italia, y que han sido identificados con las etiquetas 1, 2 y 3. Para cada tipo de vino se realizado 13 análisis químicos:
La base de datos contiene 178 muestras distintas en total.
Si no conocemos de antemano las etiquetas, es decir, los cultivos 1,2 o 3 a los que pertenece cada muestra, el problema es de clustering:
$$\textrm{Tenemos } X \in R^{n,m} \textrm{ buscamos las etiquetas } Y \in N^m$$
¿Cuántos grupos existen? ¿A que grupo pertenece cada dato?
Si conocemos los valores y las etiquetas, y se desea obtener la etiqueta de una muestra sin etiquetar, el problema es de clasificación:
$$\textrm{Tenemos } X \in R^{n \times m} \textrm{ y } Y \in N^m \textrm{ y buscamos las etiquetas de } x \in R^n$$
¿A qué grupo pertenece este nuevo dato?
In [ ]:
%%bash
cat data/Challenger.txt
In [8]:
from matplotlib import pyplot as plt
import numpy as np
# Plot of data
data = np.loadtxt("data/Challenger.txt", skiprows=1)
x = data[:,0]
y = data[:,1]
plt.figure(figsize=(16,8))
plt.plot(x, y, 'bo', ms=8)
plt.title("Exito o Falla en lanzamiento de Challenger")
plt.xlabel(r"T [${}^o F$]")
plt.ylabel(r"Bad Rings")
plt.ylim([-0.1,3.1])
plt.show()
In [9]:
# Plot of data
data = np.loadtxt("data/Challenger.txt", skiprows=1)
x = (data[:,0]-32.)/1.8
y = np.array(data[:,1]==0,int)
plt.figure(figsize=(16,8))
plt.plot(x[y==0], y[y==0], 'bo', label="Falla", ms=8)
plt.plot(x[y>0], y[y>0], 'rs', label="Exito", ms=8)
plt.ylim([-0.1, 1.1])
plt.legend(loc=0, numpoints=1)
plt.title("Exito o Falla en lanzamiento de Challenger")
plt.xlabel(r"T [${}^o C$]")
plt.ylabel(r"$y$")
plt.show()
Definimos como antes
$$\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}\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}\end{aligned}$$
In [10]:
from matplotlib import pyplot as plt
import numpy as np
def sigmoid(z):
return (1+np.exp(-z))**(-1.)
z = np.linspace(-5,5,100)
g = sigmoid(z)
fig = plt.figure(figsize=(16,8))
plt.plot(z,sigmoid(z), lw=2.0)
plt.plot(z,sigmoid(z*2), lw=2.0)
plt.plot(z,sigmoid(z-2), lw=2.0)
plt.grid("on")
plt.show()
$g(z) = (1+e^{-z})^{-1}$ y $g'(z) = g(z)(1-g(z))$.
Demostración:
$$\begin{aligned} g'(z) &= \frac{-1}{(1+e^{-z})^2} (-e^{-z}) \\ &= \frac{e^{-z}}{(1+e^{-z})^2} \\ &= \frac{1}{1+e^{-z}} \frac{e^{-z}}{1+e^{-z}} \\ &= \frac{1}{1+e^{-z}} \left(1 - \frac{1}{1+e^{-z}} \right) \\ &= g(z)(1-g(z))\end{aligned}$$
In [11]:
from matplotlib import pyplot as plt
import numpy as np
def sigmoid(z):
return (1+np.exp(-z))**(-1.)
z = np.linspace(-5,5,100)
g = sigmoid(z)
dgdz = g*(1-g)
fig = plt.figure(figsize=(16,8))
plt.plot(z, g, "k", label="g(z)", lw=2)
plt.plot(z, dgdz, "r", label="dg(z)/dz", lw=2)
plt.legend()
plt.grid("on")
plt.show()
El cálculo del gradiente es directo:
$$\begin{aligned} \frac{\partial J(\theta)}{\partial \theta_k} &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) \frac{\partial}{\partial \theta_k} h_{\theta}(x^{(i)}) \\ &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) \frac{\partial}{\partial \theta_k} g(\theta^T x^{(i)}) \\ &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) h_{\theta}(x^{(i)}) \left(1-h_{\theta}(x^{(i)})\right) \frac{\partial}{\partial \theta_k} (\theta^T x^{(i)}) \\ &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) h_{\theta}(x^{(i)}) \left(1-h_{\theta}(x^{(i)})\right) x^{(i)}_k\end{aligned}$$¿Hay alguna forma de escribir todo esto de manera matricial? Recordemos que si las componentes eran
$$\begin{aligned} \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) x^{(i)}_k = \sum_{i=1}^{m} x^{(i)}_k \left( h_{\theta}(x^{(i)}) - y^{(i)}\right)\end{aligned}$$podíamos escribirlo vectorialmente como $$X^T (X\theta - Y)$$
Luego, para
$$\begin{aligned} \frac{\partial J(\theta)}{\partial \theta_k} &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) h_{\theta}(x^{(i)}) \left(1-h_{\theta}(x^{(i)})\right) x^{(i)}_k \\ &= \sum_{i=1}^{m} x^{(i)}_k \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) h_{\theta}(x^{(i)}) \left(1-h_{\theta}(x^{(i)})\right)\end{aligned}$$podemos escribirlo vectorialmente como $$\nabla_{\theta} J(\theta) = X^T \Big[ (g(X\theta) - Y) \odot g(X\theta) \odot (1-g(X\theta)) \Big]$$ donde $\odot$ es la multiplicación elemento a elemento (element-wise).
In [12]:
import numpy as np
def sigmoid(z):
return 1./(1+np.exp(-z))
def norm2_error_logistic_regression(X, Y, theta0, tol=1E-6):
converged = False
alpha = 0.01/len(Y)
theta = theta0
while not converged:
H = sigmoid(np.dot(X, theta))
gradient = np.dot(X.T, (H-Y)*H*(1-H))
new_theta = theta - alpha * gradient
converged = np.linalg.norm(theta-new_theta) < tol * np.linalg.norm(theta)
theta = new_theta
return theta
¿Es la derivación anterior probabilísticamente correcta?
Asumamos que la pertenencia a los grupos está dado por
$$\begin{aligned} \mathbb{P}[y = 1| \ x ; \theta ] & = h_\theta(x) \\ \mathbb{P}[y = 0| \ x ; \theta ] & = 1 - h_\theta(x)\end{aligned}$$Esto es, una distribución de Bernoulli con $p=h_\theta(x)$.\ Las expresiones anteriores pueden escribirse de manera más compacta como
$$\begin{aligned} \mathbb{P}[y | \ x ; \theta ] & = (h_\theta(x))^y (1 - h_\theta(x))^{(1-y)} \\\end{aligned}$$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} \Big(h_{\theta}(x^{(i)})\Big)^{y^{(i)}} \Big(1 - h_\theta(x^{(i)})\Big)^{(1-y^{(i)})}\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.
Calculamos la log-verosimilitud:
$$\begin{aligned} l(\theta) &= \log L(\theta) \\ &= \log \prod_{i=1}^{m} (h_\theta(x^{(i)}))^{y^{(i)}} (1 - h_\theta(x^{(i)}))^{(1-y^{(i)})} \\ &= \sum_{i=1}^{m} y^{(i)}\log (h_\theta(x^{(i)})) + (1-y^{(i)}) \log (1 - h_\theta(x^{(i)}))\end{aligned}$$No existe una fórmula cerrada que nos permita obtener el máximo de la log-verosimitud. Pero podemos utilizar nuevamente el método del gradiente máximo.
Recordemos que si
$$\begin{aligned} g(z) = \frac{1}{1+e^{-z}}\end{aligned}$$Entonces
$$\begin{aligned} g'(z) &= g(z)(1-g(z))\end{aligned}$$y luego tenemos que
$$\begin{aligned} \frac{\partial}{\partial \theta_k} h_\theta(x) &= h_\theta(x) (1-h_\theta(x)) x_k\end{aligned}$$Es decir, para maximizar la log-verosimilitud obtenemos igual que para la regresión lineal:
$$\begin{aligned} \theta^{(n+1)} & = \theta^{(n)} - \alpha \nabla_{\theta} l(\theta^{(n)}) \\ \frac{\partial l(\theta)}{\partial \theta_k} &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) x^{(i)}_k\end{aligned}$$Aunque, en el caso de regresión logística, se tiene $h_\theta(x)=1/(1+e^{-x^T\theta})$
OBS: La elección de $\alpha$ es crucial para la convergencia. En particular, $0.01/m$ funciona bien.
Es decir, para maximizar la log-verosimilitud obtenemos igual que para la regresión lineal:
$$\begin{aligned} \theta^{(n+1)} & = \theta^{(n)} - \alpha \nabla_{\theta} l(\theta^{(n)}) \\ \frac{\partial l(\theta)}{\partial \theta_k} &= \sum_{i=1}^{m} \left( h_{\theta}(x^{(i)}) - y^{(i)}\right) x^{(i)}_k\end{aligned}$$Aunque, en el caso de regresión logística, se tiene $h_\theta(x)=1/(1+e^{-x^T\theta})$
OBS: La elección de $\alpha$ es crucial para la convergencia. En particular, $0.01/m$ funciona bien.
In [13]:
import numpy as np
def likelihood_logistic_regression(X, Y, theta0, tol=1E-6):
converged = False
alpha = 0.01/len(Y)
theta = theta0
while not converged:
H = sigmoid(np.dot(X, theta))
gradient = np.dot(X.T, H-Y)
new_theta = theta - alpha * gradient
converged = np.linalg.norm(theta-new_theta) < tol * np.linalg.norm(theta)
theta = new_theta
return theta
def sigmoid(z):
return 1./(1+np.exp(-z))
In [14]:
# Plot of data
data = np.loadtxt("data/Challenger.txt", skiprows=1)
x = (data[:,0]-32.)/1.8
X = np.array([np.ones(x.shape[0]), x]).T
y = np.array(data[:,1]==0,int)
theta_0 = y.mean() / X.mean(axis=0)
print "theta_0", theta_0
theta_J = norm2_error_logistic_regression(X, y, theta_0)
print "theta_J", theta_J
theta_l = likelihood_logistic_regression(X, y, theta_0)
print "theta_l",theta_l
In [15]:
# Predictions
y_pred_J = sigmoid(np.dot(X, theta_J))
y_pred_l = sigmoid(np.dot(X, theta_l))
# Plot of data
plt.figure(figsize=(16,8))
plt.plot(x[y==0], y[y==0], 'bo', label="Falla", ms=8)
plt.plot(x[y>0], y[y>0], 'rs', label="Exito", ms=8)
plt.plot(x, y_pred_J, label="Norm 2 error prediction")
plt.plot(x, y_pred_l, label="Likelihood prediction")
plt.ylim([-0.1, 1.1])
plt.legend(loc=0, numpoints=1)
plt.title("Exito o Falla en lanzamiento de Challenger")
plt.xlabel(r"T [${}^o C$]")
plt.ylabel(r"$y$")
plt.show()
In [20]:
import numpy as np
from sklearn import datasets
# Loading the data
iris = datasets.load_iris()
X = iris.data
Y = iris.target
print iris.target_names
# Print data and labels
for x, y in zip(X,Y):
print x, y
In [21]:
import numpy as np
from sklearn import datasets
# Loading the data
iris = datasets.load_iris()
names = iris.target_names
print names
X = iris.data
Y = np.array(iris.target==0, int)
# Print data and labels
for x, y in zip(X,Y):
print x, y
Para aplicar el algoritmo, utilizando el algoritmo Logistic Regression de la librería sklearn, requerimos un código como el siguiente:
In [27]:
import numpy as np
from sklearn import datasets
from sklearn.linear_model import LogisticRegression
# Loading the data
iris = datasets.load_iris()
names = iris.target_names
X = iris.data
Y = np.array(iris.target==0, int)
# Fitting the model
Logit = LogisticRegression()
Logit.fit(X,Y)
# Obtain the coefficients
print Logit.intercept_, Logit.coef_
# Predicting values
Y_pred = Logit.predict(X)
#x = X.mean(axis=0)
#Y_pred_mean = Logit.predict(x)
#print x, Y_pred_mean
In [28]:
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(Y, Y_pred)
print cm
¡Nuestra clasificación es perfecta!