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

# # Regressão Logística

# ## Introdução
# 
# Na aprendizagem de máquina supervisionada, quando a resposta que precisa ser prevista é qualitativa (discreta), é usado o processo chamado de Classificação. São problemas de classificação, por exemplo, decidir se um e-mail recebido é spam ou não, ou decidir se uma foto contém a imagem de um gato, cachorro ou um pássaro.
# 
# Para funcionar, um classificador precisa ser treinado usando dados de treinamento. Cada exemplo nos dados de treinamento possui um vetor de características $(x_1, x_2, ..., x_n)$ e uma resposta associada $y$ que é a classificação do exemplo.
# 
# Apesar do nome, Regressão Logística é um algoritmo de classificação que pode ser usado para resolver problemas com duas ou mais possíveis classes. Entretanto, para facilitar o entendimento, consideraremos que $y$ pode ter apenas dois valores, 0 ou 1, sendo 0 uma classe e 1 a outra.

# ## Função logística
# 
# Com as respostas variando apenas entre dois valores, 0 e 1, usar regressão linear não é adequado, pois ela pode prever qualquer valor contínuo que pode ir muito além dos limites de 0 e 1. Assim, mudaremos a função $f(x)$, que representa a hipótese para explicar os dados, para uma função que varie apenas entre 0 e 1.
# 
# Na regressão logística usaremos a função logística da seguinte forma:
# 
# $$f(x) = g(w^Tx)$$
# 
# $$g(z) = \frac{1}{1+e^{-z}}$$
# 
# Perceba que continua-se usando $w^Tx$ como na regressão linear, mas agora aplica-se o resultado na função logística $g(z)$.
# 
# A seguir temos a implementação das duas funções acima, $f(x)$ e $g(z)$. $f(x)$ recebe o vetor de coeficientes $w$ e as características $x$, já $g(z)$ recebe apenas um valor que é aplicado a função logística.

# In[1]:


import math

def f(x, w):
    z = 0

    for i in range(len(x)):
        z += x[i]*w[i]
    
    return g(z)

def g(z):
    return 1 / (1 + math.exp(-z))


# A função logística tem uma forma de "S" como pode ser vista no gráfico seguir:

# In[2]:


import numpy as np
from bokeh.io import show, output_notebook
from bokeh.plotting import figure
output_notebook()

x_points = np.linspace(-10, 10, 100)
y_points = [g(i) for i in x_points]

p = figure(plot_width=700, plot_height=250)
p.xaxis.axis_label = "z"
p.yaxis.axis_label = "g(z)"
p.line(x_points, y_points, line_width=2)
show(p, notebook_handle=True)


# Note que quanto menor o valor de $z$, ou seja, o valor de $w^Tx$, mais próximo de zero é o resultado da função. Por outro lado, quanto maior o valor de $z$, mais próximo de 1 é o resultado. Dessa forma, o resultado nunca ultrapassará o intervalo entre 0 e 1. 

# Na regressão logística, interpretamos o resultado de $f(x)$ como a probabilidade da resposta ser 1 para o vetor de características $x$. Assim, se $f(x) = 0,8$, então $x$ tem probabilidade $0,8$ de ser da classe 1, e $0,2$ de ser da classe 0.

# ## Fronteira de decisão
# 
# Para obter uma classificação discreta, é necessário definir um limiar em relação à probabilidade. Por exemplo, se $f(x) \geq 0,5 \Rightarrow y = 1$. Ou seja, se a probabilidade de ser 1 for maior ou igual 50%, então classificamos como 1. Adicionalmente, se $f(x) < 0,5 \Rightarrow y = 0$.
# 
# Na função logística $g(z)$, se $z \geq 0$, então $g(z) \geq 0,5$. Como a hipótese na regressão logística é $f(x) = g(w^Tx)$, então, se $w^Tx \geq 0$, isto resultará em $g(z) \geq 0,5$.
# 
# Com isso, $w^Tx$ especifica uma fronteira de decisão que separa os pontos que são classificados como 0 ou 1. Definimos matematicamente a fronteira de decisão da seguinte forma:
# 
# $$w^Tx = 0 $$
# 
# Assim,
# 
# $$w^Tx \geq 0 \Rightarrow y = 1$$
# 
# $$w^Tx < 0 \Rightarrow y = 0$$
# 
# Dependendo do tipo de função definida por $w^Tx$ a fronteira de decisão pode ter diferentes formas. 
# 
# 
# Por exemplo, usando a função $f(x) = w_0 + w_1x_1 + w_2x_2$, definimos a fronteira como:
# 
# $$w_0 + w_1x_1 + w_2x_2 = 0$$
# 
# Essa equação define uma reta.
# 
# 
# Abaixo temos o gráfico mostrando a fronteira usando a equação da reta anterior para $w^T = \begin{bmatrix} -3 & 0,8 & 1 \end{bmatrix}$.
# Ou seja, a equação é: $0,8x_1 + x_2 - 3 = 0$.
# 
# Note que no gráfico, os pontos são separados em dois grupos pela reta, sendo os de cor azul correspondentes a $y=0$ e os de cor vermelha a $y=1$.

# In[3]:


def boundary(x):
    return 3-0.8*x

x = [0.5, 1, 1.2, 1.1, 1.5, 2.6, 2.3, 3, 2.8, 2.4]
y = [0.7, 1.2, 0.5, 1.8, 1, 2.5, 2.8, 2.5, 3, 2.2]
colors = ["blue", "blue", "blue", "blue", "blue", "red", "red", "red", "red", "red"]

p = figure(plot_width=700, plot_height=250)
p.circle(x, y, radius=0.025, color=colors)
p.line([0, 3.5], [boundary(0), boundary(3.5)], line_width=2)
show(p, notebook_handle=True)


# Usando funções diferentes, podemos ter curvas mais complexas. 
# No exemplo a seguir usamos a função $f(x) = -1,2 + 3x_1^2 + 3x_2^2$, criando a fronteira $3x_1^2 + 3x_2^2 - 1,2 = 0$. Esta equação define uma elipse.

# In[4]:


import math
import itertools

def boundary(x1):
    return math.sqrt((1.2-3*x1**2)/0.3)

x = [0.5, -0.1, 0.2, -0.3, -0.3, 0, 0.7, 1, 0.8, 0, -0.5, -0.8, -0.5]
y = [0.2, -0.8, 0.5, -0.1, 0.7, 3.5, 2.0, 0.5, -1.5, -3.5, -2, 0, 2]
colors = ["blue", "blue", "blue", "blue", "blue", "red", "red", "red", "red", "red", "red", "red", "red"]

curve_x_points = np.linspace(-0.8, 0.8, 100000)
curve_x_points = [i for i in curve_x_points if ((1.2-3*i**2)/0.3) >= 0]

curve_y_points = [boundary(i) for i in curve_x_points]
curve_neg_y_points = [-boundary(i) for i in curve_x_points]

p = figure(plot_width=700, plot_height=250)
p.circle(x, y, radius=0.015, color=colors)

p.line(curve_x_points, curve_y_points, line_width=2)
p.line(curve_x_points, curve_neg_y_points, line_width=2)
show(p, notebook_handle=True)


# Dois grupos foram criados pela fronteira. Os pontos que estão fora da elipse possuem classificação $y=1$ e os que estão dentro possuem classificação $y=0$.

# ## Função de erro
# 
# Para descobrir os valores dos coeficientes em $w$ vamos, assim como nas outras regressões, minimizar uma função de erro através da descida de gradiente.
# 
# A função de erro é definida assim:
# 
# $$J(w) = \frac{1}{m} \sum_{i=1}^{m} Cost(f(x^{(i)}), y^{(i)})$$
# 
# 
# Onde $m$ é a quantidade de exemplos de treinamento.
# 
# $Cost(f(x), y)$ varia de acordo com o valor de $y$. Para $y=1$, temos:
# 
# $$Cost(f(x), y) = -\log(f(x))$$
# 
# Assim, quanto mais próximo de 1 for $f(x)$, mais próximo de 0 é o valor de $Cost$. Caso contrário, quanto mais próximo de 0 for $f(x)$, mais próximo de infinito é o valor de $Cost$.
# 
# Para $y=0$:
# 
# $$Cost(f(x), y) = -\log(1 - f(x))$$
# 
# Nesse caso, temos o oposto, quanto mais próximo de 0 for $f(x)$, mais próximo de 0 é o valor de $Cost$ e quanto mais próximo de 1 for $f(x)$, mais próximo de infinito é o valor de $Cost$. 
# 
# Com isso, se o modelo atribui probabilidades altas de forma correta paras as classificações dos exemplos, então o valor de $J(w)$ é baixo. Caso contrário, se ele atribui de forma errada probabilidades altas, o valor de $J(w)$ tende a infinito.
# 
# Os gráficos a seguir mostram como os valores de $Cost$ variam para $y=1$ (esquerda) e $y=0$ (direita).

# In[5]:


from bokeh.layouts import row

x_points = np.linspace(0.01, 0.99, 10000)
y1_points = [-math.log(i) for i in x_points]
y0_points = [-math.log(1-i) for i in x_points]

p1 = figure(plot_width=350, plot_height=250)
p1.xaxis.axis_label = "f(x)"
p1.yaxis.axis_label = "Cost(f(x), 1)"
p1.line(x_points, y1_points, line_width=2)

p0 = figure(plot_width=350, plot_height=250)
p0.xaxis.axis_label = "f(x)"
p0.yaxis.axis_label = "Cost(f(x), 0)"
p0.line(x_points, y0_points, line_width=2)

show(row(p1, p0))


# A função $Cost$ pode ser condensada em apenas uma definição:
# 
# $$Cost(f(x), y) = -y\log(f(x)) - (1-y)\log(1-f(x))$$
# 
# Note que quando $y=1$, a função se resume a $Cost(f(x), y) = -\log(f(x))$ e quando $y=0$ ela tem a forma $Cost(f(x), y) = -\log(1-f(x))$. Ou seja, as mesmas definições mostradas anteriormente.
# 
# Com isso, usaremos essa definição diretamente na função de erro $J(w)$:
# 
# $$J(w) = \frac{1}{m} \sum_{i=1}^{m} [-y^{(i)}\log(f(x^{(i)})) - (1-y^{(i)})\log(1-f(x^{(i)}))]$$
# 
# Para simplificar, multiplicamos todo o somatório por $-1$, assim trocando os sinais internos do somatório.
# 
# $$J(w) = -\frac{1}{m} \sum_{i=1}^{m} [y^{(i)}\log(f(x^{(i)})) + (1-y^{(i)})\log(1-f(x^{(i)}))]$$

# A seguir é mostrada a implementação da função de erro $J(w)$. Ela recebe as características dos dados de treinamento (parâmetro x) e as suas respostas (parâmetro y), além dos coeficientes (parâmetro w). 
# O parâmetro x é uma matriz e cada linha dessa matriz corresponde ao vetor de características de um exemplo.
# Outro ponto importante é lembrar que o primeiro elemento do vetor de características é sempre 1, por isso é feita a concatenação com `[1]` na primeira linha dentro do `for`.

# In[6]:


def cost(x, y, w):
    costs = 0

    for i in range(len(x)):
        xi = np.concatenate([[1], x[i]])
        yi = y[i]
        
        f_value = f(xi, w)

        if yi == 1:
            if (f_value == 0):
                return math.inf
            else:
                costs -= math.log(f_value)
        else:
            if (f_value == 1):
                return math.inf
            else:
                costs -= math.log(1 - f_value)

    return (1/(2*len(x)))*costs


# ## Descida de gradiente
# 
# Vamos aplicar a regressão logística no conhecido conjunto de dados sobre as flores iris. No código abaixo, carregamos e visualizamos dois tipos de flores (azul e vermelho) de acordo com duas características ($x_1$ e $x_2$).

# In[7]:


from sklearn import datasets

iris_data = datasets.load_iris()

x_chart = [i[0] for i in iris_data.data[:100]]
y_chart = [i[1] for i in iris_data.data[:100]]
colors = (["blue"] * 50) + (["red"] * 50)

p = figure(plot_width=700, plot_height=250)
p.xaxis.axis_label = "x1"
p.yaxis.axis_label = "x2"
p.circle(x_chart, y_chart, radius=0.015, color=colors)
show(p, notebook_handle=True)


# Para minimizar $J(w)$ usando a descida de gradiente é necessário calcular a derivada parcial da função:
# 
# \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*}
# 
# As derivadas são calculadas pelo código abaixo. A função recebe as características dos dados de treinamento (parâmetro x), as suas respostas (parâmetro y), os coeficientes (parâmetro w) e o índice $j$ (parâmetro j) da equação da derivada parcial que especifica qual derivada parcial está sendo calculada.

# In[8]:


def d_error(x, y, w, j):
    total = 0

    for i in range(len(x)):
        xi = np.concatenate([[1], x[i]])
        yi = y[i]

        fxi = f(xi, w)
        error = (fxi - yi)

        total += error*xi[j]
    
    return total/(len(x))


# Com isso, atualizamos os coeficientes iterativamente 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$, sendo $N+1$ o total de coeficientes.
# 
# Vamos ajustar uma reta como fronteira de decisão, ela tem a forma: $w_0 + w_1x_1 + w_2x_2 = 0$.

# In[9]:


x = [[i[0], i[1]] for i in iris_data.data[:100]]
y = iris_data.target[:100]

w = [0, 0, 0]
alpha = 0.1

for n in range(10000):
    w_temp = [0, 0, 0]
    for i in range(len(w)):
        w_temp[i] = w[i] - alpha*d_error(x, y, w, i)
        
    w = w_temp


# Os coeficientes encontrados foram:

# In[10]:


print(w)


# Dessa forma, foi encontrada a reta: $-3,3872 + 6,055x_1 -9,5045x_2 = 0$.
# Traçando ela no gráfico, percebemos claramente a divisão entre os dois tipos de flores.

# In[11]:


p = figure(plot_width=700, plot_height=250)
p.xaxis.axis_label = "x1"
p.yaxis.axis_label = "x2"
p.circle(x_chart, y_chart, radius=0.015, color=colors)
p.line([4, 8], [(-3.3872+6.055*4)/9.5045, (-3.3872+6.055*8)/9.5045], line_width=2)
show(p, notebook_handle=True)


# Vamos testar o classificador para os pontos (6, 2) e (4,3).

# In[12]:


print("Probabilidade do ponto (6,2) ter a classificação 1 (vermelho): {}".format(f([1, 6, 2], w)))
print("Probabilidade do ponto (4,3) ter a classificação 1 (vermelho): {}".format(f([1, 4, 3], w)))


# Podemos também usar uma fronteira definida por um polinômio. Por exemplo: $w_0 + w_1x_1 + w_2x_2^2 = 0$.
# Para isso, escreveremos a função `f` e a função para calcular a sua derivada de forma específica para esse exemplo, pois a implementação anterior não tratava expoentes nas variáveis.

# In[13]:


def f(x, w):
    z = w[0]*x[0] + w[1]*x[1] + w[2]*(x[2]**2) 
    return g(z)

def d_error(x, y, w, j):
    total = 0

    for i in range(len(x)):
        xi = np.concatenate([[1], x[i]])
        yi = y[i]

        fxi = f(xi, w)
        error = (fxi - yi)

        if (j == 2):
            total += error*xi[j]**2
        else:
            total += error*xi[j]
    
    return total/(len(x))


# Executando a descida de gradiente:

# In[14]:


x = [[i[0], i[1]] for i in iris_data.data[:100]]
y = iris_data.target[:100]

w = [0, 0, 0]
alpha = 0.01

for n in range(1000000):
    w_temp = [0, 0, 0]
    for i in range(len(w)):
        w_temp[i] = w[i] - alpha*d_error(x, y, w, i)
        
    w = w_temp


# Os coeficientes encontrados foram:

# In[15]:


print(w)


# Traçando a curva no gráfico, temos:

# In[16]:


x_points = np.linspace(4, 8, 100)
y_points = [math.sqrt(((w[0] + w[1]*i))/(-w[2])) for i in x_points]

p = figure(plot_width=700, plot_height=250)
p.xaxis.axis_label = "x1"
p.yaxis.axis_label = "x2"
p.circle(x_chart, y_chart, radius=0.015, color=colors)
p.line(x_points, y_points, line_width=2)
show(p, notebook_handle=True)


# Vamos testar o classificador para os pontos (6, 2) e (4,3).

# In[17]:


print("Probabilidade do ponto (6,2) ter a classificação 1 (vermelho): {}".format(f([1, 6, 2], w)))
print("Probabilidade do ponto (4,3) ter a classificação 1 (vermelho): {}".format(f([1, 4, 3], w)))