Regressão Polinomial Com Descida de Gradiente

Nos notebooks Regressão Linear Simples com Descida de Gradiente e Regressão Linear Múltipla com Descida de Gradiente foi mostrado como ajustar funções lineares a um conjunto de dados para realizar previsões.

Em ambos os casos, uma função só se adequava bem aos dados de treinamento se eles seguissem uma tendência linear. Entretanto, para dados que não seguiam essa tendência, a função não se adequava bem, o que acabaria gerando previsões ruins. Assim, para dados não lineares é possível ajustar outros tipos de funções que não sejam lineares através da Regressão Polinomial, gerando curvas ou superfícies mais complexas.

Continuaremos utilizando o exemplo dos preços das casas de acordo com o seu tamanho. Abaixo temos o gráfico representando esses dados.

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]>

A Regressão Polinomial funciona de maneira muito parecida com a regressão linear, a principal diferença é o tipo da função que é gerada. A regressão linear gera retas ou superfícies representadas por funções do tipo: $f(x) = w_0x_0 + w_1x_1 + ... + x_2x_n$. Por outro lado, a Regressão Polinomial pode gerar curvas e superfícies dos mais variados tipos, por exemplo: $f(x) = w_0x_0 + w_1x_1^2$ ou $f(x) = w_0 + w_1x_1 + w_2x_2^3$.

Vamos fazer um exemplo tentando ajustar uma função do tipo: $f(x) = w_0 + w_1x^2$ aos dados de treinamento.

O código abaixo recebe os valores de x, $w_0$ e $w_1$ e calcula o resultado da função.

In [2]:
def f(x, w0, w1):
    return w1*(x**2) + w0

A função de erro médio quadrático continua a mesma. Ela é calculada através do código abaixo:

In [3]:
def square_error(X, Y, w0, w1):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1)
        error = (yi - y)**2

        total += error

    return total/(2*len(X))

Mais uma vez utilizaremos a descida de gradiente para encontrar a função que melhor se ajusta aos dados de treinamento. Com isso, durante o algoritmo, será necessário calcular as derivadas da função de erro médio quadrático.

O código abaixo calcula a derivada parcial em relação a $w_0$

In [4]:
def d_error_w0(X, Y, w0, w1):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1)
        error = (yi - y)

        total += error

    return total/(len(X))

E a função abaixo calcula a derivada parcial em relação a $w_1$

In [5]:
def d_error_w1(X, Y, w0, w1):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1)
        error = (yi - y)

        total += error*x**2

    return total/(len(X))

Com a função e as derivadas parciais em mãos, aplicaremos o algoritmo de descida de gradiente. No código abaixo, os pesos são atualizados 1000 vezes.

In [6]:
w = [0, 0]
learning_rate = 0.1

for i in range(1000):        
    w0i = w[0] - learning_rate*d_error_w0(sizes, prices, w[0], w[1])
    w1i = w[1] - learning_rate*d_error_w1(sizes, prices, w[0], w[1])

    w[0] = w0i
    w[1] = w1i

print("Função gerada: {} + {}x^2".format(w[0], w[1]))
Função gerada: 24.615774400324003 + 5.070623894139969x^2

No gráfico a seguir observarmos que a curva encontrada não é mais uma reta.

In [11]:
import numpy as np

x_points = np.linspace(0.8, 2, 50)
y_points = [f(i, w[0], w[1]) for i in x_points]

p = figure(plot_width=700, plot_height=250)
p.circle(sizes, prices)
p.line(x_points, y_points, line_width=2)
show(p, notebook_handle=True)
Out[11]:

<Bokeh Notebook handle for In[11]>

Vamos fazer um segundo exemplo, agora tentando ajustar a função: $f(x) = w_3x^3 + w_2x^2 + w_1x + w_0$

In [12]:
def f(x, w0, w1, w2, w3):
    return w3*(x**3) + w2*(x**2) + w1*x + w0

Nesse exemplo, precisamos calcular quatro derivadas parciais, pois temos 4 coeficientes ($w_0, w_1, w_2, w_3$)

In [13]:
def d_error_w0(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error

    return total/(len(X))

def d_error_w1(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error*x

    return total/(len(X))

def d_error_w2(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error*x**2

    return total/(len(X))

def d_error_w3(X, Y, w0, w1, w2, w3):
    total = 0

    for i in range(len(X)):
        x, y = X[i], Y[i]
        yi = f(x, w0, w1, w2, w3)
        error = (yi - y)

        total += error*x**3

    return total/(len(X))

Por fim, realizando a descida de gradiente, encontramos os coeficientes.

In [14]:
w = [0, 0, 0, 0]
learning_rate = 0.1

for i in range(1000):        
    w0i = w[0] - learning_rate*d_error_w0(sizes, prices, w[0], w[1], w[2], w[3])
    w1i = w[1] - learning_rate*d_error_w1(sizes, prices, w[0], w[1], w[2], w[3])
    w2i = w[2] - learning_rate*d_error_w2(sizes, prices, w[0], w[1], w[2], w[3])
    w3i = w[3] - learning_rate*d_error_w3(sizes, prices, w[0], w[1], w[2], w[3])

    w[0] = w0i
    w[1] = w1i
    w[2] = w2i
    w[3] = w3i
    
print("Função gerada: {} + {}x^2 + {}x^3 + {}x^4".format(w[0], w[1], w[2], w[3]))
Função gerada: 15.796264315235701 + 10.422723031103592x^2 + 3.89992652251128x^3 + -1.1240198966981076x^4
In [18]:
x_points = np.linspace(0, 3, 50)
y_points = [f(i, w[0], w[1], w[2], w[3]) for i in x_points]

p = figure(plot_width=700, plot_height=250)
p.circle(sizes, prices)
p.line(x_points, y_points, line_width=2)
show(p, notebook_handle=True)
Out[18]:

<Bokeh Notebook handle for In[18]>