Regressão Linear Simples com Descida de Gradiente

A Regressão Linear é um método de aprendizagem de máquina supervisionada para prever rótulos contínuos em dados que seguem uma tendência linear.

Por exemplo, a seguir temos representados num gráfico os preços de casas (eixo Y) de acordo com o seu tamanho (eixo X). Os preços são mostrados em dólares americanos e estão divididos por mil. Os tamanhos estão em "square feet" e também estão divididos por mil. Esses dados foram coletados em 1977, assim, os preços das casas são muito inferiores aos preços atuais.

Tais dados são os dados de treinamento. Eles foram extraídos do link: http://people.sc.fsu.edu/~jburkardt/datasets/regression/x27.txt

In [1]:
from bokeh.io import show, output_notebook
from bokeh.plotting import figure
output_notebook()

sizes = [0.998, 1.500, 1.175, 1.232, 1.121, 0.988, 1.240, 1.501, 1.225, 1.552, 0.975, 1.121, 1.020,
        1.501, 1.664, 1.488, 1.376, 1.500, 1.256, 1.690, 1.820, 1.652, 1.777, 1.504, 1.831, 1.200 ]
prices = [25.9, 29.5, 27.9, 25.9, 29.9, 29.9, 30.9, 28.9, 35.9, 31.5, 31.0, 30.9, 30.0, 28.9, 36.9,
          41.9, 40.5, 43.9, 37.5, 37.9, 44.5, 37.9, 38.9, 36.9, 45.8, 41.0]

p = figure(plot_width=700, plot_height=250)
p.circle(sizes, prices)
show(p, notebook_handle=True)
Loading BokehJS ...
Out[1]:

<Bokeh Notebook handle for In[1]>

Utilizando os dados, como saber qual deve ser o preço de uma casa com o tamanho 1400 square feet?

Uma das formas de responder essa pergunta é traçar uma reta seguindo a tendência linear dos dados de treinamento e usá-la como guia para prever o valor de um ponto que não se encontra no gráfico.

In [2]:
p.line([1,2], [29.3015, 43.3078], line_width=2)
show(p, notebook_handle=True)
Out[2]:

<Bokeh Notebook handle for In[2]>

Visualmente podemos perceber que o preço apontado pela reta para uma casa de 1400 square feet é aproximadamente 35 mil.

Uma reta que segue a tendência linear de um conjunto de pontos em duas dimensões pode ser gerada usando Regressão Linear Simples. A Regressão Linear Simples modela o relacionamento entre as variáveis através de uma reta.

Na matemática, uma reta pode ser definida como uma função linear: $f(x) = w_1x + w_0$, onde $w_0$ e $w_1$ são os coeficientes da função. $w_0$ especifica o deslocamento vertical da reta e $w_1$ a sua inclinação. No gráfico acima a função da reta é: $f(x) = 14.0063x + 15.2952$.

A reta ideal encontrada pela Regressão Linear Simples é a reta na qual a média da distância de todos os pontos para a reta é a menor possível. Ou seja, a reta que minimiza o valor de $\frac{1}{2m} \sum_{i=1}^{m} (f(x_i) - y_i)^2$, onde m é a quantidade de pontos nos dados de treinamento, e ($x_i$, $y_i$) define um ponto.

Analisando o somatório, temos que $(f(x_i) - y_i)$ mede a distância entre um ponto $(x_i, y_i)$ e a reta, já que $y_i$ é o valor do ponto no eixo Y e $f(x_i)$ é o valor estimado pela reta. A diferença é elevada ao quadrado para que as distâncias sejam sempre positivas, não fazendo diferença se o ponto está acima ou abaixo da reta. Por fim, a média é calculada dividindo o somatório por $2m$. O número 2 aparece para facilitar cálculos que serão realizados em passos seguintes.

Podemos escrever o somatório anterior como uma função de $w_0$ e $w_1$. Assim, temos:

\begin{equation*} J(w_0, w_1) = \frac{1}{2m} \sum_{i=1}^{m} (w_1x + w_0 - y_i)^2 \end{equation*}

Essa função é chamada de Erro Médio Quadrático.

Vamos definir uma função em Python para calcular o Erro Médio Quadrático:

In [3]:
def mean_squared_error(xs, ys, w0, w1):
    total = 0

    for i in range(len(xs)):
        x, y = xs[i], ys[i]
        fx = w1*x + w0
        error = (fx - y)**2
        total += error

    return total/(2*len(xs))

Calculando o Erro Médio Quadrático para a função da reta exibida anteriormente, onde $w_0 = 15.2952$ e $w_1 = 14.0063$, temos:

In [4]:
print("O Erro Médio Quadrático é: {}".format(mean_squared_error(sizes, prices, 15.2952, 14.0063)))
O Erro Médio Quadrático é: 10.553925770613349

Para encontrar o valor mínimo de $J(w_0, w_1)$ utilizaremos o algoritmo da descida de gradiente.

O algoritmo inicia com dois valores aleatórios para $w_0$ e $w_1$ e em seguida, passo a passo, ele ajusta os valores de $w_0$ e $w_1$ para que $J(w_0, w_1)$ seja minimizada. O tamanho desses ajustes em cada passo é determinado por um parâmetro $\alpha$ chamado de taxa de aprendizagem. Se $\alpha$ for muito pequeno o algoritmo pode demorar a convergir, se $\alpha$ for um valor alto os ajustes podem ser grandes de uma forma que o algoritmo nunca consiga convergir.

Para descobrir a direção que $w_0$ e $w_1$ devem ser ajustados para diminuir o valor de $J(w_0, w_1)$ vamos calcular o gradiente da função, ou seja, vamos calcular as derivadas parciais de $J(w_0, w_1)$.

  • $\frac{\partial}{\partial w_0} J(w_0, w_1) = \frac{1}{m} \sum_{i=1}^{m} (f(x_i) - y_i)$
  • $\frac{\partial}{\partial w_1} J(w_0, w_1) = \frac{1}{m} \sum_{i=1}^{m} ((f(x_i) - y_i)x_i)$

As derivadas parciais são calculadas pelas funções abaixo:

In [5]:
def d_w0(xs, ys, w0, w1):
    total = 0

    for i in range(len(xs)):
        x, y = xs[i], ys[i]
        fx = w1*x + w0
        error = (fx - y)
        total += error

    return total/(len(xs))


def d_w1(xs, ys, w0, w1):
    total = 0

    for i in range(len(xs)):
        x, y = xs[i], ys[i]
        fx = w1*x + w0
        error = (fx - y)
        total += error*x

    return total/(len(xs))

Com isso, o algoritmo se resume a repetir até a convergência os seguintes passos:

  • $w_0 := w_0 - \alpha \frac{\partial}{\partial w_0} J(w_0, w_1)$
  • $w_1 := w_1 - \alpha \frac{\partial}{\partial w_1} J(w_0, w_1)$

É importante que as atualizações de $w_0$ e $w_1$ sejam feitas de forma simultânea.

A cada iteração gravaremos os valores de $w_0$, $w_1$ e do Erro Médio Quadrático para analisarmos as mudanças na reta posteriormente. Além disso, o algoritmo terminará quando o Erro Médio Quadrático não diminuir pelo menos 0.00001.

In [6]:
import random

w0 = random.uniform(-2, 2)
w1 = random.uniform(-2, 2)
alpha = 0.1

coefficients = [(w0, w1)]
errors = [mean_squared_error(sizes, prices, w0, w1)]

while True:
    w0i = w0 - alpha*d_w0(sizes, prices, w0, w1)
    w1i = w1 - alpha*d_w1(sizes, prices, w0, w1)

    w0 = w0i
    w1 = w1i
    
    coefficients.append((w0, w1))
    errors.append(mean_squared_error(sizes, prices, w0, w1))
    
    if (errors[-2] - errors[-1] < 0.00001):
        break

print("Mean squared error: {}".format(mean_squared_error(sizes, prices, w0, w1)))
print("f(x) = {}x + {}".format(w1, w0))
Mean squared error: 10.556045385998793
f(x) = 14.251632539788755x + 14.948262924183263

Usando a função encontrada, vamos visualizar a reta iniciando em $x=1$ e terminando em $x=2$:

In [7]:
p_final = figure(plot_width=700, plot_height=250)
p_final.circle(sizes, prices)
p_final.line([1,2], [w1 + w0, w1*2 + w0], line_width=2)
show(p_final, notebook_handle=True)
Out[7]:

<Bokeh Notebook handle for In[7]>

Como gravamos a evolução dos valores de $w_0$ e $w_1$ podemos plotar gráficos para visualizarmos quais foram as retas geradas nas iterações. A função a seguir exibe as retas nas 12 primeiras iterações.

In [8]:
from bokeh.layouts import gridplot

plots_grid = [[]]
for i in range(12):
    w0, w1 = coefficients[i]
    
    p = figure(plot_width=250, plot_height=200)
    p.circle(sizes, prices)
    p.line([1,2], [w1 + w0, w1*2 + w0], line_width=2)
    
    if (len(plots_grid[-1]) == 3):
        plots_grid.append([])
    
    plots_grid[-1].append(p)
    
p = gridplot(plots_grid)
show(p)

Logo nas primeiras iterações a reta já é muito semelhante ao resultado final.

Outra forma de analisar a evolução da solução ao longo das iterações é através dos valores dos erros. Vamos plotar uma curva com todos os seus valores.

In [9]:
error_points = [[], []]
for i, error in enumerate(errors):
    error_points[0].append(i)
    error_points[1].append(error)

error_chart = figure(plot_width=800, plot_height=300)
error_chart.line(error_points[0], error_points[1], line_width=2)
show(error_chart)

O Erro Médio Quadrático cai rapidamente nas primeiras iterações, confirmando o fato das retas nas primeiras iterações já serem muito parecidas com o resultado final.