import numpy as np
import matplotlib.pyplot as plt

Regressão Linear

Muitas vezes, usamos equações para representar comportamentos físicos ou da vida. Nesta disciplina, já vimos equações que representam movimentos, dinâmicas de animais, produções de palavras, e assim por diante.

Um problema que acontece com modelos é que, às vezes, temos uma boa pista de que ele tem uma forma geral, mas não conhecemos seus parâmetros.

Num objeto em aceleração constante, por exemplo, podemos saber que sua trajetória obedece:

\[ s = s_0 + v_0t + \frac{at^2}{2}, \]

mas talvez não consigamos saber, ao menos em princípio, quais são os valores de \(s_0\), \(v_0\), e de \(a\).

Quando usamos dados medidos para encontrar parâmetros de um modelo, dizemos que estamos fazendo uma regressão.

Regressão Linear

Até o momento, a função que executamos para o problema de aproximação de dados:

\[ f(x) = \hat{y} = ax + b. \]

Quando fazemos a regressão linear, estamos assumindo que existem várias medições \((x_i, y_i)\). Cada uma dessas medições é diferente do valor estimado, isto é:

\[ e_i = y_i - \hat{y}_i \]

A média do quadrado dessa diferença pode ser calculada ao longo das \(n\) medições que fizermos:

\[ \text{EQM} = \frac{1}{n} \sum_{i=1}^n ||e_i||^2. \]

Isso é uma maneira absolutamente válida de pensar. Porém, vamos agora transportar essa formulação para usar matrizes, e em depois verificaremos como isso pode nos ajudar.

Veja, poderíamos escrever nossa função \(f(x)\) assumindo que a entrada é um vetor de medidas \(x_i\). Daí, ficaríamos com:

\[ \begin{array}{rl} f(\boldsymbol{x}) &= a \boldsymbol{x} + b = a \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} + \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} b = \begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_n \end{bmatrix} \end{array} \]

Nessa formulaçao, a diferença entre as predições e os valores medidos pode ser expressa em um vetor:

\[ e = y - \hat{y} = \begin{bmatrix} y_1 - \hat{y}_1 \\ y_1 - \hat{y}_1 \\ \vdots \\ y_n - \hat{y}_n \end{bmatrix} = \begin{bmatrix} e_1 \\ e_2 \\ \vdots \\ e_n \end{bmatrix} \]
Mostre que o EQM pode ser calculado como: \(\text{EQM} = \frac{1}{n} e^T e\)
\[ \begin{array}{rl} e &= \begin{bmatrix} e_1 \\ e_2 \\ \vdots \\ e_n \end{bmatrix} \\ e^T e &= \begin{bmatrix} e_1 & e_2 & \cdots & e_n \end{bmatrix} \begin{bmatrix} e_1 \\ e_2 \\ \vdots \\ e_n \end{bmatrix} \\ &= e_1^2 + e_2^2 + \cdots + e_n^2 \\ &= ||e||^2 \end{array} \]

EXERCÍCIO

Analise o código-exemplo abaixo.

a) Modifique-o para que ele imprima na tela o EQM da estimação

b) Modifique o valor de a_estimado e b_estimado de forma a reduzir o EQM ao máximo. Qual é o menor EQM que você consegue encontrar?

# Exemplo

# Começamos com a, b "reais" (que tentaremos identificar)
a_real = 2.5
b_real = -3.2

# Vamos medir nosso fenômeno em várias posições x:
n = 100 # Esse é o número de medições
x_medido = np.linspace(-10,10,n).reshape( (n,1) )

# O fenômeno "real" é:
y_real = a_real * x_medido + b_real

# O fenômeno "medido" tem um pouco de erro:
alpha = 4
phi = np.random.normal(size=(n,1))
y_medido = y_real + alpha * phi

# Vamos "chutar" um modelo estimado:
a_estimado = 1.0
b_estimado = -1.0

# Com os parâmetros estimados do modelo, criamos saídas estimadas:
y_estimado = x_medido * a_estimado + b_estimado

plt.figure(figsize=(6,3))
plt.scatter(x_medido, y_medido, c='r', label='Pontos medidos')
plt.plot(x_medido, y_estimado, c='b', label='Modelo estimado')
plt.legend()
plt.show()
No description has been provided for this image

Regressão linear com pseudoinversa

Nesse momento, temos uma formulação como:

\[ \begin{bmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_n \end{bmatrix} = \begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} \]

Vamos nos referir a essas matrizes usando a notação típica de Machine Learning:

\[ \hat{y} = X w^T \]

Vamos analisar rapidamente as matrizes que estamos construindo.

A matriz de estimativas \(\hat{y}=\begin{bmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_n \end{bmatrix}\). Ela é uma matriz \(n \times 1\), isto é, há \(1\) resultado para cada uma das medidas que fizermos.

A matriz de entradas \(X=\begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{bmatrix}\). Ela é uma matriz \(n \times 2\), onde \(n\) é o número de medidas coletadas e temos duas colunas.

Por fim, a matriz \(w^T = \begin{bmatrix} a \\ b \end{bmatrix}\) contém os "pesos" - ou, os parâmetros lineares do modelo - que serão aplicados a cada um dos resultados colocados na matriz \(X\) de forma a gerar as saídas correspondentes.

Nosso problema de estimativa de parâmetros pode ser refraseado da seguinte forma: gostaríamos de estimar os elementos de \(w^T\) de tal forma que a estimativa \(\hat{y}\) se aproxime dos valores medidos \(y\) minimizando o erro calculado por \(e^Te\).

Pseudoinversa

Nosso modelo gera estimativas para \(y\) à partir de medições:

\[ \hat{y} = X w^T \]

Se o modelo fosse perfeito, ele iria gerar exatamente os valores medidos para \(y\):

\[ y = X w^T \]

sabemos que, se \(X\) fosse inversível, encontrar \(w^T\) à partir das medições seria muito simples. Porém, \(X\) não é sequer uma matriz quadrada.

Vamos usar então a seguinte estratégia. Vamos multiplicar os dois lados da equação, à esquerda, por \(X^T\):

\[ X^T y = X^T X w^T \]

A matriz $ X^T X$ é quadrada (\(n \times n\)). Então, podemos multiplicar os dois lados da equação por sua inversa, chegando a:

\[ \begin{array}{rl} X^T y &= (X^T X) w^T \\ (X^T X)^{-1} X^T y &= (X^T X)^{-1} (X^T X) w^T \\ (X^T X)^{-1} X^T y &= w^T \\ \end{array} \]

A matriz \((X^T X)^{-1} X^T\) é chamada de pseudo-inversa de \(X\) e usamos a notação \(X^+\) para nos referir a ela.

Veja como agora podemos encontrar \(w^T\) em função de nossas medições:

\[ w^T = X^+ y \]

EXERCÍCIO

O código abaixo calcula a_estimado e b_estimado usando a pseudo-inversa. Modifique-o para que ele calcule também o EQM da aproximação.

# Exemplo

# Começamos com a, b "reais" (que tentaremos identificar)
a_real = 2.5
b_real = -3.2

# Vamos medir nosso fenômeno em várias posições x:
n = 100 # Esse é o número de medições
x_medido = np.linspace(-10,10,n).reshape( (n,1) )

# O fenômeno "real" é:
y_real = a_real * x_medido + b_real

# O fenômeno "medido" tem um pouco de erro:
alpha = 4
phi = np.random.normal(size=(n,1))
y_medido = y_real + alpha * phi

# Colocando as estimativas na forma matricial:
X = np.hstack ( (x_medido, np.ones_like(x_medido)) )
pinvX = np.linalg.inv( X.T @ X ) @ X.T
wT = pinvX @ y_medido

# Com os parâmetros estimados do modelo, criamos saídas estimadas:
a_estimado = wT[0,0]
b_estimado = wT[1,0]
y_estimado = a_estimado * x_medido + b_estimado

plt.figure(figsize=(6,3))
plt.scatter(x_medido, y_medido, c='r', label='Pontos medidos')
plt.plot(x_medido, y_estimado, c='b', label='Modelo estimado')
plt.title(f'A={a_estimado:.2f}, B={b_estimado:.2f}')
plt.legend()
plt.show()
No description has been provided for this image

Regressão para polinômios

Até o momento, usamos a formulação:

\[ \begin{bmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_n \end{bmatrix} = \begin{bmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} \]

Porém, talvez queiramos usar uma função mais complicada que \(f(x)=ax+b\). Talvez queiramos, por exemplo, modelar alguma função polinomial. Nesse caso, vamos ter que usar matrizes de dados com mais colunas. Se usarmos a forma: \(f(x) = ax^2 + bx + c\), teremos:

\[ \begin{bmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_n \end{bmatrix} = \begin{bmatrix} x^2_1 & x_1 & 1 \\ x^2_2 & x_2 & 1 \\ \vdots & \vdots \\ x^2_n & x_n & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} \]

Veja como nosso modelo nos permite calcular o valor dessa função em \(x_k=3\) usando:

\[ f(3) = \begin{bmatrix} 3^2 & 3 & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} 9 & 3 & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} \]

EXERCÍCIO

O código-exemplo abaixo implementa a estimação de parâmetros de \(y=ax^2 + bx + c\). Modifique-o para que ele passe a estimar os parâmetros da função \(y=a\cos(x) + bx + c\).

# Exemplo

# Começamos com a, b "reais" (que tentaremos identificar)
a_real = 1.1
b_real = 0.2
c_real = 3

# Vamos medir nosso fenômeno em várias posições x:
n = 100 # Esse é o número de medições
x_medido = np.linspace(-10,10,n).reshape( (n,1) )

# O fenômeno "real" é:
y_real = a_real * x_medido**2 + b_real * x_medido + c_real

# O fenômeno "medido" tem um pouco de erro:
alpha = 4
phi = np.random.normal(size=(n,1))
y_medido = y_real + alpha * phi

# Colocando as estimativas na forma matricial:
X = np.hstack ( (x_medido**2, x_medido, np.ones_like(x_medido)) )
pinvX = np.linalg.inv( X.T @ X ) @ X.T
wT = pinvX @ y_medido

# Com os parâmetros estimados do modelo, criamos saídas estimadas:
a_estimado = wT[0,0]
b_estimado = wT[1,0]
c_estimado = wT[2,0]
y_estimado = a_estimado * x_medido**2 + b_estimado * x_medido + c_estimado

plt.figure(figsize=(6,3))
plt.scatter(x_medido, y_medido, c='r', label='Pontos medidos')
plt.plot(x_medido, y_estimado, c='b', label='Modelo estimado')
plt.title(f'A={a_estimado:.2f}, B={b_estimado:.2f}, C={c_estimado:.2f}')
plt.legend()
plt.show()
No description has been provided for this image

Exercício

O código abaixo faz o download de uma base de dados com medidas de anos de experiência (\(x\)) e salário (\(y\), em dólares por ano).

  1. Tente encontrar um modelo linear \(y = ax + b\) para esses dados usando a estratégia da pseudo-inversa. De acordo com o que você encontrou, qual é o salário esperado para alguém acabando de entrar no mercado, e quanto o salário aumenta por ano de experiência? Qual é o EQM encontrado?

  2. Evolua seu modelo para um polinômio de segundo grau. Observando os resultados (os coeficientes, o gráfico e o EQM), como esse modelo se compara com o inicial?

  3. Por fim, vamos tentar um modelo do tipo $y = a e^x + b. Observando os resultados (os coeficientes, o gráfico e o EQM), como esse modelo se compara com o inicial?

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

URL = "https://raw.githubusercontent.com/debakarr/machinelearning/refs/heads/master/Part%202%20-%20Regression/Simple%20Linear%20Regression/Salary_Data.csv"
df = pd.read_csv(URL)
# print(df.head())
x = df['YearsExperience'].to_numpy()
y = df['Salary'].to_numpy()

plt.figure( figsize=(6,3) )
plt.scatter(x, y, c='r')
plt.xlabel('Years of Experience')
plt.ylabel('Salary')
plt.show()
No description has been provided for this image