Machine Learning: Regresión Lineal

Ecuación Normal y Descenso del Gradiente


Contenido



  1. Introducción
  2. Paquetes Python
  3. Prep. de Datos
  4. Variables
  5. Modelo
  6. Ecuación Normal
  7. D. Gradiente
  8. Scikit Learn
  9. Modelo
  10. Visualizaciones
  11. Predicciones
  12. Multivariada
  13. Redes Sociales

Introducción

En este proyecto haremos un recorrido muy completo por la regresión lineal simple y multivariada. Utilizaremos la ecuación normal y el descenso del gradiente para encontrar los parámetros buscados y trabajaremos sobre el conjunto de datos conocido como Boston Housing.
El conjunto de datos está incluido en el paquete scikit learn (Machine Learning in Python).

Paquetes Python

In [2]:
import datetime
import pandas as pd
import numpy as np
import random

from sklearn.datasets import load_boston
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm

import holoviews as hv
from holoviews import opts
import plotly
import plotly.graph_objects as go
#hv.extension('plotly')

from IPython.display import HTML
from IPython.core.display import display

from bokeh import plotting
from bokeh.io import output_notebook, show
from bokeh.plotting import figure, output_file, show
from bokeh.models import HoverTool
import seaborn as sns
Loading BokehJS ...

Preparación de Datos

In [4]:
boston = load_boston()
bostonH = pd.DataFrame(boston.data)
bostonH.columns = boston.feature_names
bostonH['PRECIO'] = boston.target
bostonH['X0'] = 1
bostonH = bostonH[[
    'X0', 'CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
    'TAX', 'PTRATIO', 'B', 'LSTAT', 'PRECIO'
]]
#bostonH.head().iloc[:, 0:5]

Descripción de Datos

  • CRIM: per capita crime rate by town.
  • ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
  • INDUS: proportion of non-retail business acres per town.
  • CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).
  • NOX: nitrogen oxides concentration (parts per 10 million).
  • RM: average number of rooms per dwelling.
  • AGE: proportion of owner-occupied units built prior to 1940.
  • DIS: weighted mean of distances to five Boston employment centres.
  • RAD: index of accessibility to radial highways.
  • TAX: full-value property-tax rate per $10,000.
  • PRATIO: pupil-teacher ratio by town.
  • BLACK: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.
  • LSTAT: lower status of the population (percent).
  • MEDV (PRECIO): median value of owner-occupied homes in $1000s.

Relación Entre Variables

In [5]:
sns.pairplot(bostonH.iloc[:][['RM', 'AGE', 'LSTAT', 'PRECIO']])
Out[5]:
<seaborn.axisgrid.PairGrid at 0x7f773c1e8fa0>

Para la regresión lineal simple, elegimos como variable independiente a "RM" (número de habitaciones por casa) y como objetivo "PRECIO" (precio de la casa) es decir, queremos elaborar un modelo para predecir el precio de la casa en función del número de habitaciones.
Realizamos un scatter plot para revisar la relación a detalle.

In [6]:
def vis_1():
    ''' 
    Autor: Francisco Jaramillo
    https://www.pakin.lat/
    
    '''

    source = plotting.ColumnDataSource(bostonH)

    p = figure(
        title="",
        x_axis_label='Número de Habitaciones (RM)',
        y_axis_label='Precio',
        plot_width=900,
        plot_height=400,
        sizing_mode='scale_width',
    )

    plot1 = p.x(x='RM', y='PRECIO', size=9, source=source, color='red')

    p.add_tools(
        HoverTool(
            renderers=[plot1],
            tooltips=[
                ("Cant. habitaciones", "@RM"),
                ("Precio", "@PRECIO"),
            ],
        ))

    p.toolbar.active_drag = None
    p.left[0].formatter.use_scientific = False
    p.title.text_font_size = '3vw'
    p.xaxis.axis_label_text_font_size = "2vw"
    p.yaxis.axis_label_text_font_size = "2vw"
    p.xaxis.major_label_text_font_size = "1.3vw"
    p.yaxis.major_label_text_font_size = "1.8vw"

    show(p)

Modelo

Variables de entrada (input feature: número de habitaciones)


\begin{equation} \left(x_{1}, x_{2},...x_{i}\right); \ i = 1,2,...,m \end{equation}

Variables salida (target: precio)
\begin{equation} \left(y_{1}, y_{2},...y_{i}\right); \ i = 1,2,...,m \end{equation}

Conjunto de entrenamiento
\begin{equation} \left(x_{i} , y_{i}\right); \ i = 1,2,3,...,m \end{equation}

$X$ e $Y$ son los espacios de entrada y salida respectivamente y $ X, Y \in \mathbb{R} $ .
Cuando la variable objetivo, la que queremos predecir es contínua, entonces estamos frente a un problema de regresión.
En la regresión lineal simple, queremos que dado un conjunto de entrenamiento $ \left( x_{i} , y_{i}\right) $ que servirá para entrenar una función $ h: X \rightarrow Y $, la función $ h(x) $ sea un buen predictor de los elementos de $ Y $.
Donde: $ h_{\theta}\left(x\right) = \theta_{0} + \theta_{1}x$ es una hipótesis.

Tomemos como ejemplo la función hipótesis:


\begin{equation} h_{\theta}\left(x\right) = 5x + 1 \end{equation}

Queremos minimizar todas las distancias verticales entre los puntos de la recta y los puntos del conjunto de entrenamiento; esto implica encontrar la recta $ h_{\theta}\left(x\right) $ que minimice estas distancias y esto a su vez implica encontrar el valor de los parámetros ($ \theta_{0}, \theta_{1} $).

In [8]:
xs = list(range(3, 10))
ys = [1 + 5 * x for x in xs]
In [9]:
def vis_2():
    ''' 
    Autor: Francisco Jaramillo
    https://www.pakin.lat/
    
    '''

    source = plotting.ColumnDataSource(bostonH)

    p = figure(
        title="",
        x_axis_label='Número de Habitaciones (RM)',
        y_axis_label='Precio',
        plot_width=900,
        plot_height=400,
        sizing_mode='scale_width',
    )

    plot1 = p.x(x='RM', y='PRECIO', size=9, source=source, color='red')

    plot2 = p.line(x=xs, y=ys, line_width=2, color='black')

    p.segment(x0=7.393,
              x1=7.393,
              y0=17.8,
              y1=37.965,
              color="blue",
              line_width=1)
    p.segment(x0=6.144,
              x1=6.144,
              y0=36.2,
              y1=31.72,
              color="blue",
              line_width=1)

    p.add_tools(
        HoverTool(
            renderers=[plot1],
            tooltips=[
                ("Cant. habitaciones", "@RM"),
                ("Precio", "@PRECIO"),
            ],
        ))

    p.toolbar.active_drag = None
    p.left[0].formatter.use_scientific = False
    p.title.text_font_size = '3vw'
    p.xaxis.axis_label_text_font_size = "2vw"
    p.yaxis.axis_label_text_font_size = "2vw"
    p.xaxis.major_label_text_font_size = "1.3vw"
    p.yaxis.major_label_text_font_size = "1.8vw"

    show(p)

Ecuación Normal

La ecuación normal:

$$\theta = (X^{T}X)^{-1} X^{T}Y$$

Nos da la solución óptima a nuestro problema, en donde:

$$\theta = \begin{bmatrix} \theta_{0} \\ \theta_{1} \end{bmatrix}$$
$$X = \begin{bmatrix} 1 & x_{1} \\ 1 & x_{2} \\ 1 & x_{3} \\ . & . \\ . & . \\ 1 & x_{m} \end{bmatrix}$$
$$Y = \begin{bmatrix} y_{1} \\ y_{2} \\ y_{3} \\ . \\ . \\ y_{m} \end{bmatrix}$$

Implementación de la Ecuación Normal

In [11]:
m = len(bostonH.iloc[:]['X0'])

X = np.matrix(bostonH.iloc[:][['X0', 'RM']])

Y = np.matrix(bostonH.iloc[:][['PRECIO']])
In [12]:
en = np.linalg.inv(np.transpose(X).dot(X)).dot(np.transpose(X)).dot(Y)
en
Out[12]:
matrix([[-34.67062078],
        [  9.10210898]])

Descenso del Gradiente

En el descenso del gradiente, si utilizamos todos los puntos del conjunto de entrenamiento (m), entonces es llamado Descenso del Gradiente por Lote (Batch Gradient Descent).
Si utilizamos una muestra aleatoria de tamaño n < m, entonces es llamado Descenso del Gradiente por Mini Lote (Mini Batch Gradient Descent).
Si utilizamos sólo un punto (n = 1) elegido aleatoriamente en cada iteración, entonces es llamado Descenso del Gradiente Estocástico (Stochastic Gradient Descent).


Función de costo (función de error cuadrático)


\begin{equation} \begin{split} J & = \frac{1}{2m}\sum_{i = 1}^{m} {\left(\hat{y}_{i} - y_{i}\right)}^{2} = \\ & J = \frac{1}{2m}\sum_{i = 1}^{m} {\left[h_{\theta}(x_{i}) - y_{i}\right]}^{2} = \\ & J(\theta_{0}, \theta_{1}) = \frac{1}{2m}\sum_{i = 1}^{m} {\left[\left(\theta_{0} + \theta_{1}x_{i}\right) - y_{i}\right]}^{2} \end{split} \end{equation}

Derivadas de la Función de Costo


\begin{equation} \frac{\partial}{\partial \theta_{0}}J\left(\theta_{0}, \theta_{1}\right) =\frac{1}{m}\sum_{i = 1}^{m} {\left[h_{\theta}(x_{i}) - y_{i}\right]} \end{equation}
\begin{equation} \frac{\partial}{\partial \theta_{1}}J(\theta_{0}, \theta_{1}) =\frac{1}{m}\sum_{i = 1}^{m} {\left[h_{\theta}\left(x_{i}\right) - y_{i}\right]}x_{i} \end{equation}

Ecuaciones para Descenso del Gradiente


\begin{equation} \theta_0:= \theta_0 - \alpha\frac{1}{m}\sum_{i = 1}^{m} {\left[\left(\theta_{0} + \theta_{1}x_{i}\right) - y_{i}\right]} \end{equation}
\begin{equation} \theta_1:= \theta_1 - \alpha\frac{1}{m}\sum_{i = 1}^{m} {\left[\left(\theta_{0} + \theta_{1}x_{i}\right) - y_{i}\right]}x_{i} \end{equation}

Implementación del Descenso el Gradiente

Función de Costo

In [13]:
def funCosto(x, y, m, theta):
    ''' 
    Autor: Francisco Jaramillo
    https://www.pakin.lat/
    
    '''
    '''Devuelve el valor de la función de Costo
       x: variable independiente (matriz)
       y: variable dependiente (matriz)
       m: número de muestras (escalar)
    '''
    return np.sum(np.square((x.dot(theta) - y))) / (2 * m)

Implementación del Descenso de Gradiente por Lote

In [14]:
def DescGrad(x, y, m, alpha, iterations):
    
    ''' 
    Autor: Francisco Jaramillo
    https://www.pakin.lat/
    
    '''
    
    '''Implementa el descenso del gradiente por lote.
       Calcula y guarda el valor de la función de costo J en cada iteración.
       Guarda una serie del número de iteraciones n_iter.
       Devuelve los parámteros theta y un dataframe con el número iteraciones
       y los valores de la función costo.
    '''
    
    J = []
    n_iter = []
    theta = np.zeros([2, 1])

    for _ in range(iterations):

        theta[0][0] = theta[0][0] - alpha * np.sum((x.dot(theta) - y)) / (m)

        theta[1][0] = theta[1][0] - alpha * np.sum(
            np.transpose(
                (x.dot(theta) - y)) * (np.transpose(np.transpose(x)[1]))) / (m)
        
        J.append(funCosto(X, Y, m, theta))
        n_iter.append(_)
    
    
    df = pd.DataFrame([n_iter, J])
    df = df.transpose()
    #DGEst.reset_index(inplace=True)
    df.columns = ['N_ITER', 'J']
    df['N_ITER'] = df['N_ITER'] + 1

    return theta, df
In [15]:
dgL = DescGrad(X, Y, m, alpha = 0.05, iterations = 8500)
dgL[0]
Out[15]:
array([[-34.50660814],
       [  9.07634126]])
In [16]:
def vis_3():
    
    ''' 
    Autor: Francisco Jaramillo
    https://www.pakin.lat/
    
    '''
    
    source = plotting.ColumnDataSource(dgL[1])

    p = figure(
        title="",
        x_axis_label='Número de Iteraciones',
        y_axis_label='Función de Costo',
        plot_width=900,
        plot_height=400,
        sizing_mode='scale_width',
    )

    plot1 = p.line(x='N_ITER', y='J', width=2, source=source, color='red')

    #plot2 = p.line(x='N_ITER', y='J2', width=4, source=source, color='blue')

    #plot3 = p.line(x='N_ITER', y='J3', width=8, source=source, color='green')

    p.add_tools(
        HoverTool(
            renderers=[plot1],
            tooltips=[
                ("Número Iteraciones", "@N_ITER"),
                ("Función de Costo", "@J"),
            ],
        ))

    p.toolbar.active_drag = None
    p.left[0].formatter.use_scientific = False
    p.title.text_font_size = '3vw'
    p.xaxis.axis_label_text_font_size = "2vw"
    p.yaxis.axis_label_text_font_size = "2vw"
    p.xaxis.major_label_text_font_size = "1.3vw"
    p.yaxis.major_label_text_font_size = "1.8vw"

    show(p)

Implementación del Descenso de Gradiente por Mini Lote

In [18]:
def muesAl(r):
    
    ''' 
    Autor: Francisco Jaramillo
    https://www.pakin.lat/
    
    '''
    
    '''Devuelve una muestra aleatoria del conjunto de entrenamiento
       n: tamaño de la muestra en función de la proporción r
       inds: toma n índices aleatorios
    '''
    
    n = int(np.ceil(bostonH.index.size*r))

    inds = random.sample(range(bostonH.index.size), n)

    x = np.matrix(bostonH.iloc[inds][['X0', 'RM']])

    y = np.matrix(bostonH.iloc[inds][['PRECIO']])
    
    return x, y, n
In [19]:
X_ml, Y_ml, n = muesAl(r=0.7)

dgML = DescGrad(X_ml, Y_ml, n, alpha=0.05, iterations=15000)

dgML[0]
Out[19]:
array([[-29.61871161],
       [  8.2363519 ]])
In [20]:
def vis_4():
    ''' 
    Autor: Francisco Jaramillo
    https://www.pakin.lat/
    
    '''

    source = plotting.ColumnDataSource(dgML[1])

    p = figure(
        title="",
        x_axis_label='Número de Iteraciones',
        y_axis_label='Función de Costo',
        plot_width=900,
        plot_height=400,
        sizing_mode='scale_width',
    )

    plot1 = p.line(x='N_ITER', y='J', width=2, source=source, color='blue')

    p.add_tools(
        HoverTool(
            renderers=[plot1],
            tooltips=[
                ("Número Iteraciones", "@N_ITER"),
                ("Función de Costo", "@J"),
            ],
        ))

    p.toolbar.active_drag = None
    p.left[0].formatter.use_scientific = False
    p.title.text_font_size = '3vw'
    p.xaxis.axis_label_text_font_size = "2vw"
    p.yaxis.axis_label_text_font_size = "2vw"
    p.xaxis.major_label_text_font_size = "1.3vw"
    p.yaxis.major_label_text_font_size = "1.8vw"

    show(p)