#!/usr/bin/env python
# coding: utf-8

# # Regressão Linear Múltipla com Descida de Gradiente
# 
# No notebook [Regressão Linear Simples com Descida de Gradiente](https://nbviewer.jupyter.org/github/yurimalheiros/ai-notebooks/blob/master/ml/gdlinearregression.ipynb) foi mostrado como funciona a regressão linear
# para prever uma resposta (preço de uma casa) a partir de uma variável de entrada (tamanho de uma casa).
# 
# Entretanto, existem diversos casos nos quais temos duas ou mais variáveis de entrada para prever uma resposta. Para isso precisamos de uma abordagem que estenda a Regressão Linear Simples para suportar múltiplas variáveis de entrada.

# Neste notebook, além dos dados dos tamanhos das casas para prever os preços, usaremos também o número de quartos. A seguir temos um gráfico tridimensional mostrando cada casa dos dados de treinamento como um ponto de acordo com as três características (tamanho, número de quartos e preço).

# In[1]:


import plotly
import plotly.graph_objs as go
plotly.offline.init_notebook_mode()

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]

bedrooms = [4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 2, 3, 4, 3, 3, 3, 3, 3, 4, 3, 4, 3, 4, 3]

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]

plot = go.Scatter3d(x=sizes, y=bedrooms, z=prices, mode='markers')
layout = go.Layout(width=900, height=400,
                   margin=go.layout.Margin(l=10, r=10, b=10, t=10, pad=10),
                   scene=dict(xaxis=dict(title="tamanho"), yaxis=dict(title="quartos"), zaxis=dict(title="preço")))

data = [plot]
plotly.offline.iplot({'data': data, 'layout': layout})


# A Regressão Linear Múltipla é uma generalização da Regressão Linear Simples, que podemos usar para qualquer quantidade de variáveis.
# 
# Diferente da Regressão Linear Simples, na qual foi gerada uma reta para prever o valor de uma casa de acordo com um tamanho, na Regressão Linear Múltipla, como temos três ou mais dimensões, geraremos um plano (ou hiperplano).

# Para entender os próximos passos definiremos algumas notações.
# 
# Cada característica de um exemplo dos dados de treinamento será representada por $x_i$. Nesse caso, $x_1$ é o tamanho da casa e $x_2$ é a quantidade de quartos. Na regressão múltipla podemos ter $N$ características.
# 
# 
# Com isso, a função que representa o plano para esse exemplo é:
# 
# \begin{equation*}
# f(x) = w_0 + w_1x_1 + w_2x_2
# \end{equation*}
# 
# Generalizando a função para exemplos com $N$ características, temos:
# 
# \begin{equation*}
# f(x) = w_0 + w_1x_1 + ... + w_nx_n
# \end{equation*}

# O vetor $x$ vai representar todas as características de uma casa, mas com uma ressalva, ele terá $N+1$ elementos, sendo o primeiro sempre $x_0 = 1$. Isto nos ajudará em cálculos posteriores. Assim:
# 
# $x = \begin{bmatrix}
#      x_0 \\
#      x_1 \\
#      ... \\
#      x_n
#      \end{bmatrix}$
#      
# Além disso, temos o vetor $w$ que contém todos os coeficientes da função do plano.
# 
# $w = \begin{bmatrix}
#      w_0 \\
#      w_1 \\
#      ... \\
#      w_n
#      \end{bmatrix}$
#      
# 
# Representando os dados das casas e os coeficientes como vetores, podemos simplificar a função escrevendo da seguinte forma:
# 
# $f(x) = w^Tx$
# 
# Onde $w^T$ é a matriz transposta de $w$.
# 
# Efetuando a multiplicação $w^Tx$ temos: $w_0x_0 + w_1x_1 + ... + w_nx_n$

# Abaixo temos a implementação para calcular o resultado da função de acordo com os vetores $x$ e $w$.

# In[2]:


def f(x, w):
    result = 0
    for i in range(len(x)):
        result += x[i]*w[i]
    
    return result

print(f([1,3,-1], [1,1,2]))


# Além disso, precisamos especificar a qual exemplo pertence uma característica.
# Para isso usaremos a notação $x_j^{(i)}$ para representar a característica $x_j$ da i-ésima casa.
# Por exemplo, o tamanho da quinta casa nos dados de treinamento é representado por $x_1^{(5)}$, já a quantidade de quartos da terceira casa é representada por $x_2^{(3)}$. Dessa forma, o vetor de características da i-ésima casa será representado por $x^{(i)}$.

# Como na Regressão Linear Simples, a função resultante ideal na Regressão Linear Múltipla é a que minimiza a função de erro médio quadrático $\frac{1}{2m} \sum_{i=1}^{m} (f(x^{(i)}) - y^{(i)})^2$, onde $m$ é a quantidade de exemplos nos dados de treinamento.
# 
# Escrevendo ela em função de $w$, temos:
# 
# \begin{equation*}
# J(w) = \frac{1}{2m} \sum_{i=1}^{m} (w^Tx^{(i)} - y^{(i)})^2
# \end{equation*}
# 
# O código para esse cálculo pode ser visto abaixo. Na função implementada são passados todos os dados dos tamanhos, números de quartos e preços das casas, além do vetor $w$ com os coeficientes da função.

# In[3]:


def squared_error(sizes, bedrooms, prices, w):
    total = 0

    for i in range(len(sizes)):
        x1, x2, y = sizes[i], bedrooms[i], prices[i]
        yi = f([x1, x2], w)
        error = (yi - y)**2

        total += error

    return total/(2*len(sizes))


# Novamente, para encontrar o valor mínimo de $J(w)$ utilizaremos o algoritmo da descida de gradiente.
# 
# Na Regressão Linear Simples tínhamos duas regras para atualizar $w_0$ e $w_1$ utilizando derivadas parciais.
# Entretanto, podemos generalizar o cálculo das derivadas parciais para qualquer quantidade de coeficientes.
# 
# \begin{equation*}
# \frac{\partial}{\partial w_j} J(w) = \frac{1}{m} \sum_{i=1}^{m} (f(x^{(i)}) - y^{(i)})x_j^{(i)}
# \end{equation*}

# Como convencionamos que $x_0^{(i)}$ é sempre igual a 1, então $\frac{\partial}{\partial w_0} J(w)$ vai ser igual ao que calculamos na Regressão Linear Simples, assim como $\frac{\partial}{\partial w_1} J(w)$.

# As derivadas podem ser calculadas usando o código a seguir. Na função são passados todos os dados dos tamanhos, números de quartos e preços das casas, o vetor $w$ e o índice $j$ para especificar qual derivada parcial está sendo calculada.

# In[4]:


def d_error(sizes, bedrooms, prices, w, j_index):
    total = 0

    for i in range(len(sizes)):
        x1, x2, y = sizes[i], bedrooms[i], prices[i]
        xij = [1, sizes[i], bedrooms[i]][j_index]
        yi = f([x1, x2], w)
        error = (yi - y)

        total += error*xij

    return total/(len(sizes))


# Com isso, podemos resumir a atualização dos pesos nas iterações do algoritmo de descida de gradiente para Regressão Múltipla assim:
# 
# \begin{equation*}
# w_j := w_j - \alpha \frac{\partial}{\partial w_j} J(w)
# \end{equation*}
# 
# ou seja
# 
# \begin{equation*}
# w_j := w_j - \alpha \frac{1}{m} \sum_{i=1}^{m} (f(x^{(i)}) - y^{(i)})x_j^{(i)}
# \end{equation*}
# 
# onde $j$ varia de 0 até $N$.

# Vamos agora fazer atualizações dos coeficientes até o erro médio quadrático não diminuir pelo menos 0.000001.

# In[5]:


import random

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

coefficients = [(w0, w1, w2)]
errors = [squared_error(sizes, bedrooms, prices, coefficients[-1])]

while True:
    w = coefficients[-1]
    w0i = w0 - alpha*d_error(sizes, bedrooms, prices, w, 0)
    w1i = w1 - alpha*d_error(sizes, bedrooms, prices, w, 1)
    w2i = w2 - alpha*d_error(sizes, bedrooms, prices, w, 2)

    w0 = w0i
    w1 = w1i
    w2 = w2i
    
    coefficients.append((w0, w1, w2))
    errors.append(squared_error(sizes, bedrooms, prices, coefficients[-1]))
    
    if (abs(errors[-2] - errors[-1]) < 0.000001):
        break

print("Mean squared error: {}".format(squared_error(sizes, bedrooms, prices, coefficients[-1])))
print("f(x) = {}x_2 + {}x_1 + {}".format(w2, w1, w0))


# Usando a função encontrada, vamos visualizar o plano que ela representa:

# In[6]:


import numpy as np

plot = go.Scatter3d(x=sizes, y=bedrooms, z=prices, mode='markers',
                   marker=dict(size=12, line=dict(color='rgba(217, 217, 217, 0.14)', width=0.5),
                               opacity=0.8))

x = np.linspace(1, 2, 50)
y = np.linspace(2, 4, 50)

z = []
for i in range(len(x)):
    z.append([])
    for j in range(len(y)):
        z[-1].append(f([x[i], y[j]], coefficients[-1]))

plane = go.Surface(x=x, y=y, z=z, showscale=False)

camera = dict(
    up=dict(x=0, y=0, z=0),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=2, y=-2, z=0.5)
)

layout = go.Layout(width=900, height=400,
                   margin=go.layout.Margin(l=10, r=10, b=10, t=10, pad=10),
                   scene=dict(xaxis=dict(title="tamanho"), yaxis=dict(title="quartos"), zaxis=dict(title="preço")))

data = [plot, plane]

fig = go.Figure(data=data, layout=layout)
fig['layout'].update(scene=dict(camera=camera))

plotly.offline.iplot(fig)