Regressão Linear Simples: Mínimos Quadrados, Máxima Verossimilhança e Inferência Bayesiana

library(tidyverse)
library(MASS)
library(patchwork)

## Tema global para todos os graficos
theme_set(theme_bw(base_size = 14))

1 Objetivos

Este documento apresenta a estimação dos parâmetros de um modelo de regressão linear simples sob três abordagens: mínimos quadrados, máxima verossimilhança e inferência bayesiana. Para cada uma delas, derivamos as estimativas pontuais dos parâmetros e as respectivas medidas de incerteza (variâncias, distribuições posteriores). Ao final, todas as soluções são verificadas numericamente por meio de um exemplo com dados simulados.

2 Definição do problema

O modelo de regressão linear simples relaciona uma variável resposta \(y\) a uma variável preditiva numérica \(x\) por meio de

\[ y_i = \beta_0 + \beta_1 x_i + w_i, \quad i = 1, 2, \ldots, n \]

onde \(\beta_0\) é o intercepto, \(\beta_1\) é a inclinação, e \(w_i\) é o termo de erro. Cada valor fixo de \(x_i\) determina uma distribuição de probabilidade para \(y_i\) com média \(\beta_0 + \beta_1 x_i\) e variância \(\sigma^2\):

\[ E[y_i] = \beta_0 + \beta_1 x_i, \quad \text{Var}[y_i] = \sigma^2 \]

Os erros \(w_i = y_i - \beta_0 - \beta_1 x_i\) são assumidos permutáveis, com \(E[w_i] = 0\) e \(\text{Var}[w_i] = \sigma^2\). Além disso, admitimos que \(y_i | x_i \sim N(\beta_0 + \beta_1 x_i, \sigma^2)\) e que as observações são independentes. O problema consiste em estimar os parâmetros desconhecidos \(\beta_0\), \(\beta_1\) e \(\sigma^2\), obter medidas de incerteza e fazer inferência.

3 Desenvolvimento

3.1 Notação e estatísticas auxiliares

Seguindo a convenção de Kinas & Andrade, definimos as médias amostrais com \(n\) no denominador:

\[ \bar{x} = \frac{\sum_{i=1}^n x_i}{n}, \quad \bar{y} = \frac{\sum_{i=1}^n y_i}{n}, \quad \overline{xx} = \frac{\sum_{i=1}^n x_i^2}{n}, \quad \overline{yy} = \frac{\sum_{i=1}^n y_i^2}{n}, \quad \overline{xy} = \frac{\sum_{i=1}^n x_i y_i}{n} \]

As variâncias e covariância amostrais (com \(n\) no denominador):

\[ S_{xx} = \overline{xx} - \bar{x}^2, \quad S_{yy} = \overline{yy} - \bar{y}^2, \quad S_{xy} = \overline{xy} - \bar{x}\bar{y} \]

O coeficiente de correlação e a variância residual (corrigida por \(n-2\) graus de liberdade):

\[ r = \frac{S_{xy}}{\sqrt{S_{xx} S_{yy}}}, \quad S_e^2 = \frac{n \, S_{yy}(1 - r^2)}{n - 2} \]

O uso de \(n\) no denominador das variâncias é uma simplificação notacional. Essa escolha não interfere nas estimativas finais, pois as fórmulas derivadas já incorporam as correções necessárias. Por exemplo, \(n \cdot S_{xx} = \sum(x_i - \bar{x})^2\) recupera a soma de quadrados sem divisão por qualquer fator. Essa compensação aparece sistematicamente nas expressões para \(S_e^2\), \(S_{\beta_0}\), \(S_{\beta_1}\) e \(V_\beta\).

3.2 Estimativas de mínimos quadrados

3.2.1 Derivação das estimativas

O método dos mínimos quadrados (MQ) consiste em encontrar os valores de \(\beta_0\) e \(\beta_1\) que minimizam a soma de quadrados dos resíduos:

\[ SQR(\beta_0, \beta_1) = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2 \]

Derivando em relação a \(\beta_0\):

\[ \frac{\partial \, SQR}{\partial \beta_0} = -2 \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i) = 0 \implies \sum y_i = n\beta_0 + \beta_1 \sum x_i \]

Derivando em relação a \(\beta_1\):

\[ \frac{\partial \, SQR}{\partial \beta_1} = -2 \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i) x_i = 0 \implies \sum x_i y_i = \beta_0 \sum x_i + \beta_1 \sum x_i^2 \]

Da primeira equação, \(\beta_0 = \bar{y} - \beta_1 \bar{x}\). Substituindo na segunda e resolvendo para \(\beta_1\):

\[ b_1 = \frac{\sum x_i y_i - n\bar{x}\bar{y}} {\sum x_i^2 - n\bar{x}^2} = \frac{S_{xy}}{S_{xx}}, \quad b_0 = \bar{y} - b_1 \bar{x} \]

3.2.2 Esperança dos estimadores

Para verificar que os estimadores são não-viesados, expressamos \(b_1\) como função linear dos \(y_i\). Definindo \(c_i = (x_i - \bar{x})/\sum(x_j - \bar{x})^2\), podemos escrever:

\[ b_1 = \frac{\sum(x_i - \bar{x})(y_i - \bar{y})} {\sum(x_i - \bar{x})^2} = \sum_{i=1}^n c_i \, y_i \]

onde usamos o fato de que \(\sum c_i = 0\) (o que elimina o termo \(\bar{y}\)). Aplicando a esperança e usando \(E[y_i] = \beta_0 + \beta_1 x_i\):

\[ E[b_1] = \sum c_i (\beta_0 + \beta_1 x_i) = \beta_0 \underbrace{\sum c_i}_{0} + \beta_1 \underbrace{\sum c_i x_i}_{1} = \beta_1 \]

onde \(\sum c_i x_i = \sum(x_i - \bar{x})x_i / \sum(x_i - \bar{x})^2 = 1\). Para o intercepto:

\[ E[b_0] = E[\bar{y} - b_1 \bar{x}] = (\beta_0 + \beta_1 \bar{x}) - E[b_1]\bar{x} = \beta_0 \]

Ambos os estimadores são, portanto, não-viesados.

3.2.3 Variâncias dos estimadores

Como as observações \(y_i\) são independentes com variância \(\sigma^2\) e \(b_1 = \sum c_i y_i\):

\[ \text{Var}[b_1] = \sum c_i^2 \, \text{Var}[y_i] = \sigma^2 \sum c_i^2 = \frac{\sigma^2}{\sum(x_i - \bar{x})^2} = \frac{\sigma^2}{n \, S_{xx}} \]

Para o intercepto, usamos \(b_0 = \bar{y} - b_1 \bar{x}\) e a independência entre \(\bar{y}\) e \(b_1\) (condicional em \(x\), são funções lineares de \(y\) com covariância nula — o que pode ser verificado diretamente). Assim:

\[ \text{Var}[b_0] = \text{Var}[\bar{y}] + \bar{x}^2 \, \text{Var}[b_1] - 2\bar{x}\,\text{Cov}[\bar{y}, b_1] \]

Os três termos são:

\[ \text{Var}[\bar{y}] = \frac{\sigma^2}{n}, \quad \text{Cov}[\bar{y}, b_1] = \text{Cov}\!\left[\frac{1}{n}\sum y_i, \sum c_j y_j\right] = \frac{\sigma^2}{n}\sum c_i = 0 \]

Logo:

\[ \text{Var}[b_0] = \frac{\sigma^2}{n} + \frac{\bar{x}^2 \sigma^2} {n \, S_{xx}} = \frac{\sigma^2 (\overline{xx})}{n \, S_{xx}} \]

onde usamos \(1/n + \bar{x}^2/(n \, S_{xx}) = (S_{xx} + \bar{x}^2)/(n \, S_{xx}) = \overline{xx}/(n \, S_{xx})\). A covariância entre os estimadores é:

\[ \text{Cov}[b_0, b_1] = \text{Cov}[\bar{y} - b_1\bar{x},\; b_1] = -\bar{x}\,\text{Var}[b_1] = \frac{-\sigma^2 \bar{x}}{n \, S_{xx}} \]

3.2.4 Distribuição conjunta e matriz de covariâncias

Sob a suposição de normalidade dos erros, os estimadores \(b_0\) e \(b_1\), sendo combinações lineares de variáveis normais independentes, têm distribuição normal conjunta:

\[ \begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N_2 \left( \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix}, \; \sigma^2 V_\beta \right), \quad V_\beta = \frac{1}{n^2 S_{xx}} \begin{pmatrix} n\overline{xx} & -n\bar{x} \\ -n\bar{x} & n \end{pmatrix} \]

Como \(\sigma^2\) é desconhecido, substituímos por \(S_e^2\). Com essa substituição, as distribuições marginais de \(b_0\) e \(b_1\) deixam de ser normais e passam a ser \(t\) de Student com \(n - 2\) graus de liberdade.

3.2.5 Estimador da variância

O estimador não-viesado de \(\sigma^2\) é

\[ S_e^2 = \frac{\sum_{i=1}^n (y_i - b_0 - b_1 x_i)^2}{n - 2} = \frac{SQR}{n - 2} \]

O divisor \(n - 2\) (ao invés de \(n\)) reflete a perda de dois graus de liberdade na estimação de \(\beta_0\) e \(\beta_1\). A distribuição de \(S_e^2\) está ligada à qui-quadrado:

\[ \frac{(n-2) S_e^2}{\sigma^2} \sim \chi^2_{n-2} \]

De onde decorre imediatamente a esperança e a variância de \(S_e^2\). Usando \(E[\chi^2_k] = k\) e \(\text{Var}[\chi^2_k] = 2k\):

\[ E\!\left[\frac{(n-2) S_e^2}{\sigma^2}\right] = n - 2 \implies E[S_e^2] = \sigma^2 \]

\[ \text{Var}\!\left[\frac{(n-2) S_e^2}{\sigma^2}\right] = 2(n-2) \implies \text{Var}[S_e^2] = \frac{2\sigma^4}{n - 2} \]

3.3 Estimativas de máxima verossimilhança

3.3.1 Função de verossimilhança

Sob a suposição de normalidade, a função de verossimilhança para \((\beta_0, \beta_1, \sigma^2)\) é

\[ L(\beta_0, \beta_1, \sigma^2) = \left(\frac{1}{2\pi\sigma^2}\right)^{n/2} \exp\!\left[-\frac{1}{2\sigma^2} \sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2\right] \]

e a log-verossimilhança:

\[ \ell(\beta_0, \beta_1, \sigma^2) = -\frac{n}{2}\log(2\pi) -\frac{n}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2 \]

3.3.2 Derivação das estimativas

As derivadas parciais de \(\ell\) igualadas a zero fornecem:

\[ \frac{\partial \ell}{\partial \beta_0} = \frac{1}{\sigma^2}\sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i) = 0 \]

\[ \frac{\partial \ell}{\partial \beta_1} = \frac{1}{\sigma^2}\sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i) x_i = 0 \]

\[ \frac{\partial \ell}{\partial \sigma^2} = -\frac{n}{2\sigma^2} + \frac{1}{2\sigma^4}\sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2 = 0 \]

Resolvendo:

\[ b_1 = \frac{S_{xy}}{S_{xx}}, \quad b_0 = \bar{y} - b_1 \bar{x}, \quad \hat{\sigma}^2_{MV} = \frac{1}{n} \sum_{i=1}^{n}(y_i - b_0 - b_1 x_i)^2 \]

As estimativas de \(\beta_0\) e \(\beta_1\) são idênticas às de mínimos quadrados. Isso decorre do fato de que, para \(\sigma^2\) fixo, maximizar a verossimilhança equivale a minimizar \(\sum(y_i - \beta_0 - \beta_1 x_i)^2\).

3.3.3 Viés do estimador de variância

O estimador de MV da variância utiliza \(n\) no denominador:

\[ \hat{\sigma}^2_{MV} = \frac{SQR}{n} \]

Para obter sua esperança, usamos o resultado da seção anterior: \((n-2)S_e^2/\sigma^2 \sim \chi^2_{n-2}\), onde \(S_e^2 = SQR/(n-2)\). Como \(SQR = n \, \hat{\sigma}^2_{MV}\):

\[ E[\hat{\sigma}^2_{MV}] = E\!\left[\frac{SQR}{n}\right] = \frac{(n-2)}{n} E[S_e^2] = \frac{n-2}{n}\sigma^2 \]

Portanto, o EMV é viesado, subestimando \(\sigma^2\) em média. A versão corrigida é \(S_e^2 = \frac{n}{n-2}\hat{\sigma}^2_{MV}\).

3.3.4 Variâncias assintóticas via informação de Fisher

As segundas derivadas da log-verossimilhança fornecem a matriz de informação de Fisher. Para o bloco \((\beta_0, \beta_1)\):

\[ \frac{\partial^2 \ell}{\partial \beta_0^2} = -\frac{n}{\sigma^2}, \quad \frac{\partial^2 \ell}{\partial \beta_1^2} = -\frac{\sum x_i^2}{\sigma^2}, \quad \frac{\partial^2 \ell}{\partial \beta_0 \partial \beta_1} = -\frac{\sum x_i}{\sigma^2} \]

A informação de Fisher observada (negativo das segundas derivadas) é

\[ I_O(\beta_0, \beta_1) = \begin{pmatrix} n/\sigma^2 & n\bar{x}/\sigma^2 \\ n\bar{x}/\sigma^2 & n\overline{xx}/\sigma^2 \end{pmatrix} \]

Nenhuma entrada depende de \(y\), de modo que \(I_E = I_O\). A inversa dessa matriz é exatamente \(\sigma^2 V_\beta\), recuperando as variâncias obtidas por mínimos quadrados. Isso demonstra que, neste caso, as variâncias assintóticas (via informação de Fisher) coincidem com as variâncias exatas. Esta coincidência é específica do modelo normal linear; em modelos mais gerais, a informação de Fisher fornece apenas uma aproximação assintótica.

Para \(\sigma^2\), a segunda derivada é

\[ \frac{\partial^2 \ell}{\partial (\sigma^2)^2} = \frac{n}{2(\sigma^2)^2} - \frac{1}{(\sigma^2)^3}\sum(y_i - \beta_0 - \beta_1 x_i)^2 \]

Tomando a esperança e invertendo, obtém-se a variância assintótica \(\text{Var}[\hat{\sigma}^2_{MV}] \approx 2\sigma^4/n\).

3.4 Abordagem bayesiana

3.4.1 A priori de Jeffreys: definição geral

Na abordagem bayesiana, os parâmetros são tratados como variáveis aleatórias com distribuições a priori. Quando não há informação prévia substancial, é desejável utilizar uma priori que seja “minimamente informativa”. A priori de Jeffreys atende a esse critério por ser invariante a reparametrizações: se mudarmos a escala do parâmetro (e.g. de \(\sigma^2\) para \(\log \sigma^2\)), a priori se transforma de maneira coerente.

Para um vetor de parâmetros \(\boldsymbol{\theta}\), a priori de Jeffreys é definida como

\[ p(\boldsymbol{\theta}) \propto \sqrt{\det\!\left[I_E(\boldsymbol{\theta})\right]} \]

onde \(I_E(\boldsymbol{\theta})\) é a matriz de informação de Fisher esperada, cujo elemento \((j, k)\) é

\[ [I_E(\boldsymbol{\theta})]_{jk} = -E\!\left[ \frac{\partial^2 \ell(\boldsymbol{\theta})} {\partial \theta_j \, \partial \theta_k} \right] \]

No caso escalar (\(\theta\) unidimensional), a expressão se reduz a

\[ p(\theta) \propto \sqrt{I_E(\theta)} = \sqrt{-E\!\left[ \frac{\partial^2 \ell(\theta)}{\partial \theta^2} \right]} \]

3.4.2 Priori de Jeffreys para o modelo de regressão

Para o modelo \(y_i \sim N(\beta_0 + \beta_1 x_i, \sigma^2)\), o vetor de parâmetros é \(\boldsymbol{\theta} = (\beta_0, \beta_1, \sigma^2)\). A informação de Fisher é bloco-diagonal entre \((\beta_0, \beta_1)\) e \(\sigma^2\), o que permite tratar cada bloco separadamente.

Bloco de \((\beta_0, \beta_1)\): A sub-matriz de informação de Fisher é

\[ I_E(\beta_0, \beta_1) = \frac{1}{\sigma^2} \begin{pmatrix} n & n\bar{x} \\ n\bar{x} & n\overline{xx} \end{pmatrix} \]

O determinante dessa matriz é

\[ \det[I_E] = \frac{1}{\sigma^4} \left(n \cdot n\overline{xx} - (n\bar{x})^2\right) = \frac{n^2}{\sigma^4}\left(\overline{xx} - \bar{x}^2\right) = \frac{n^2 S_{xx}}{\sigma^4} \]

Todas as quantidades: \(n\), \(S_{xx}\), \(\overline{xx}\), \(\bar{x}\), dependem exclusivamente dos valores observados de \(x\) (tratados como fixos no modelo de regressão) e de \(\sigma^2\). Nenhuma depende de \(\beta_0\) ou \(\beta_1\). Portanto, \(\sqrt{\det[I_E]}\) é constante como função de \((\beta_0, \beta_1)\), o que resulta em \(p(\beta_0, \beta_1) \propto 1\).

A priori de Jeffreys para os coeficientes de regressão é uniforme (imprópria).

Bloco de \(\sigma^2\): A segunda derivada da log-verossimilhança em relação a \(\sigma^2\) é

\[ \frac{\partial^2 \ell}{\partial (\sigma^2)^2} = \frac{n}{2(\sigma^2)^2} - \frac{1}{(\sigma^2)^3}\sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2 \]

Tomando a esperança (usando \(E\!\left[\sum(y_i - \beta_0 - \beta_1 x_i)^2\right] = n\sigma^2\)):

\[ I_E(\sigma^2) = -E\!\left[\frac{\partial^2 \ell}{\partial (\sigma^2)^2}\right] = -\frac{n}{2(\sigma^2)^2} + \frac{n\sigma^2}{(\sigma^2)^3} = \frac{n}{2(\sigma^2)^2} \]

Aplicando a definição:

\[ p(\sigma^2) \propto \sqrt{I_E(\sigma^2)} = \sqrt{\frac{n}{2(\sigma^2)^2}} \propto \frac{1}{\sigma^2} \]

Em resumo, a priori de Jeffreys para o modelo de regressão linear simples é

\[ p(\beta_0, \beta_1, \sigma^2) \propto \frac{1}{\sigma^2} \]

3.4.3 Derivação das distribuições posteriores

3.4.3.1 Posterior conjunta

Pelo teorema de Bayes:

\[ p(\beta_0, \beta_1, \sigma^2 \mid y) \propto L(\beta_0, \beta_1, \sigma^2) \cdot p(\beta_0, \beta_1, \sigma^2) \]

A posterior conjunta fica:

\[ p(\boldsymbol{\beta}, \sigma^2 \mid y) \propto (\sigma^2)^{-n/2 - 1} \exp\!\left[ -\frac{(n-2)S_e^2 + Q_\beta}{2\sigma^2} \right] \]

O termo \(Q_\beta\) é uma forma quadrática:

\[ Q_\beta = (\boldsymbol{\beta} - \mathbf{b})^T V_\beta^{-1} (\boldsymbol{\beta} - \mathbf{b}) \]

onde \(\mathbf{b} = (b_0, b_1)^T\) e \(V_\beta^{-1} = \begin{pmatrix} n & n\bar{x} \\ n\bar{x} & n\overline{xx} \end{pmatrix}\).

3.4.3.2 Posterior marginal de \(\sigma^2\) (Gama-Inversa)

Integrando em \((\beta_0, \beta_1)\):

\[ p(\sigma^2 \mid y) \propto (\sigma^2)^{-n/2 - 1} \exp\!\left[-\frac{(n-2)S_e^2}{2\sigma^2}\right] \int_{\mathbb{R}^2} \exp\!\left[-\frac{Q_\beta}{2\sigma^2}\right] d\boldsymbol{\beta} \]

O resultado e o núcleo de uma Gama-Inversa:

\[ \boxed{\; \sigma^2 \mid y \sim \text{Gama-Inv}\!\left(\frac{n-2}{2},\;\frac{(n-2)S_e^2}{2}\right) \;} \]

3.4.3.3 Posterior marginal de \(\boldsymbol{\beta}\) (Student bivariada)

Integrando a posterior conjunta em \(\sigma^2\):

\[ p(\boldsymbol{\beta} \mid y) \propto \int_0^\infty (\sigma^2)^{-n/2 - 1} \exp\!\left[-\frac{(n-2)S_e^2 + Q_\beta}{2\sigma^2}\right] d\sigma^2 \]

O resultado é o núcleo de uma \(t\) de Student bivariada:

\[ \boxed{\; \boldsymbol{\beta} \mid y \sim t_2\!\left(n - 2,\; \mathbf{b},\; S_e^2 V_\beta\right) \;} \]

3.4.3.4 Posteriores marginais univariadas

As marginais de uma \(t\) multivariada são \(t\) univariadas. Para cada \(\beta_j\), a localização é \(b_j\) e a escala é a raiz do elemento diagonal correspondente de \(S_e^2 V_\beta\):

\[ \boxed{\; \begin{aligned} \beta_0 \mid y &\sim \text{St}\!\left(n - 2,\; b_0,\; S_{\beta_0}\right) \\ \beta_1 \mid y &\sim \text{St}\!\left(n - 2,\; b_1,\; S_{\beta_1}\right) \end{aligned} \;} \]

onde

\[ S_{\beta_0} = S_e\sqrt{\frac{\overline{xx}}{n\,S_{xx}}}, \quad S_{\beta_1} = \frac{S_e}{\sqrt{n\,S_{xx}}} \]

Esses parâmetros de escala são exatamente as raízes quadradas das variâncias dos estimadores de MQ derivadas na Seção 3.2.3. De fato, lembrando que

\[ \text{Var}[b_0] = \frac{\sigma^2 \overline{xx}}{n\,S_{xx}}, \quad \text{Var}[b_1] = \frac{\sigma^2}{n\,S_{xx}} \]

e substituindo \(\sigma^2\) por sua estimativa \(S_e^2\), os erros-padrão frequentistas são

\[ \text{EP}[b_0] = \sqrt{\widehat{\text{Var}}[b_0]} = S_e\sqrt{\frac{\overline{xx}}{n\,S_{xx}}} = S_{\beta_0}, \quad \text{EP}[b_1] = \sqrt{\widehat{\text{Var}}[b_1]} = \frac{S_e}{\sqrt{n\,S_{xx}}} = S_{\beta_1} \]

Portanto, os parâmetros de escala das posteriores bayesianas são idênticos aos erros-padrão da inferência clássica — os mesmos valores que a função lm() do R reporta na coluna Std. Error. A distribuição \(t\) de Student aparece nos dois paradigmas, mas com interpretações distintas: na abordagem frequentista, descreve a variabilidade amostral do estimador; na bayesiana, expressa a incerteza posterior sobre o parâmetro.

3.5 Comparação teórica entre as três abordagens

Para este modelo com priori de Jeffreys, as três abordagens conduzem a estimativas pontuais equivalentes. As diferenças residem na interpretação e na quantificação da incerteza:

• Mínimos quadrados fornece estimativas pontuais e variâncias exatas sob normalidade. A inferência usa a distribuição \(t\) de Student como distribuição amostral.
• Máxima verossimilhança chega às mesmas estimativas pontuais. A teoria assintótica (informação de Fisher) fornece uma via geral para variâncias, útil em modelos sem solução analítica. Neste caso, as variâncias assintóticas coincidem com as exatas.
• Inferência bayesiana produz distribuições posteriores completas. Sob a priori de Jeffreys, as posteriores marginais de \(\beta_0\) e \(\beta_1\) são \(t\) de Student com os mesmos parâmetros de localização e escala da inferência clássica — mas a interpretação é probabilística.

4 Aplicação: exemplo com dados simulados

Nesta seção, verificamos numericamente todas as soluções derivadas acima, utilizando dados gerados a partir de valores conhecidos dos parâmetros. Isso permite comparar as estimativas obtidas por cada método com os valores verdadeiros.

4.1 Simulação dos dados

set.seed(2024)
n <- 100

## Parametros verdadeiros
beta0_true <- 10
beta1_true <- 0.5
sigma2_true <- 20

## Simulacao
x <- rnorm(n, mean = 150, sd = 15)
y <- beta0_true + beta1_true * x + rnorm(n, mean = 0, sd = sqrt(sigma2_true))

Fazer:

• Revisar as expressões das distribuições posteriores na abordagem bayesiana
• Ajustar o modelo linear para as 3 abordagens com os dados simulados acima
• Comparar todas as estimativas
• Verificar qual é a distribuição conjunta de \(\beta_0, \beta_1\) condicional a \(\sigma^2\).
  • Amostre desta conjunta
  • Verifique a relação entre \(\beta_0, \beta_1\)

tibble(x = x, y = y) |>
  ggplot(aes(x = x, y = y)) +
  geom_point(alpha = 0.6) +
  geom_abline(
    intercept = beta0_true, slope = beta1_true,
    colour = "red", linetype = "dashed"
  ) +
  labs(x = "x", y = "y")

Figura 1: Diagrama de dispersao dos dados simulados.

4.2 Estatísticas auxiliares

xbar <- mean(x)
ybar <- mean(y)
Sxx <- mean(x^2) - xbar^2
Syy <- mean(y^2) - ybar^2
Sxy <- mean(x * y) - xbar * ybar
r2 <- (Sxy / sqrt(Sxx * Syy))^2
Se2 <- n * Syy * (1 - r2) / (n - 2)
## SQR <- sum((y - b0 - b1*x)^2)/(n-2)
Se <- sqrt(Se2)

tibble(
  quantidade = c("barX", "barY", "Sxx", "Syy", "Sxy", "r2", "Se2", "Se"),
  valor = round(c(xbar, ybar, Sxx, Syy, Sxy, r2, Se2, Se), 4)
)

# A tibble: 8 × 2
  quantidade   valor
  <chr>        <dbl>
1 barX       149.   
2 barY        85.0  
3 Sxx        233.   
4 Syy         83.0  
5 Sxy        121.   
6 r2           0.752
7 Se2         21.0  
8 Se           4.58 

4.3 Mínimos quadrados

## Estimativas pontuais
b1 <- Sxy / Sxx
b0 <- ybar - b1 * xbar

## Variancias dos estimadores (usando Se2 no lugar de sigma2)
var_b0 <- Se2 * mean(x^2) / (n * Sxx)
var_b1 <- Se2 / (n * Sxx)
cov_b0b1 <- -Se2 * xbar / (n * Sxx)
var_Se2 <- 2 * Se2^2 / (n - 2) # 2\sigma2^2/(n-2)

## Verificacao via lm()
modelo <- lm(y ~ x)

## Compara resultados
tibble(
  parametro = c("beta0", "beta1"),
  verdadeiro = c(beta0_true, beta1_true),
  estimativa_MQ = round(c(b0, b1), 4),
  estimativa_lm = round(coef(modelo), 4),
  EP_manual = round(sqrt(c(var_b0, var_b1)), 4),
  EP_lm = round(sqrt(diag(vcov(modelo))), 4)
)

# A tibble: 2 × 6
  parametro verdadeiro estimativa_MQ estimativa_lm EP_manual EP_lm
  <chr>          <dbl>         <dbl>         <dbl>     <dbl> <dbl>
1 beta0           10           7.97          7.97       4.49  4.49
2 beta1            0.5         0.518         0.518      0.03  0.03

Note que

## Matriz de covariancia dos betas (V_beta)
summary(modelo)$cov.unscaled

             (Intercept)             x
(Intercept)  0.959819552 -6.386371e-03
x           -0.006386371  4.294051e-05

## e
summary(modelo)$sigma^2

[1] 21.00383

## é o "escalonamento" para a vcov (usando Se2)
summary(modelo)$sigma^2 * summary(modelo)$cov.unscaled

            (Intercept)             x
(Intercept)  20.1598896 -0.1341382775
x            -0.1341383  0.0009019154

## chega na vcov
vcov(modelo)

            (Intercept)             x
(Intercept)  20.1598896 -0.1341382775
x            -0.1341383  0.0009019154

## Compara
tibble(
  quantidade = c("Se2 (manual)", "Se2 (lm)", "Var(Se2)"),
  valor = round(c(Se2, summary(modelo)$sigma^2, var_Se2), 4)
)

# A tibble: 3 × 2
  quantidade   valor
  <chr>        <dbl>
1 Se2 (manual) 21.0 
2 Se2 (lm)     21.0 
3 Var(Se2)      9.00

4.4 Máxima verossimilhança

## Negativo da log-verossimilhanca
nll <- function(par, n, y, x) {
  sigma2 <- par[1]
  beta0 <- par[2]
  beta1 <- par[3]
  (n / 2) * log(2 * pi) +
    (n / 2) * log(sigma2) +
    (1 / (2 * sigma2)) * sum((y - beta0 - beta1 * x)^2)
}

## Otimizacao
mle <- optim(
  par = c(20, 10, 0.5),
  fn = nll, n = n, x = x, y = y,
  lower = c(1e-6, -Inf, -Inf),
  upper = c(Inf, Inf, Inf),
  method = "L-BFGS-B",
  hessian = TRUE
)

## Estimativas de MV
(sigma2_mv <- mle$par[1]) # viesada

[1] 20.5833

(sigma2_corr <- sigma2_mv * n / (n - 2)) # corrigida

[1] 21.00337

tibble(
  parametro = c("beta0", "beta1", "sigma2_mv", "sigma2_corr"),
  verdadeiro = c(beta0_true, beta1_true, sigma2_true, sigma2_true),
  estimativa = round(c(mle$par[2], mle$par[3], sigma2_mv, sigma2_corr), 4)
)

# A tibble: 4 × 3
  parametro   verdadeiro estimativa
  <chr>            <dbl>      <dbl>
1 beta0             10        7.97 
2 beta1              0.5      0.518
3 sigma2_mv         20       20.6  
4 sigma2_corr       20       21.0  

## Compara resultados
tibble(
    parametro = c("beta0", "beta1"),
    verdadeiro = c(beta0_true, beta1_true),
    estimativa_MQ = round(c(b0, b1), 4),
    estimativa_lm = round(coef(modelo), 4),
    EP_manual = round(sqrt(c(var_b0, var_b1)), 4),
    EP_lm = round(sqrt(diag(vcov(modelo))), 4)
)

# A tibble: 2 × 6
  parametro verdadeiro estimativa_MQ estimativa_lm EP_manual EP_lm
  <chr>          <dbl>         <dbl>         <dbl>     <dbl> <dbl>
1 beta0           10           7.97          7.97       4.49  4.49
2 beta1            0.5         0.518         0.518      0.03  0.03

tibble(
    parametro = c("beta0", "beta1", "sigma2"),
    verdadeiro = c(beta0_true, beta1_true, sigma2_true),
    MQ = round(c(b0, b1, Se2), 4),
    lm = round(c(coef(modelo), summary(modelo)$sigma^2), 4),
    MV = round(c(mle$par[2], mle$par[3], sigma2_corr), 4),
    EP_MQ = round(c(sqrt(var_b0), sqrt(var_b1), sqrt(var_Se2)), 4),
    EP_lm = round(c(sqrt(diag(vcov(modelo))), NA), 4),
    EP_MV = round(sqrt(diag(solve(mle$hessian)))[c(2, 3, 1)], 4)
)

# A tibble: 3 × 8
  parametro verdadeiro     MQ     lm     MV EP_MQ EP_lm  EP_MV
  <chr>          <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl>  <dbl>
1 beta0           10    7.97   7.97   7.97   4.49  4.49 4.44  
2 beta1            0.5  0.518  0.518  0.518  0.03  0.03 0.0297
3 sigma2          20   21.0   21.0   21.0    3.00 NA    2.91  

NOte que:

• As estimativas pontuais de MQ, lm e MV devem ser praticamente idênticas para os betas, diferindo apenas em sigma2 antes da correção de viés.
• Os erros-padrão de MQ e lm coincidem por construção. Já os de MV (via Hessiana) são uma aproximação assintótica.
• lm() não reporta erro-padrão para \(\sigma^2\) diretamente (por isso o NA), mas poderia ser calculada com

sqrt(2 * summary(modelo)$sigma^4 / (n - 2))

[1] 3.000548

4.5 Abordagem bayesiana

4.5.1 Posteriores marginais analíticas

g <- n - 2

## Parametros de escala das posteriores marginais
Sb0 <- Se * sqrt(mean(x^2) / (n * Sxx))
Sb1 <- Se / sqrt(n * Sxx)
dp_b0 <- Sb0 * sqrt(g / (g - 2))
dp_b1 <- Sb1 * sqrt(g / (g - 2))

tibble(
  parametro = c("beta0", "beta1"),
  verdadeiro = c(beta0_true, beta1_true),
  localizacao = round(c(b0, b1), 4),
  escala = round(c(Sb0, Sb1), 4),
  dp_posterior = round(c(dp_b0, dp_b1), 4),
  ICr_2.5 = round(qt(0.025, g) * c(Sb0, Sb1) + c(b0, b1), 4),
  ICr_97.5 = round(qt(0.975, g) * c(Sb0, Sb1) + c(b0, b1), 4)
)

# A tibble: 2 × 7
  parametro verdadeiro localizacao escala dp_posterior ICr_2.5 ICr_97.5
  <chr>          <dbl>       <dbl>  <dbl>        <dbl>   <dbl>    <dbl>
1 beta0           10         7.97    4.49       4.54    -0.940   16.9  
2 beta1            0.5       0.518   0.03       0.0303   0.458    0.577

## Simulando das posteriores marginais INDEPENDENTES (via rt)
set.seed(42)
m <- 5000
beta0_ind <- rt(m, df = g) * Sb0 + b0
beta1_ind <- rt(m, df = g) * Sb1 + b1

## Correlacao entre as amostras independentes (~0)
cor(beta0_ind, beta1_ind)

[1] -0.001799802

p_ind <- tibble(beta0 = beta0_ind, beta1 = beta1_ind) |>
    ggplot(aes(x = beta0, y = beta1)) +
    geom_point(alpha = 0.1, size = 0.5, colour = "grey60") +
    geom_density_2d(colour = "black", linewidth = 0.6) +
    geom_point(
        data = tibble(beta0 = beta0_true, beta1 = beta1_true),
        colour = "red", size = 3, shape = 4, stroke = 1.5
    ) +
    labs(
        x = expression(beta[0]), y = expression(beta[1]),
        title = "Marginais independentes"
    )
p_ind

## Para sigma2: Gama-Inversa via transformacao da Gama
## Se tau = 1/sigma2 ~ Gama(alpha, lambda), entao
## f(sigma2) = dgamma(1/sigma2, alpha, lambda) / sigma2^2
a_gi <- g / 2
b_gi <- g * Se2 / 2
d_sigma2 <- \(s2) dgamma(1 / s2, shape = a_gi, rate = b_gi) / s2^2
## Amostra
tau_ind <- rgamma(m, shape = g / 2, rate = g * Se2 / 2)
sigma2_ind <- 1 / tau_ind
d_sigma2_ind <- d_sigma2(sigma2_ind)

ggplot(tibble(x = sigma2_ind, y = d_sigma2_ind), aes(x = x, y = y)) +
    geom_line() +
    geom_vline(xintercept = sigma2_true, linetype = "dashed", colour = "red") +
    labs(x = expression(sigma^2), y = "Densidade")

4.5.2 Simulação da posterior conjunta

set.seed(42)
m <- 5000

## Extrair elementos do lm()
est <- coef(modelo)
v_beta <- summary(modelo)$cov.unscaled

## Amostrar tau (precisao) e sigma2
prec <- rgamma(m, g / 2, Se^2 * g / 2)
sig2_post <- 1 / prec

## Amostrar (beta0, beta1) condicionalmente em sigma2
## OBS.: Vamos amostrar de uma Normal bivariada
beta_conj <- map_dfr(sig2_post, \(s2) {
  samp <- mvrnorm(1, est, s2 * v_beta)
  tibble(beta0 = samp[1], beta1 = samp[2])
})

## Densidades posteriores analiticas
## Para beta0 e beta1: Student nao-padronizada
## f(x) = dt((x - mu) / sigma, df = g) / sigma
## onde mu = localizacao, sigma = escala (EP), g = graus de liberdade
d_beta0 <- \(val) dt((val - b0) / Sb0, df = g) / Sb0
d_beta1 <- \(val) dt((val - b1) / Sb1, df = g) / Sb1

## Para sigma2: Gama-Inversa via transformacao da Gama
## Se tau = 1/sigma2 ~ Gama(alpha, lambda), entao
## f(sigma2) = dgamma(1/sigma2, alpha, lambda) / sigma2^2
a_gi <- g / 2
b_gi <- g * Se2 / 2
d_sigma2 <- \(s2) dgamma(1 / s2, shape = a_gi, rate = b_gi) / s2^2

p1 <- ggplot(beta_conj, aes(x = beta0)) +
  geom_histogram(aes(y = after_stat(density)),
    bins = 40, fill = "grey80", colour = "white"
  ) +
  stat_function(fun = d_beta0, colour = "steelblue", linewidth = 1) +
  geom_vline(xintercept = beta0_true, linetype = "dashed", colour = "red") +
  labs(x = expression(beta[0]), y = "Densidade")

p2 <- ggplot(beta_conj, aes(x = beta1)) +
  geom_histogram(aes(y = after_stat(density)),
    bins = 40, fill = "grey80", colour = "white"
  ) +
  stat_function(fun = d_beta1, colour = "steelblue", linewidth = 1) +
  geom_vline(xintercept = beta1_true, linetype = "dashed", colour = "red") +
  labs(x = expression(beta[1]), y = "Densidade")

p3 <- ggplot(tibble(sigma2 = sig2_post), aes(x = sigma2)) +
  geom_histogram(aes(y = after_stat(density)),
    bins = 40, fill = "grey80", colour = "white"
  ) +
  stat_function(fun = d_sigma2, colour = "steelblue", linewidth = 1) +
  geom_vline(xintercept = sigma2_true, linetype = "dashed", colour = "red") +
  labs(x = expression(sigma^2), y = "Densidade")

p1 + p2 + p3 +
  plot_annotation(
    caption = "Linha vermelha: valor verdadeiro. Curva azul: posterior analitica."
  )

Figura 2: Distribuicoes posteriores marginais simuladas para beta0, beta1 e sigma2.

p_conj <- beta_conj |>
    ggplot(aes(x = beta0, y = beta1)) +
    geom_point(alpha = 0.1, size = 0.5, colour = "grey60") +
    geom_density_2d(colour = "black", linewidth = 0.6) +
    geom_point(
        data = tibble(beta0 = beta0_true, beta1 = beta1_true),
        colour = "red", size = 3, shape = 4, stroke = 1.5
    ) +
    labs(
        x = expression(beta[0]), y = expression(beta[1]),
        title = "Posterior conjunta"
    )
p_conj

p_ind + p_conj

Figura 3: Amostra posterior conjunta de (beta0, beta1).

Figura 4: Amostra posterior conjunta de (beta0, beta1).

NotaDensidades analíticas nos gráficos e propriedades da \(t\) de Student

As curvas azuis sobrepostas aos histogramas no gráfico anterior são as densidades posteriores analíticas, avaliadas ponto a ponto. Abaixo explicamos a construção de cada uma.

Posteriores de \(\beta_0\) e \(\beta_1\) (Student não-padronizada). Se \(X \sim \text{St}(g, \mu, \sigma)\) — uma \(t\) de Student com \(g\) graus de liberdade, localização \(\mu\) e escala \(\sigma\) — sua densidade é obtida a partir da \(t\) padrão por mudança de variável:

\[ f_X(x) = \frac{1}{\sigma} \, f_t\!\left(\frac{x - \mu}{\sigma};\, g\right) \]

onde \(f_t(t; g)\) é a densidade da \(t\) padrão com \(g\) graus de liberdade (a função dt() do R). Ou seja, padronizamos o argumento com \((x - \mu)/\sigma\) e dividimos por \(\sigma\) (o jacobiano da transformação). No código, isso corresponde a:

## Para beta0: mu = b0, sigma = Sb0, g graus de liberdade
\(val) dt((val - b0) / Sb0, df = g) / Sb0

Posterior de \(\sigma^2\) (Gama-Inversa). Se \(\sigma^2 \sim \text{Gama-Inv}(\alpha, \lambda)\), então \(\tau = 1/\sigma^2 \sim \text{Gama}(\alpha, \lambda)\). Para obter a densidade de \(\sigma^2\) a partir da de \(\tau\), aplicamos a transformação \(\sigma^2 = 1/\tau\) com jacobiano \(|d\tau/d\sigma^2| = 1/(\sigma^2)^2\):

\[ f_{\sigma^2}(s) = f_\tau(1/s) \cdot \frac{1}{s^2} \]

No código, usamos dgamma() avaliada em \(1/s\) e multiplicamos por \(1/s^2\):

## alpha = g/2, lambda = g * Se2 / 2
\(s2) dgamma(1 / s2, shape = a_gi, rate = b_gi) / s2^2

Escala vs. desvio padrão na \(t\) de Student. O parâmetro \(\sigma\) na parametrização \(\text{St}(g, \mu, \sigma)\) é a escala, não o desvio padrão. A variância da distribuição é

\[ \text{Var}[X] = \sigma^2 \cdot \frac{g}{g-2}, \quad g > 2 \]

e portanto o desvio padrão é

\[ \text{DP}[X] = \sigma\sqrt{\frac{g}{g-2}} \]

O fator \(g/(g-2) > 1\) reflete as caudas mais pesadas da \(t\) em relação à normal. Quando \(g \to \infty\), esse fator converge para 1 e a \(t\) se aproxima da normal, onde escala e desvio padrão coincidem. Para \(n\) pequeno, a diferença é apreciável: com \(g = 8\), por exemplo, \(\sqrt{g/(g-2)} \approx 1.15\), inflando o DP em 15% em relação à escala.

4.5.3 Outras prioris

library(rjags)

## Dados para o JAGS
dados_jags <- list(
    y = y,
    x = x,
    n = n
)

## Modelo
modelo_jags <- "model {
  ## Verossimilhanca
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1 * x[i]
  }
  ## Prioris
  beta0 ~ dnorm(0, 0.001)
  beta1 ~ dnorm(0, 0.001)
  tau ~ dgamma(0.001, 0.001)
  ## Transformacao
  sigma2 <- 1 / tau
}"

## Ajuste
jags_fit <- jags.model(
    textConnection(modelo_jags),
    data = dados_jags,
    n.chains = 3,
    n.adapt = 1000,
    quiet = TRUE
)

## Burn-in
update(jags_fit, 1000)

## Amostragem
post_jags <- coda.samples(
    jags_fit,
    variable.names = c("beta0", "beta1", "sigma2"),
    n.iter = 5000
)

summary(post_jags)

Iterations = 1001:6000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD Naive SE Time-series SE
beta0   8.0438 4.17342 0.034076       0.413397
beta1   0.5171 0.02792 0.000228       0.002762
sigma2 21.4153 3.10988 0.025392       0.026761

2. Quantiles for each variable:

           2.5%     25%     50%     75%   97.5%
beta0   0.07913  5.1833  7.9774 10.8165 16.3249
beta1   0.46132  0.4988  0.5176  0.5362  0.5707
sigma2 16.16991 19.1894 21.1339 23.3268 28.2304

plot(post_jags)

## Extrair amostras do JAGS como data.frame
post_df <- post_jags |>
  map_dfr(as_tibble)

## Grafico comparativo
p1 <- ggplot() +
  geom_density(
    data = post_df, aes(x = beta0, colour = "JAGS"),
    linewidth = 0.8
  ) +
  stat_function(
    fun = d_beta0, aes(colour = "Analitica"),
    linewidth = 0.8
  ) +
  geom_vline(xintercept = beta0_true, linetype = "dashed") +
  scale_colour_manual(values = c("Analitica" = "steelblue", "JAGS" = "tomato")) +
  labs(x = expression(beta[0]), y = "Densidade", colour = NULL)

p2 <- ggplot() +
  geom_density(
    data = post_df, aes(x = beta1, colour = "JAGS"),
    linewidth = 0.8
  ) +
  stat_function(
    fun = d_beta1, aes(colour = "Analitica"),
    linewidth = 0.8
  ) +
  geom_vline(xintercept = beta1_true, linetype = "dashed") +
  scale_colour_manual(values = c("Analitica" = "steelblue", "JAGS" = "tomato")) +
  labs(x = expression(beta[1]), y = "Densidade", colour = NULL)

p3 <- ggplot() +
  geom_density(
    data = post_df, aes(x = sigma2, colour = "JAGS"),
    linewidth = 0.8
  ) +
  stat_function(
    fun = d_sigma2, aes(colour = "Analitica"),
    linewidth = 0.8
  ) +
  geom_vline(xintercept = sigma2_true, linetype = "dashed") +
  scale_colour_manual(values = c("Analitica" = "steelblue", "JAGS" = "tomato")) +
  labs(x = expression(sigma^2), y = "Densidade", colour = NULL)

p1 + p2 + p3 +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

4.6 Comparação final

tibble(
  parametro = c("beta0", "beta1", "sigma2"),
  verdadeiro = c(beta0_true, beta1_true, sigma2_true),
  MQ = round(c(b0, b1, Se2), 4),
  MV_corrigida = round(c(mle$par[2], mle$par[3], sigma2_corr), 4),
  Bayes_media = round(c(
    mean(beta_conj$beta0),
    mean(beta_conj$beta1),
    mean(sig2_post)
  ), 4)
)

# A tibble: 3 × 5
  parametro verdadeiro     MQ MV_corrigida Bayes_media
  <chr>          <dbl>  <dbl>        <dbl>       <dbl>
1 beta0           10    7.97         7.97        8.04 
2 beta1            0.5  0.518        0.518       0.517
3 sigma2          20   21.0         21.0        21.5  

As três abordagens fornecem estimativas pontuais numericamente equivalentes, todas próximas dos valores verdadeiros utilizados na simulação. A concordância ilustra o resultado teórico: sob a priori de Jeffreys e o modelo normal linear, as soluções de MQ, MV e a inferência bayesiana convergem para os mesmos números.

4.7 Comparação entre as três abordagens

Reunimos abaixo, lado a lado, os resultados obtidos por cada método.

4.7.1 Estimativas pontuais de \(\beta_0\) e \(\beta_1\)

As três abordagens fornecem exatamente as mesmas estimativas pontuais:

	MQ	MV	Bayes (moda/mediana posterior)
\(\beta_0\)	\(b_0 = \bar{y} - b_1\bar{x}\)	\(b_0 = \bar{y} - b_1\bar{x}\)	\(b_0 = \bar{y} - b_1\bar{x}\)
\(\beta_1\)	\(b_1 = S_{xy}/S_{xx}\)	\(b_1 = S_{xy}/S_{xx}\)	\(b_1 = S_{xy}/S_{xx}\)

A coincidência entre MQ e MV é uma propriedade do modelo normal: para \(\sigma^2\) fixo, maximizar a verossimilhança equivale a minimizar a soma de quadrados dos resíduos. A coincidência com a abordagem bayesiana decorre da priori de Jeffreys ser uniforme nos betas, de modo que a moda da posterior é a moda da verossimilhança. Além disso, como a posterior é uma \(t\) de Student (simétrica), moda, mediana e média coincidem.

4.7.2 Estimativa de \(\sigma^2\)

	MQ	MV	Bayes (média posterior)
Fórmula	\(S_e^2 = \dfrac{SQR}{n-2}\)	\(\hat{\sigma}^2_{MV} = \dfrac{SQR}{n}\)	\(E[\sigma^2 \mid y] = \dfrac{(n-2)S_e^2/2}{(n-2)/2 - 1} = \dfrac{(n-2)S_e^2}{n-4}\)
Viés	Não-viesado	Viesado por fator \((n-2)/n\)	—

O estimador de MQ (\(S_e^2\)) é não-viesado por construção. O de MV (\(\hat{\sigma}^2_{MV} = SQR/n\)) subestima \(\sigma^2\) em média, e a correção usual é multiplicar por \(n/(n-2)\), recuperando \(S_e^2\). Na abordagem bayesiana, a média da Gama-Inversa com forma \(\alpha = (n-2)/2\) e escala \(\lambda = (n-2)S_e^2/2\) é \(\lambda/(\alpha - 1)\), que existe apenas para \(n > 4\).

4.7.3 Variâncias dos estimadores / incerteza

	MQ (exata)	MV (assintótica)	Bayes (posterior)
\(\text{Var}[b_0]\)	\(\dfrac{\sigma^2 \overline{xx}}{n\,S_{xx}}\)	\([I_O^{-1}]_{11} = \dfrac{\sigma^2 \overline{xx}}{n\,S_{xx}}\)	\(\dfrac{g}{g-2} S_{\beta_0}^2 = \dfrac{g}{g-2} \dfrac{S_e^2\,\overline{xx}}{n\,S_{xx}}\)
\(\text{Var}[b_1]\)	\(\dfrac{\sigma^2}{n\,S_{xx}}\)	\([I_O^{-1}]_{22} = \dfrac{\sigma^2}{n\,S_{xx}}\)	\(\dfrac{g}{g-2} S_{\beta_1}^2 = \dfrac{g}{g-2} \dfrac{S_e^2}{n\,S_{xx}}\)
\(\text{Var}[\sigma^2]\)	\(\dfrac{2\sigma^4}{n-2}\)	\(I_O(\sigma^2)^{-1} \approx \dfrac{2\sigma^4}{n}\)	\(\dfrac{2[(n-2)S_e^2/2]^2}{[(n-2)/2-1]^2[(n-2)/2-2]}\)

Na coluna de MQ e MV, \(\sigma^2\) é o valor verdadeiro (desconhecido); na prática, substitui-se por \(S_e^2\). A variância assintótica de MV para \(\sigma^2\) difere ligeiramente da exata (denominador \(n\) ao invés de \(n-2\)), convergindo para o mesmo valor quando \(n\) é grande. Na abordagem bayesiana, a variância posterior de \(\sigma^2\) vem da fórmula da Gama-Inversa com forma \(\alpha = (n-2)/2\) e escala \(\lambda = (n-2)S_e^2/2\).

4.7.4 Distribuições

	MQ / Frequentista	MV (assintótica)	Bayes (posterior)
\(\beta_0\)	\(\dfrac{b_0 - \beta_0}{\text{EP}[b_0]} \sim t_{n-2}\)	\(b_0 \;\dot\sim\; N\!\left(\beta_0,\, [I_O^{-1}]_{11}\right)\)	\(\beta_0 \mid y \sim \text{St}\!\left(n-2,\; b_0,\; S_{\beta_0}\right)\)
\(\beta_1\)	\(\dfrac{b_1 - \beta_1}{\text{EP}[b_1]} \sim t_{n-2}\)	\(b_1 \;\dot\sim\; N\!\left(\beta_1,\, [I_O^{-1}]_{22}\right)\)	\(\beta_1 \mid y \sim \text{St}\!\left(n-2,\; b_1,\; S_{\beta_1}\right)\)
\(\sigma^2\)	\(\dfrac{(n-2)S_e^2}{\sigma^2} \sim \chi^2_{n-2}\)	\(\hat{\sigma}^2_{MV} \;\dot\sim\; N\!\left(\sigma^2,\, I_O(\sigma^2)^{-1}\right)\)	\(\sigma^2 \mid y \sim \text{Gama-Inv}\!\left(\tfrac{n-2}{2},\;\tfrac{(n-2)S_e^2}{2}\right)\)

Algumas observações sobre as distribuições:

• Na abordagem frequentista (MQ), as distribuições são amostrais: a quantidade pivotal \((b_j - \beta_j)/\text{EP}[b_j]\) segue uma \(t\) de Student, onde \(\beta_j\) é o valor fixo e desconhecido e \(b_j\) é aleatório.
• Na abordagem de MV, as distribuições são assintóticas (indicadas pelo \(\dot\sim\)): normais centradas no valor verdadeiro. Para \(n\) grande, a normal assintótica aproxima a \(t\) com \(n-2\) graus de liberdade. Neste modelo particular, as variâncias assintóticas coincidem com as exatas para os betas (pois \(I_E = I_O\)), mas diferem ligeiramente para \(\sigma^2\).
• Na abordagem bayesiana, as distribuições são posteriores: descrevem a incerteza sobre o parâmetro dado os dados observados. O parâmetro é aleatório e os dados são fixos — a inversão conceitual em relação à abordagem frequentista.

4.7.5 Intervalos

	Frequentista (95%)	Bayesiano (ICr 95%)
\(\beta_j\)	\(b_j \pm t_{0.975,\,n-2} \cdot \text{EP}[b_j]\)	\(b_j \pm t_{0.975,\,n-2} \cdot S_{\beta_j}\)
\(\sigma^2\)	\(\left[\dfrac{(n-2)S_e^2}{\chi^2_{0.975}},\;\dfrac{(n-2)S_e^2}{\chi^2_{0.025}}\right]\)	Percentis da \(\text{Gama-Inv}\!\left(\tfrac{n-2}{2},\;\tfrac{(n-2)S_e^2}{2}\right)\)

Para os betas, os intervalos de confiança (frequentistas) e de credibilidade (bayesianos) são numericamente idênticos sob a priori de Jeffreys — embora suas interpretações sejam distintas. O intervalo frequentista afirma que, em repetidas amostras, 95% dos intervalos construídos dessa forma conteriam o verdadeiro \(\beta_j\). O intervalo bayesiano afirma que, dados os dados observados, a probabilidade posterior de \(\beta_j\) estar no intervalo é 0.95.

Para \(\sigma^2\), o intervalo frequentista é baseado na distribuição qui-quadrado, enquanto o bayesiano usa os percentis da Gama-Inversa. Numericamente, ambos são muito próximos para \(n\) moderado, mas não são idênticos.

4.7.6 Síntese

A equivalência numérica entre as três abordagens neste modelo é uma consequência de três fatores: (i) a priori de Jeffreys é não-informativa, (ii) a verossimilhança normal domina a priori, e (iii) existe solução analítica fechada. Em modelos mais complexos — sem solução fechada, com prioris informativas, ou com verossimilhanças não-normais — as abordagens tipicamente divergem, e cada uma oferece vantagens distintas.