Regressão Linear na Unha: Construindo um Estimador do Zero

Aula 1 da matéria do Dr. Boldt. O notebook dele quase não tem célula de markdown, é ele codando direto, então boa parte do trabalho deste post é meu: recontar o "por quê" de cada bloco de código, com o Bishop do lado como referência. Bora.

O dataset: 442 pacientes de verdade

Diferente da playlist da especialização do Andrew Ng, aqui o dataset já vem embutido no próprio scikit-learn: o load_diabetes, com 442 pacientes reais, 10 variáveis de entrada (idade, sexo, índice de massa corporal, pressão arterial média e seis medidas de sangue) e um alvo, uma medida numérica de quanto a diabetes progrediu depois de um ano.

from sklearn.datasets import load_diabetes
data = load_diabetes()
X = data.data   # (442, 10)
y = data.target # (442,)

Um detalhe que pega muita gente desprevenida: as 10 colunas de X não estão nas unidades originais. O scikit-learn já entrega esse dataset com cada coluna centralizada em zero e reescalada. Por isso o índice de massa corporal (a coluna que eu vou usar o post inteiro) varia entre -0.09 e 0.17 em vez de ficar na faixa normal de 15 a 40. Guarda essa informação, ela explica por que os coeficientes que vão aparecer mais pra frente parecem números gigantes: eles estão compensando uma variável de entrada minúscula.

O modelo mais simples possível

$y = a \cdot x + b$

A mesma cara de sempre: uma reta, onde a é a inclinação (o quanto o alvo muda pra cada unidade de índice de massa corporal) e b é o intercepto. O professor testa dois chutes na mão, sem otimizar nada ainda, só pra ver a cara da reta:

def modelo_linear(X, a, b):
    return a*X + b

ypred = modelo_linear(X2, 1000, 100)  # primeiro chute
ypred = modelo_linear(X2, 1500, 150)  # segundo chute

Mexe nos sliders abaixo (aqui no blog eu uso a notação padrão w/b, é exatamente o mesmo a/b do professor) numa amostra real de 45 pacientes do conjunto de treino e vê você mesmo qual reta cai melhor em cima dos pontos:

Como saber se um chute é bom

"Cai melhor" precisa virar número. O professor implementa três métricas na unha, a mesma trinca que já apareceu na playlist da especialização (se você já leu aquela, pode passar os olhos rápido):

• MAE (erro absoluto médio): a média de |previsto - real|. Fácil de interpretar (tá em unidade de paciente), mas trata todo erro do mesmo jeito.
• MSE (erro quadrático médio): a média de (previsto - real)². Eleva ao quadrado, então um erro grande pesa desproporcionalmente mais que vários erros pequenos.
• RMSE: a raiz do MSE, pra voltar pra unidade original depois de ter quadrado tudo.

def mae(y, ypred):
    return np.sum(np.abs(y - ypred)) / len(y)

def mse(y, ypred):
    return np.sum((y - ypred)**2) / len(y)

def rmse(y, ypred):
    return np.sqrt(mse(y, ypred))

Saída no dataset inteiro (442 pacientes): chute (1000, 100) → MAE 64.93, MSE 6614.14, RMSE 81.33. Chute (1500, 150) → MAE 53.93, MSE 4580.80, RMSE 67.68. O segundo chute é melhor nas três métricas.

O Bishop chama isso de função de erro, e define ela quase igual ao MSE, só que sem dividir pela quantidade de pontos e com um fator 1/2 de bônus:

$E(\mathbf{w}) = \frac{1}{2}\sum_{n=1}^{N}\{y(x_n,\mathbf{w}) - t_n\}^2$

O 1/2 existe só por conveniência: quando você deriva essa função pra achar o mínimo, o quadrado desce multiplicando e cancela o 1/2 exatamente. Detalhe estético, não muda pra onde o mínimo aponta. E a versão "raiz" que ele usa pra comparar datasets de tamanhos diferentes,

$E_{\text{RMS}} = \sqrt{2E(\mathbf{w}^*)/N}$

é a mesma RMSE que a gente acabou de calcular, com a mesma divisão por N escondida atrás do fator 2 (porque a função dele não divide por N, só a nossa MSE já divide).

O modelo idiota que vira régua

Antes de qualquer regressão de verdade, o segundo notebook (aula02a) faz uma coisa esperta: separa 20% dos dados pra teste (train_test_split implementado na mão, random_state=42 pra dar sempre o mesmo resultado) e define o modelo mais preguiçoso que existe, chutar sempre a média:

def modelo_media(y):
    media = np.mean(y)
    return np.ones_like(y) * media

Saída (353 pacientes de treino): MSE do chute (1500, 150) = 4646.20. MSE do "sempre a média" = 5978.59.

Por que eu preciso disso? Porque um MSE de 4646 sozinho não diz nada, é só um número. Comparado com o baseline da média, ele fica menor, então o modelo tá aprendendo alguma coisa. Isso é literalmente a definição de R²:

def r2_score(y, ypred):
    ss_res = mse(y, ypred)
    ss_tot = mse(y, modelo_media(y))
    return 1 - (ss_res / ss_tot)

$R^2 = 1 - \frac{\text{MSE do modelo}}{\text{MSE do baseline}}$

Com o chute (1500, 150): R² = 1 - 4646.20/5978.59 = 0.2229. Só uma variável explicando 22% da variação já é um começo, considerando que a diabetes depende de um monte de coisa que não tá nem nessas 10 colunas. E o R² do próprio baseline, por definição, é sempre zero.

A API do scikit-learn, construída por dentro

Aqui que a aula fica interessante de verdade. Em vez de só chamar LinearRegression() pronto, o professor ensina a construir um estimador que segue a mesma convenção do scikit-learn, herdando de duas classes-base:

from sklearn.base import BaseEstimator, RegressorMixin

class AverageRegressor(BaseEstimator, RegressorMixin):
    def fit(self, X, y):
        self.media_ = np.mean(y)
        return self
    def predict(self, X):
        return np.ones(shape=(X.shape[0],)) * self.media_

BaseEstimator dá de graça um monte de infraestrutura (comparação de parâmetros, clonagem, integração com Pipeline e GridSearchCV, coisas que vão aparecer só lá na frente da matéria). RegressorMixin adiciona o método .score() padrão. E repara na convenção do underscore: self.media_, terminando em _, marca "isso foi aprendido dos dados", diferente de um parâmetro configurado na mão. É o mesmo detalhe que eu já tinha visto de fora, usando SGDRegressor na outra playlist, só que agora eu tô escrevendo o _ com a própria mão.

O próximo passo mostra algo sutil: um regressor com coeficientes fixos, passados no construtor, sem aprender nada:

class LinearRegressor(BaseEstimator, RegressorMixin):
    def __init__(self, a, b):
        self.a_ = a
        self.b_ = b
    def fit(self, X, y):
        return self  # não aprende nada, os coeficientes já vieram prontos
    def predict(self, X):
        return (self.a_*X + self.b_).reshape(X.shape[0],)

E depois, a versão que inicializa aleatório:

class LinearRegressor(BaseEstimator, RegressorMixin):
    def fit(self, X, y):
        self.a_ = np.random.rand()
        self.b_ = np.random.rand()
        return self
    # ...

Saída: MSE do chute (1500, 150) do professor = 4646.20. MSE do coeficiente aleatório = 29993.30.

Isso rende um R² de -4.02 para o chute aleatório. Negativo: pior que simplesmente chutar a média toda vez. É a prova de que "aprender" não é mágica, é sair de um lugar ruim (aleatório) na direção de um lugar melhor, e um número aleatório sozinho não tem nenhum motivo pra ser bom.

Do chute aleatório ao gradiente

A última versão da classe substitui "aleatório e parado" por "aleatório e caminhando":

class LinearRegressor(BaseEstimator, RegressorMixin):
    def __init__(self, max_iter=1000, learning_rate=0.001):
        self.max_iter = max_iter
        self.learning_rate = learning_rate

    def fit(self, X, y):
        self.coefs_ = np.random.rand(X.shape[1])
        self.intercept_ = np.random.rand()
        for i in range(self.max_iter):
            y_pred = self.predict(X)
            error = y - y_pred
            self.coefs_ += X.T @ error * self.learning_rate
            self.intercept_ += error.sum() * self.learning_rate
        return self

    def predict(self, X):
        return (X @ self.coefs_ + self.intercept_).reshape(X.shape[0],)

Se você já passou pela playlist da especialização, essa conta te é familiar: é o mesmo gradiente descendente batch de sempre, coeficiente += X.T @ erro * taxa, só que sem dividir pela quantidade de exemplos (o que só empurra o efeito pra dentro da learning_rate, um pouco menor compensa a falta da divisão). Esse formato específico, atualizar o peso proporcionalmente ao erro vezes a entrada, tem nome próprio na literatura: regra delta (ou regra de Widrow-Hoff), um dos algoritmos de aprendizado mais antigos que existem, décadas antes do termo "gradiente descendente" pegar.

Roda ao vivo abaixo, na mesma amostra de 45 pacientes, e repara como ele sai de longe e converge pra perto do ajuste ideal:

Nessa amostra de 45 pontos (menor que os 353 do treino inteiro, então o ajuste ideal aqui não é exatamente igual ao número que a aula reporta), com alpha=0.5 clicando em "Rodar 2000" ele já pousa bem perto do fundo do vale. Testa também um alpha maior, tipo 1.5: o b começa a oscilar antes de convergir, porque a escala de a (multiplicado por um x minúsculo, entre -0.09 e 0.17) e a escala de b (que soma direto) são bem diferentes, o mesmo problema de escala que a normalização resolve, só que aqui apareceu escondido dentro do próprio dataset já normalizado.

No dataset de treino inteiro (353 pacientes), rodando essa mesma classe com uma única feature (o índice de massa corporal):

Saída: coefs_ = [529.65], intercept_ = 154.82, MSE = 4302.16.

Isso já bate o melhor chute manual (4646.20 → 4302.16) e rende R² = 1 - 4302.16/5978.59 = 0.2804. Sem eu escolher a e b na mão, o gradiente achou um ajuste melhor sozinho.

De uma variável pra dez

O dataset tem 10 colunas, e até agora eu só usei uma (o índice de massa corporal). A mesma classe, sem nenhuma mudança de código, aceita X_train inteiro:

regressor = LinearRegressor()
regressor.fit(X_train, y_train)  # agora com as 10 colunas

Saída: MSE = 3027.12 (contra 4302.16 usando só uma variável). RMSE cai de 65.59 para 55.02. R² sobe de 0.2804 para 0.4937.

Quase o dobro de variância explicada, só por deixar o modelo enxergar as outras 9 variáveis (idade, sexo, pressão, as seis medidas de sangue) que ele simplesmente não tinha acesso antes. Ninguém mudou o algoritmo, só a quantidade de informação que ele recebeu.

Fechando

Duas coisas ficaram em aberto de propósito. Primeiro: resolver isso tudo iterando, passo a passo, funciona, mas regressão linear tem solução exata, fechada, numa única conta, e é isso que a próxima aula (aula02b, equação normal) resolve. Segundo: eu não toquei nos outros modelos que a aula compara (aula02c), fica pro próximo post também.

Aplicação Prática

Uso o mesmo dataset de diabetes que já apareceu o post inteiro, sem introduzir nada de novo: o objetivo aqui é só visualizar o quão bom ficou o ajuste de 10 variáveis que acabei de mostrar.

regressor = LinearRegressor(max_iter=1000, learning_rate=0.001)
regressor.fit(X_train, y_train)
y_pred = regressor.predict(X_train)

Reproduzi essa mesma classe (com uma semente aleatória diferente, então os coeficientes batem só na segunda casa decimal com os do professor, o MSE final ficou em 3027.43, praticamente idêntico ao 3027.12 reportado na aula) pra extrair os pares de valor real contra valor previsto de 40 pacientes do treino. Cada ponto é um paciente: quanto mais perto da linha pontilhada (a previsão perfeita, previsto = real), melhor o modelo acertou aquele caso.

Dá pra ver a nuvem de pontos seguindo a diagonal, mas com bastante espalhamento, exatamente o que um R² de 0.49 significa na prática: o modelo capturou boa parte do padrão, mas ainda erra bastante caso a caso. Nada mau pra um modelo linear simples, sem nenhuma feature nova, sem regularização, só 10 números crus e 1000 passos de gradiente.