Regressão Linear
A regressão linear é o algoritmo mais antigo deste curso. Legendre publicou os mínimos quadrados em 1805 para ajustar órbitas de cometas, num apêndice intitulado Sur la Méthode des moindres quarrés e Gauss alegou usá-los desde 1795 — mas nenhum dos dois chamava aquilo de regressão. Essa palavra chegou oitenta anos depois, vinda de um estudo sobre altura humana e descrevia um fenômeno, não um método. Dois séculos adiante, o modelo continua sendo a primeira coisa a tentar em qualquer alvo contínuo: rápido, interpretável e uma baseline constrangedoramente difícil de superar.
Esta aula pega uma ideia e a segue até o fim: ajuste uma reta, olhe o que ela errou e use esses erros para melhorar. Os erros são os resíduos; usá-los de forma sistemática é o gradiente.
De onde vem o nome
Em 1886 Francis Galton publicou Regression towards Mediocrity in Hereditary Stature. Ele havia coletado as alturas de 928 filhos adultos e de seus pais; procurava medir a força da hereditariedade. O que encontrou o intrigou.
Pais muito altos de fato tinham filhos altos — mas, em média, os filhos ficavam mais perto da média populacional do que seus pais. Pais muito baixos tinham filhos baixos que, de novo, ficavam mais perto do meio. Galton chamou isso primeiro de "reversão", depois de "regressão à mediocridade". O nome grudou na técnica que ele usou para medir o efeito — e é por isso que um algoritmo sobre ajustar retas carrega uma palavra que significa "andar para trás".
O painel da esquerda é o efeito em si. A linha tracejada é a identidade — como os dados seriam se os filhos simplesmente repetissem os pais. A reta ajustada é mais achatada. Um valor médio dos pais 3 polegadas acima da média produz um filho previsto apenas cerca de 2 polegadas acima. A seta vermelha é a diferença e essa diferença é a "regressão".
A inclinação que Galton relatou era de cerca de \(2/3\) e vale ver de onde ela sai:
Correlação e inclinação não são o mesmo número. A inclinação carrega as unidades; a correlação é o que sobra depois de removê-las. As duas coincidem apenas quando ambas as variáveis têm a mesma dispersão — que é exatamente a situação em que o simulador abaixo coloca você, para observar a geometria sem a contabilidade das unidades.
Um número dessa conta merece um segundo olhar: por que os filhos (2,52) variam mais que a média dos pais (1,80)? Não porque os pais sejam um grupo mais uniforme. A média dos pais é a média de duas pessoas e tirar média encolhe a dispersão: dois indivíduos com \(s \approx 2{,}5\) resultam em \(2{,}5/\sqrt{2} \approx 1{,}77\), que é essencialmente o 1,80 que Galton relatou. A assimetria é um artefato de como a variável \(x\) foi construída, não um fato sobre hereditariedade — e perceber isso é o mesmo hábito mental de que trata o resto desta seção.
Regressão à média não é uma força
O painel da direita da figura é a parte que quase ninguém vê. Se você regredir pais sobre filhos, essa reta também é mais achatada — um filho alto prevê pais mais próximos da média. As duas direções "regridem", o que é impossível para uma atração causal: filhos não podem estar puxando os pais em direção à média.
O que acontece de fato é aritmética. Toda medida é sinal mais ruído. Uma observação extrema é extrema em parte porque o sinal era alto e em parte porque o ruído ajudou. O sinal persiste na segunda medida; a sorte não. Então a segunda medida cai mais perto do meio — qualquer que seja a variável que você chame de "primeira".
Essa é uma das formas mais confiáveis de se enganar com dados:
- os alunos com as piores notas recebem reforço e melhoram — então o reforço funciona? Eles teriam melhorado de qualquer jeito;
- o atleta que aparece na capa da revista tem uma temporada seguinte pior — o Sports Illustrated jinx, que é apenas uma temporada excepcional seguida de uma temporada normal;
- uma clínica trata os pacientes com pressão mais alta e vê a pressão cair. Parte da queda é o tratamento. Parte é regressão à média e só um grupo de controle separa as duas.
O modelo
Prever um alvo contínuo como uma soma ponderada de atributos:
- \(w_j\) é a variação em \(\hat{y}\) por unidade de variação em \(x_j\), mantendo os outros atributos fixos;
- \(w_0\) (intercepto, ou viés) é a previsão quando todos os atributos são zero.
Em forma matricial, com uma coluna inicial de 1s absorvida em \(X\): \(\hat{y} = Xw\).
Por que não inverter \(X\)?
Escrito como \(\hat{y} = Xw\), o reflexo natural é \(w = X^{-1}y\). Não funciona — e ver exatamente por que é o caminho mais rápido para entender o que os mínimos quadrados de fato fazem.
Tome três observações e dois parâmetros:
Primeiro obstáculo: \(X\) não é quadrada. Ela é \(3 \times 2\), enquanto apenas matriz quadrada tem inversa. No caso geral \(X\) é \(n \times d\) com \(n \gg d\) — milhares de linhas, um punhado de colunas.
Segundo obstáculo, o verdadeiro: não existe solução exata de qualquer forma. \(y = Xw\) são três equações com duas incógnitas — um sistema sobredeterminado. Um \(w\) exato existe apenas se \(y\) por acaso estiver no espaço-coluna de \(X\), o que aqui não acontece:
Acrescentar \(y\) elevou o posto, então \(y\) aponta para algum lugar que as colunas de \(X\) não alcançam. Nenhum \(w\) satisfaz \(Xw = y\). Isso não é azar — é a situação normal — e a razão inteira de a regressão existir. Se houvesse solução exata, todo resíduo seria zero e a reta passaria por todos os pontos.
O que fazemos no lugar
Como \(y\) está fora de alcance, contentamo-nos com o ponto alcançável mais próximo dele — minimizar \(\lVert y - Xw\rVert^2\). Igualar o gradiente a zero dá as equações normais; o truque está no que acontece com os formatos:
\(X^\top X\) é \(2 \times 2\) — quadrada — e invertível sempre que \(X\) tem posto-coluna cheio. Nunca invertemos \(X\); invertemos \(X^\top X\). São duas equações, pequenas o bastante para resolver no papel:
Os valores ajustados são \(1{,}5,\ 4,\ 6{,}5\) e os resíduos \(0{,}5,\ -1,\ 0{,}5\). Eles não se anulam — não podem — mas \(\lVert e \rVert^2 = 1{,}5\) é o menor alcançável. E repare no que as equações normais dizem nessa notação: \(X^\top e = 0\), que é exatamente o \(\sum e_i = 0\) mais o \(\sum e_i x_i = 0\) que você reencontrará adiante.
O seu instinto é o caso particular
Agora descarte a terceira observação, ficando com duas equações e duas incógnitas:
Aqui \(X^{-1}\) existe e \(w = X^{-1}y\) funciona perfeitamente — o resíduo é exatamente zero, porque uma reta por dois pontos os interpola. E a fórmula de mínimos quadrados devolve exatamente a mesma resposta, como a álgebra exige:
Então \(w = X^{-1}y\) não está errado — é o caso exatamente determinado de uma fórmula mais geral. Essa fórmula geral, \(X^{+} = (X^\top X)^{-1}X^\top\), chama-se pseudo-inversa de Moore–Penrose: o que a "inversa" vira quando a matriz não é quadrada.
Ninguém calcula essa inversa
A fórmula é como a solução se escreve, não como ela é calculada. Formar \(X^\top X\) eleva o número de condição ao quadrado, então as bibliotecas resolvem o sistema direto por decomposição QR ou SVD — mais rápido e muito mais estável. O LinearRegression do scikit-learn chama scipy.linalg.lstsq, que nunca constrói inversa alguma.
Aquelas quatro palavras — mantendo os outros atributos fixos — pesam mais do que parecem. Leia este modelo ajustado para preços de apartamentos:
3200 não significa "apartamentos maiores custam R$ 3200 a mais por m² nesta cidade". Significa: entre apartamentos de mesma idade, um m² a mais está associado a R$ 3200 a mais. Se área e idade forem correlacionadas nos seus dados — digamos, prédios mais novos são maiores — então a relação bruta entre área e preço mistura os dois efeitos — e o coeficiente os separa deliberadamente.
Os coeficientes não são importâncias
Um reflexo comum é ordenar atributos por \(|w_j|\). Isso não significa nada, a menos que os atributos compartilhem uma escala: meça a área em km² em vez de m² e o coeficiente cresce um milhão de vezes sem que nada no modelo mude. Para comparar magnitudes, padronize os atributos antes — e mesmo assim "coeficiente grande" quer dizer "resposta íngreme", não "importante".
Mínimos quadrados
Precisamos de uma regra para escolher \(w\). Os mínimos quadrados ordinários — OLS, de ordinary least squares, como aparecerá daqui em diante — escolhem a reta que minimiza a soma dos quadrados dos resíduos (SSE, de sum of squared errors): as distâncias verticais ao quadrado entre os dados e a reta ajustada.
Por que distâncias verticais e por que ao quadrado? As duas escolhas merecem um instante, porque as duas poderiam ter sido diferentes.
Verticais, porque o trabalho do modelo é prever \(y\) a partir de \(x\). Um erro é "o quanto minha previsão de \(y\) errou", medido ao longo do eixo \(y\). (Minimizar a distância perpendicular é um método diferente e perfeitamente válido — é o que a PCA faz — mas responde a outra pergunta e é por isso que as duas retas de regressão da figura de Galton diferem.)
Ao quadrado, por três razões que por acaso se alinham:
- é suave e diferenciável em toda parte, então o cálculo funciona e existe resposta em forma fechada;
- pune um erro grande mais do que vários pequenos, que costuma ser o que se quer;
- sob ruído gaussiano, é máxima verossimilhança — a reta de mínimos quadrados é a reta mais provável.
A razão 2 é também a sua fraqueza: elevar ao quadrado deixa o ajuste hipersensível a outliers. Minimizar \(\sum|e_i|\) dá os desvios absolutos mínimos, mais robustos, ao preço de perder a forma fechada.
Arraste os pontos abaixo e veja a reta perseguir o mínimo. Depois arraste um ponto para longe dos demais e veja o quanto uma única observação consegue movê-la:
Resolvendo exatamente
Iguale o gradiente a zero, \(\nabla_w J = -2X^\top(y - Xw) = 0\) e você obtém as equações normais:
Para a regressão simples — um atributo — isso se reduz a duas fórmulas que vale memorizar:
A segunda diz algo útil: a reta ajustada sempre passa por \((\bar{x}, \bar{y})\). Faça o que fizer, ela pivota em torno do centro de massa dos dados.
Feito à mão
Cinco pontos, pequenos o bastante para conferir cada passo:
| \(x_i\) | \(y_i\) | \(x_i - \bar{x}\) | \(y_i - \bar{y}\) | produto | \((x_i-\bar{x})^2\) |
|---|---|---|---|---|---|
| 1 | 2 | −2 | −2,2 | 4,4 | 4 |
| 2 | 4 | −1 | −0,2 | 0,2 | 1 |
| 3 | 5 | 0 | 0,8 | 0,0 | 0 |
| 4 | 4 | 1 | −0,2 | −0,2 | 1 |
| 5 | 6 | 2 | 1,8 | 3,6 | 4 |
| 8,0 | 10 |
Com \(\bar{x}=3\) e \(\bar{y}=4{,}2\):
Agora os resíduos \(e_i = y_i - \hat{y}_i\):
| \(x_i\) | \(y_i\) | \(\hat{y}_i\) | \(e_i\) |
|---|---|---|---|
| 1 | 2 | 2,6 | −0,6 |
| 2 | 4 | 3,4 | 0,6 |
| 3 | 5 | 4,2 | 0,8 |
| 4 | 4 | 5,0 | −1,0 |
| 5 | 6 | 5,8 | 0,2 |
Duas coisas valem para esses números e valem para todo ajuste OLS com intercepto:
Confira: \(-0{,}6+0{,}6+0{,}8-1{,}0+0{,}2 = 0\). ✓
Não são coincidências, elas são as equações normais — uma por parâmetro. Os resíduos são obrigados a ser ortogonais à coluna do intercepto e a cada coluna de atributo. Isso tem uma leitura geométrica: \(\hat{y}\) é a projeção de \(y\) sobre o espaço gerado pelas colunas de \(X\) e o vetor de resíduos é o que sobra, perpendicular a esse espaço. Mínimos quadrados é um problema de ângulo reto disfarçado.
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(X_train, y_train)
model.coef_, model.intercept_
y_pred = model.predict(X_test)
Quando a forma fechada tem dificuldade
Inverter \(X^\top X\) custa \(O(d^3)\) e falha de vez quando os atributos são perfeitamente colineares — a matriz é singular e não há resposta única. Para problemas muito largos ou malcondicionados, passamos para a rota iterativa, que é onde esta aula termina.
Resíduos
Tudo o que o modelo não conseguiu explicar está nos resíduos e lê-los é o hábito de diagnóstico mais útil em regressão.
O gráfico de resíduos põe \(e_i\) no eixo vertical contra o previsto \(\hat{y}_i\) no horizontal. Um gráfico saudável parece não ter nada: uma faixa sem estrutura em torno de zero. Qualquer padrão é o modelo avisando qual suposição acabou de ser quebrada.
Tente diagnosticar os quatro casos abaixo antes de ler o veredicto — o painel de cima é o que você normalmente olharia e a questão é justamente que ele não basta:
| O que você vê | O que significa | O que fazer |
|---|---|---|
| Faixa sem estrutura | A forma linear serve | Nada |
| Curva em U ou ∩ | Relação não linear que a reta não acompanha | Adicionar termos polinomiais, ou transformar \(x\) |
| Funil (dispersão cresce) | Heteroscedasticidade — o ruído depende de \(\hat{y}\) | Transformar \(y\) (muitas vezes \(\log\)), ou mínimos quadrados ponderados |
| Ponto extremo isolado | Outlier ou observação de alta alavancagem | Investigar; nunca apagar em silêncio |
| Ondas / deriva ao longo do índice | Erros correlacionados, tipicamente tempo | Modelo de série temporal; as barras de erro do OLS são inválidas |
O R² não enxerga nada disso
No simulador os quatro conjuntos são ajustados pelo mesmo procedimento e vários têm R² respeitável. O quarteto de Anscombe leva o argumento ao extremo: quatro conjuntos com médias, variâncias, correlação, reta de regressão e R² idênticos, que não se parecem em nada. Estatísticas-resumo comprimem; gráficos de resíduos não. Plote-os.
Outlier e alavancagem são coisas diferentes. Um outlier tem resíduo grande — o modelo errou nele. Um ponto de alta alavancagem fica longe de \(\bar{x}\) ao longo do eixo \(x\) e consegue arrastar a reta inteira para si. O caso perigoso é um ponto com alta alavancagem e puxão: ele move a reta com tanta eficácia que o próprio resíduo acaba pequeno, escondendo o estrago. É o quarto caso do simulador. Por isso "basta descartar os resíduos grandes" é um mau conselho — os piores infratores não têm resíduo grande.
Medindo a qualidade
Toda métrica abaixo é um jeito diferente de resumir o mesmo vetor de resíduos \(e_i = y_i - \hat{y}_i\) num único número. Elas discordam porque comprimem de formas diferentes, então vale ler os nomes por extenso antes das fórmulas.
| Métrica | Significa | Fórmula | Unidade |
|---|---|---|---|
| MAE | Mean Absolute Error — erro absoluto médio | \(\frac{1}{n}\sum \lvert e_i \rvert\) | a mesma de \(y\) |
| MSE | Mean Squared Error — erro quadrático médio | \(\frac{1}{n}\sum e_i^2\) | \(y\) ao quadrado |
| RMSE | Root Mean Squared Error — raiz do erro quadrático médio | \(\sqrt{\text{MSE}}\) | a mesma de \(y\) |
| R² | coeficiente de determinação | \(1 - \frac{\sum e_i^2}{\sum (y_i - \bar{y})^2}\) | nenhuma (é razão) |
As siglas não são arbitrárias: cada uma soletra a própria receita, lida da direita para a esquerda. RMSE é Error → Squared → Mean → Root: pegue cada erro, eleve ao quadrado, tire a média e então extraia a raiz. Ler o nome de trás para frente é escrever o código.
Calculando as quatro nos cinco pontos
Os resíduos do exemplo feito à mão eram \(-0{,}6,\ 0{,}6,\ 0{,}8,\ -1{,}0,\ 0{,}2\).
Leia esses números no contexto: os valores do alvo iam de 2 a 6, então errar cerca de 0,7 em média é um erro relevante, não um detalhe de arredondamento. MAE e RMSE são citáveis para quem não é técnico — "o modelo erra uns 0,7" — enquanto o MSE é 0,48 unidades ao quadrado, que não é uma frase que alguém consiga dizer em voz alta. É exatamente para isso que o RMSE existe: desfaz o quadrado e devolve o número à escala daquilo que você está prevendo.
MAE ou RMSE?
Não são intercambiáveis e a diferença está justamente em como tratam um erro grande. Dois modelos, cinco previsões cada:
| erros | MAE | RMSE | |
|---|---|---|---|
| Modelo A | 2, 2, 2, 2, 2 | 2,0 | 2,0 |
| Modelo B | 0, 0, 0, 0, 10 | 2,0 | 4,47 |
MAE idêntico e RMSE mais que o dobro no B. O MAE diz que os dois modelos são igualmente bons; o RMSE diz que o B é bem pior, porque eleva ao quadrado antes de tirar a média e o único erro de 10 contribui com 100 à soma. Nenhum dos dois está certo no abstrato — a pergunta é o que o seu problema custa. Errar o horário de uma entrega em 10 minutos uma vez costuma ser pior que errar 2 minutos cinco vezes (use RMSE); numa previsão de demanda o erro total pode ser o que importa (use MAE). Escolha antes de olhar os resultados, não depois.
O R² e por que esse nome
O R² compara o seu modelo com o mais preguiçoso possível. Ele é feito de duas somas. A primeira é a SSE (sum of squared errors, a soma dos quadrados dos erros) que a sua reta ainda comete. A segunda é a SST (total sum of squares, a soma total de quadrados), que é o erro de ignorar \(x\) por completo.
A SST, a soma total de quadrados, é o erro que você cometeria ignorando \(x\) por completo e prevendo sempre a média \(\bar{y}\). A SSE, a soma de quadrados dos erros, é o que a sua reta ainda erra. Então o R² responde: da variação que havia para explicar, que fração o modelo removeu? Daí "coeficiente de determinação" — o quanto de \(y\) fica determinado pelos atributos.
Nos cinco pontos: \(R^2 = 1 - 2{,}40/8{,}80 = 0{,}727\), ou seja, a reta removeu cerca de 73% da dispersão original. E como \(r_{xy} = 0{,}853\) e \(0{,}853^2 = 0{,}727\), note que na regressão simples o R² é mesmo a correlação ao quadrado — é de onde vem o símbolo. Deixa de valer no instante em que você acrescenta um segundo atributo.
O R² pode ser negativo em dados de teste, o que surpreende quem o aprendeu como "uma porcentagem". Significa que o modelo foi pior que a reta horizontal \(\hat{y}=\bar{y}\). Sempre ajuste essa baseline constante (DummyRegressor) e reporte-a: é desconfortável com que frequência um modelo elaborado mal a supera.
Avalie em dados separados
Toda métrica acima só tem sentido em dados que o modelo não viu. O R² no conjunto de treino nunca diminui quando você adiciona um atributo, mesmo uma coluna de puro ruído — o que o torna inútil para escolher entre modelos. Esse é o tema de Validação & Vazamento de Dados.
Suposições por trás das inferências
As previsões do OLS exigem muito pouco. Interpretar coeficientes, intervalos de confiança e p-valores se apoia nas suposições clássicas:
- Linearidade — a relação verdadeira é aproximadamente linear nos atributos;
- Independência — os resíduos não são correlacionados entre si (cuidado com séries temporais);
- Homoscedasticidade — a variância dos resíduos é constante ao longo da faixa de \(\hat{y}\);
- Normalidade dos resíduos — necessária para intervalos e p-valores exatos, não para o ajuste em si;
- Sem multicolinearidade severa — atributos altamente correlacionados tornam os coeficientes individuais instáveis.
A multicolinearidade merece uma frase própria, porque seu sintoma é contraintuitivo. Com dois atributos correlacionados a \(r = 0{,}99\), os dados mal conseguem distinguir seus coeficientes: \(X^\top X\) fica quase singular, então, se a verdade é \((w_1, w_2) = (1, 1)\), o par \((2, 0)\) ajusta com cerca de 2% de diferença e \((3, -1)\) com cerca de 8%. Reajuste numa amostra nova e as estimativas individuais se deslocam uma unidade inteira ou mais — o bastante para trocar de sinal — enquanto a soma delas fica cravada perto de 2 e as previsões junto. O modelo está bem; a interpretação de qualquer um dos coeficientes isolado não vale nada.
Da resposta exata ao gradiente
Temos uma fórmula que resolve o problema de uma vez. Então por que o resto do aprendizado de máquina se dá ao trabalho de iterar?
Porque \((X^\top X)^{-1}\) deixa de estar disponível. Custa \(O(d^3)\), então um milhão de atributos é inviável. Não existe quando as colunas são colineares. E não tem equivalente algum para os modelos da segunda metade do curso — uma rede neural não tem equações normais. O que generaliza é a ideia por baixo.
Repare em como chegamos à forma fechada: escrevemos \(J(w)\), tomamos o gradiente e perguntamos onde ele se anula. O gradiente descendente mantém os dois primeiros passos e desiste do terceiro. Em vez de resolver \(\nabla J = 0\) algebricamente, ele desce a ladeira até lá.
Para a regressão simples, escreva a perda como média, para que a escala não dependa de \(n\):
Derive em relação a cada parâmetro, pela regra da cadeia, lembrando que \(e_i = y_i - w_0 - w_1x_i\):
Duas coisas saem daí e as duas merecem uma pausa.
Primeira: iguale ambas a zero e você recupera exatamente as duas identidades que o exemplo feito à mão satisfazia: \(\sum e_i = 0\) e \(\sum e_i x_i = 0\). A forma fechada e o gradiente são a mesma afirmação, abordada por dois caminhos.
Segunda: veja do que o gradiente é feito — uma soma ponderada pelos resíduos. Pontos que o modelo já prevê bem contribuem com \(e_i \approx 0\) e quase não votam. Os pontos em que ele erra feio dominam o passo. A regra de atualização
lê-se, então, em português claro: mova cada parâmetro na direção para onde os erros apontam, numa quantidade proporcional ao tamanho desses erros. É a frase que continuará verdadeira, sem mudar uma vírgula, quando o modelo for uma rede de 100 camadas e \(\nabla J\) vier da retropropagação.
A figura abaixo torna a geometria concreta. À esquerda está \(J(w_0, w_1)\) como uma paisagem — cada ponto é uma reta candidata, a cor é o erro dela e a cruz branca é o ótimo em forma fechada. Clique para soltar um palpite inicial e descer:
Agora ligue centrar x e rode de novo a partir de um ponto parecido. Mesmos dados, mesmo algoritmo, comportamento radicalmente diferente.
Sem centrar, as curvas de nível formam um vale longo e estreito: \(w_0\) e \(w_1\) ficam fortemente acoplados, porque mexer na inclinação de uma reta cujos valores de \(x\) estão todos longe de zero também balança a altura dela. O gradiente aponta atravessado ao vale e não ao longo dele, então o caminho ziguezagueia e se arrasta. Centrado, o vale vira uma tigela e os mesmos passos vão quase direto à resposta.
É por isso que escalonar vem antes de ajustar
A forma fechada não liga a nada disso — \((X^\top X)^{-1}X^\top y\) devolve a mesma reta ajustada nos dois casos. A rota iterativa liga e muito. Essa é a ligação entre pré-processamento e otimização: padronizar atributos não é arrumação cosmética, é remodelar a superfície que o otimizador tem de percorrer.
Escolher \(\eta\), fazer isso em lotes em vez do conjunto inteiro e acrescentar a \(J\) penalidades que encolhem os pesos, são o tema da próxima aula: Gradiente Descendente & Regularização.
Material de aula
Roteiros
Quatro materiais de apoio, listados na página de handouts:
- Regressão linear: da reta ao gradiente — o roteiro de aula: Galton, mínimos quadrados à mão, resíduos e o gradiente, com simuladores ao vivo e cinco etapas de laboratório;
- Regressão linear em profundidade — o estudo aprofundado: a geometria da projeção, a álgebra das equações normais, alavancagem e influência, inferência — e onde cada suposição realmente aperta.
- Regressão linear do zero — notebook no Colab que constrói a conta inteira, de ŷ = Xw ao muro de memória e ao gradiente descendente à mão;
- Para casa: Python que roda na página — nove células executáveis via Pyodide, só biblioteca padrão, fechando com oito perguntas e respostas comentadas.
Referências
- Legendre, A.-M. "Nouvelles méthodes pour la détermination des orbites des comètes." Courcier (1805) — apêndice Sur la Méthode des moindres quarrés**. texto completo
- Gauss, C. F. "Theoria motus corporum coelestium in sectionibus conicis solem ambientium." Perthes & Besser (1809). texto completo
- Galton, F. "Regression towards Mediocrity in Hereditary Stature." J. Anthropological Institute 15 (1886). texto completo
- Hanley, J. A. ""Transmuting" Women into Men: Galton's Family Data on Human Stature." The American Statistician 58 (2004). DOI
- Anscombe, F. J. "Graphs in Statistical Analysis." The American Statistician 27 (1973). DOI
- Cook, R. D. "Detection of Influential Observation in Linear Regression." Technometrics 19 (1977). DOI
Bibliografia completa do curso na página de referências.