1. Completando o quadrado

O procedimento de completar o quadrado consiste em reescrever uma expressão quadrática, como a soma de um quadrado perfeito mais uma constante.

  1. Mostre que, para \(a > 0\), \[ax^2 - 2bx = a\!\left(x - \frac{b}{a}\right)^{\!2} - \frac{b^2}{a}.\]

  2. Deduza a identidade \[\sum_{i=1}^n (y_i - \theta)^2 = \sum_{i=1}^n(y_i - \bar{y})^2 + n(\bar{y} - \theta)^2.\]

Soluções:

  1. \(a(x - b/a)^2 - b^2/a = ax^2 - 2bx + b^2/a - b^2/a = ax^2 - 2bx\).

  2. \(\sum_i(y_i-\theta)^2 = \sum_i\big[(y_i-\bar{y})+(\bar{y}-\theta)\big]^2 = \sum_i(y_i-\bar{y})^2 + 2(\bar{y}-\theta)\underbrace{\sum_i(y_i-\bar{y})}_{=\,0} + n(\bar{y}-\theta)^2\).

2. Resolvendo integrais: o truque da constante normalizadora

Muitas integrais que aparecerão ao longo do curso podem ser resolvidas sem nenhuma técnica de integração, apenas observando que o integrando é o núcleo de uma densidade de probabilidade conhecida, ou seja o integrando é, a menos de uma constante, a expressão de uma função de densidade de probabilidades.

Suponha que, ao olhar para um integrando \(g(x)\), você reconheça que existe uma constante \(\gamma\) tal que \(\gamma \cdot g(x)\) é uma densidade de probabilidade válida, isto é,

\[\int \gamma \, g(x)\, dx = 1.\]

Então, para calcular \(\int g(x)\,dx\), basta multiplicar e dividir o integrando por \(\gamma\) e colocar a constante \(\gamma^{-1}\) para fora da integral:

\[\int g(x)\, dx = \int \gamma^{-1}\big[\gamma\, g(x)\big]\, dx = \gamma^{-1}\int \gamma\, g(x)\, dx = \gamma^{-1}\cdot 1 = \gamma^{-1}.\]

Note-se que \(\gamma\) é uma constante para a integral, ou seja, pode depender de outros parâmetros, mas não da variável de integração.

Em resumo: a integral é igual ao inverso da constante que falta para que o integrando seja uma função de densidade de probabilidades.

Exemplo 1:

Para resolver a integral \(\displaystyle\int_0^{+\infty} {\rm e}^{-0,25x} \, dx\) notamos que o integrando é o núcleo de uma distribuição de probabilidades exponencial de parâmetro 0,25. Portanto \[ \displaystyle\int_0^{+\infty} e^{-0{,}25x} \, dx = \displaystyle\int_0^{+\infty} \frac{0{,}25}{0{,}25}\,e^{-0{,}25x} \, dx = \frac{1}{0{,}25}\underbrace{\displaystyle\int_0^{+\infty} 0{,}25\, e^{-0{,}25x} \, dx}_{=\,1} = \frac{1}{0{,}25} \cdot 1 = 4. \]

Em cada questão a seguir, encontre essa constante \(\gamma\) e resolva a integral.

  1. \(\displaystyle\int_0^1 x^3 (1-x)^2\, dx\)

  2. \(\displaystyle\int_0^\infty \lambda^{4}\, e^{-3\lambda}\, d\lambda\)

  3. \(\displaystyle\int_0^\infty x^{a-1}\, e^{-x/\theta}\, dx\), para \(a>0\) e \(\theta>0\)

  4. \(\displaystyle\int_{-\infty}^{\infty} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) dx\)

  5. \(\displaystyle\int_{-\infty}^{\infty} \exp\left(-\frac{1}{2}\big(n x^2 - 2n\bar{x}\, x\big)\right) dx\)

  6. \(\displaystyle\int_{-\infty}^{\infty} \exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n (x_i-\mu)^2\right) d\mu\)

  7. \(\displaystyle\int_0^\infty \lambda^{\sum_i x_i}\, e^{-n\lambda}\cdot \lambda^{a-1} e^{-b\lambda}\, d\lambda\)

  8. \(\displaystyle\int_0^1 p^{\sum_i x_i + a - 1} (1-p)^{n - \sum_i x_i + b - 1}\, dp\)

  9. \(\displaystyle\int_{\mathbb{R}^p} \exp\left(-\frac12(\boldsymbol\beta-\mathbf{b})^\top \Sigma^{-1} (\boldsymbol\beta-\mathbf{b})\right) d\boldsymbol\beta\)

  10. \(\displaystyle\int_0^\infty x^{\alpha-1} (1+x)^{-(\alpha+\beta)}\, dx\)

3. Função de verossimilhança

A função verossimilhança responde a uma pergunta bem específica: dado que observamos estes dados, quão plausível é cada possível valor do parâmetro? A função de verossimilhança é construída a partir do modelo probabilístico usado para gerar os dados, porém, em vez de perguntar qual a probabilidade dos dados para um determinado valor do parâmetro, perguntamos como essa probabilidade varia conforme mudamos o parâmetro, mantendo os dados fixos.

Em quase todos os casos, os passos são os mesmos:

O exemplo a seguir ilustra esses passos em detalhe.

Exemplo 2:

Sejam \(Y_1,\dots,Y_n \mid \beta \overset{iid}{\sim} \text{Gama}(\alpha,\beta)\), com \(\alpha\) conhecido (forma) e \(\beta\) desconhecido (taxa), com densidade

\[f(y;\alpha,\beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}\, y^{\alpha-1} e^{-\beta y}, \qquad y>0.\]

Para determinar a verossimilhança neste caso, seguimos os seguintes passos:

Passo 1: Verossimilhança.

Como \(Y_1,\dots,Y_n\) são independentes, a densidade conjunta da amostra é o produto das densidades marginais:

\[f(\mathbf y\mid\beta) = \prod_{i=1}^n \frac{\beta^\alpha}{\Gamma(\alpha)}\, y_i^{\alpha-1} e^{-\beta y_i}.\]

A função de verossimilhança é essa mesma expressão, mas encarada como função de \(\beta\), com os dados \(\mathbf y\) fixos (observados):

\[L(\beta;\mathbf y) \equiv f(\mathbf y\mid\beta).\]

Passo 2: Separar o que depende de \(\beta\) do que não depende.

Separando o produtório em três partes:

\[L(\beta;\mathbf y) = \beta^{n \alpha} \left(\frac{1}{\Gamma(\alpha)}\right)^n \cdot \left(\prod_{i=1}^n y_i^{\alpha-1}\right) \cdot \exp\left(-\beta\sum_{i=1}^n y_i\right).\]

Passo 3: Simplificar.

Os termos \(\prod_i y_i^{\alpha-1}\) e \(\Gamma(\alpha)^{n}\) não dependem de \(\beta\), portanto podem ser absorvidos em uma constante \(k(\mathbf y)\) que não interfere na forma da verossimilhança como função de \(\beta\):

\[L(\beta;\mathbf y) = \beta^{n\alpha}\, k(\mathbf y)\, e^{-\beta\sum_i y_i}.\]

Passo 4: Forma final, como função de \(\beta\).

\[\boxed{L(\beta;\mathbf y) \propto \beta^{\,n\alpha}\, e^{-\beta\sum_i y_i}}\]

Note que essa expressão tendo \(\beta\) como variavel, tem exatamente a forma de um núcleo de densidade Gama\((n\alpha+1, \sum_i y_i)\).

Ilustração numérica e gráfica.

A log-verossimilhança (a menos de constante), \[\boxed{\ell(\beta;\mathbf y) = n\alpha\log\beta - \beta\sum_i y_i}\] é côncava em \(\beta\) e seu ponto de máximo é obtido derivando e igualando a zero:

\[\frac{d\ell}{d\beta} = \frac{n\alpha}{\beta} - \sum_i y_i = 0 \quad\Longrightarrow\quad \hat\beta_{MV} = \frac{n\alpha}{\sum_i y_i}.\]

Para visualizar a função de verossimilhança, simulamos uma amostra de tamanho \(n=15\) de uma \(\text{Gama}(\alpha=3, \beta=2)\) e calculamos \(L(\beta;\mathbf y)\) para vários valores de \(\beta\).

set.seed(42)
alphaF <- 3
betaT <- 2
n <- 15
y <- rgamma(n, shape = alphaF, rate = betaT)
##
sy <- sum(y)
n_alpha <- n * alphaF
betaMV <- n_alpha/sy
##
cat("Amostra simulada (arredondada a 2 casas decimais):\n")
## Amostra simulada (arredondada a 2 casas decimais):
round(y, 2)
##  [1] 2.57 0.84 1.29 0.53 1.17 0.78 3.35 1.20 0.54 0.39 2.51 2.14 2.04 1.66 2.13
cat("\nn =", n, "  soma(y_i) =", round(sy, 3),
    "  n*alpha =", n_alpha,
    "\nEstimador de máxima verossimilhança: beta_hat =", round(betaMV, 4))
## 
## n = 15   soma(y_i) = 23.152   n*alpha = 45 
## Estimador de máxima verossimilhança: beta_hat = 1.9436

Os gráficos abaixo mostram a verossimilhança \(L(\beta;\mathbf y)\) e a log-verossimilhança \(\ell(\beta;\mathbf y)\) (normalizadas para terem máximo igual a 1 e a 0, respectivamente, apenas para facilitar a visualização), com a linha vertical tracejada indicando \(\hat\beta_{MV}\).

beta.grid <- seq(0.5, 4.5, length.out = 101)

logL <- n_alpha * log(beta.grid) - beta.grid * sy
logLrel <- logL - max(logL)          
Lrel <- exp(logLrel)                 

op <- par(mfrow = c(1,2), mar = c(4,4,3,1), mgp = c(2,1,0))

plot(beta.grid, Lrel, type = "l", lwd = 2, col = "blue",
     xlab = expression(beta),
     ylab = expression(L(beta~";"~y) / L(hat(beta)[MV]~";"~y)),
     main = "Verossimilhança relativa")
abline(v = betaMV, lty = 2, col = "red")

plot(beta.grid, logLrel, type = "l", lwd = 2, col = "blue",
     xlab = expression(beta),
     ylab = expression(l(beta~";"~y) - l(hat(beta)[MV]~";"~y)),
     main = "Log-verossimilhança relativa")
abline(v = betaMV, lty = 2, col = "red")

par(op)

Repare que a curva é unimodal, com o pico exatamente em \(\hat\beta_{MV}=n\alpha/\sum_i y_i\), próximo do valor verdadeiro \(\beta=2\) usado na simulação.

Verossimilhanças de alguns modelos:

Para cada modelo abaixo, escreva a expressão da função de verossimilhança, simplificando o que for possível.

  1. \(Y_1,\dots,Y_n \mid p \overset{iid}{\sim} \text{Bernoulli}(p)\).

  2. \(Y_1,\dots,Y_n \mid \lambda \overset{iid}{\sim} \text{Poisson}(\lambda)\).

  3. \(Y_1,\dots,Y_n \mid \theta \overset{iid}{\sim} \text{Exponencial}(\theta)\) (taxa \(\theta\)).

  4. \(Y_1,\dots,Y_n \mid \mu \overset{iid}{\sim} N(\mu,\sigma_0^2)\), com \(\sigma_0^2\) conhecido.

  5. \(Y_1,\dots,Y_n \mid \mu,\sigma^2 \overset{iid}{\sim} N(\mu,\sigma^2)\), ambos desconhecidos.

  6. Uma única observação \(Y\mid n,\theta \sim \text{Binomial}(n,\theta)\), com \(n\) conhecido.

  7. \(Y_1,\dots,Y_n \mid \theta \overset{iid}{\sim} \text{Uniforme}(0,\theta)\).

  8. \(Y_1,\dots,Y_n \mid \theta \overset{iid}{\sim} \text{Geométrica}(\theta)\), com suporte em \(\{1,2,3,\dots\}\) e função de probabilidade \(P(Y=y)=\theta(1-\theta)^{y-1}\).

  9. Uma moeda é lançada repetidamente até que a primeira cara ocorra, e o número de lançamentos necessários é registrado. Esse experimento é repetido de forma independente \(n\) vezes, gerando \(y_1,\dots,y_n\). Proponha um modelo probabilístico para os dados e escreva a verossimilhança do parâmetro de interesse.

  10. O número de acessos a um site é registrado a cada minuto, ao longo de \(n\) minutos, resultando em \(y_1,\dots,y_n\) acessos por minuto. Proponha um modelo para os dados e escreva a verossimilhança do parâmetro de interesse.

  11. Uma empresa quer estimar a proporção de peças defeituosas em sua linha de produção. Para isso, retira uma amostra de \(n\) peças e conta quantas são defeituosas. Proponha um modelo para o número de peças defeituosas e escreva a verossimilhança do parâmetro de interesse.

  12. O tempo de vida (em horas) de \(n\) componentes eletrônicos idênticos é medido até a queima de cada um. Proponha um modelo para esses tempos e escreva a verossimilhança do parâmetro de interesse.

4. Séries de Taylor e aproximações quadráticas

A série de Taylor permite aproximar uma função suave \(f(\theta)\), na vizinhança de um ponto \(\theta_0\), por um polinômio construído a partir de suas derivadas em \(\theta_0\):

\[f(\theta) \approx f(\theta_0) + f'(\theta_0)(\theta-\theta_0) + \frac{1}{2}f''(\theta_0)(\theta-\theta_0)^2 + \dots\]

Truncando na segunda ordem, obtemos a aproximação quadrática, a mais usada em inferência bayesiana:

\[f(\theta) \approx f(\theta_0) + f'(\theta_0)(\theta-\theta_0) + \frac{1}{2}f''(\theta_0)(\theta-\theta_0)^2.\]

No curso, essa ideia será aplicada à log-verossimilhança \(\ell(\theta) = \log L(\theta;\mathbf y)\), expandida em torno do estimador de máxima verossimilhança \(\hat\theta\). Como \(\hat\theta\) é ponto de máximo, \(\ell'(\hat\theta)=0\), e o termo linear desaparece, restando

\[\ell(\theta) \approx \ell(\hat\theta) - \frac{1}{2} I(\hat\theta)\,(\theta-\hat\theta)^2, \qquad I(\hat\theta) = -\ell''(\hat\theta),\]

em que \(I(\hat\theta)\) é a informação observada de Fisher. Exponenciando os dois lados, a verossimilhança perto de \(\hat\theta\) se comporta como o núcleo de uma densidade Normal:

\[L(\theta;\mathbf y) \approx L(\hat\theta;\mathbf y)\, \exp\left(-\frac{1}{2}I(\hat\theta)(\theta-\hat\theta)^2\right),\]

o que sugere a aproximação normal

\[\theta \,\dot\sim\, N\!\left(\hat\theta,\ I(\hat\theta)^{-1}\right).\]

É essa aproximação que será usada mais adiante para construir intervalos de credibilidade (ou confiança) aproximados, do tipo \(\hat\theta \pm z_{1-\alpha/2}\sqrt{I(\hat\theta)^{-1}}\).

Exemplo 3:

Sejam \(Y_1,\dots,Y_n \mid \lambda \overset{iid}{\sim} \text{Poisson}(\lambda)\). A log-verossimilhança (a menos de constante) é

\[\ell(\lambda;\mathbf y) = -n\lambda + \Big(\sum_{i=1}^n y_i\Big) \log\lambda.\]

Passo 1: Derivadas.

\[\ell'(\lambda) = -n + \frac{\sum_i y_i}{\lambda}, \qquad \ell''(\lambda) = -\frac{\sum_i y_i}{\lambda^2}.\]

Passo 2: Estimador de máxima verossimilhança. Igualando \(\ell'(\lambda)=0\):

\[\hat\lambda = \frac{\sum_i y_i}{n} = \bar y.\]

Passo 3: Informação observada. Substituindo \(\hat\lambda\) em \(-\ell''(\lambda)\), e usando \(\sum_i y_i = n\hat\lambda\):

\[I(\hat\lambda) = -\ell''(\hat\lambda) = \frac{\sum_i y_i}{\hat\lambda^2} = \frac{n}{\hat\lambda}.\]

Passo 4: Aproximação quadrática e normal resultante.

\[\ell(\lambda;\mathbf y) \approx \ell(\hat\lambda;\mathbf y) - \frac{1}{2}\,\frac{n}{\hat\lambda}\,(\lambda-\hat\lambda)^2 \quad\Longrightarrow\quad \boxed{\lambda \,\dot\sim\, N\!\left(\hat\lambda,\ \frac{\hat\lambda}{n} \right)}\]

Essa é a aproximação que fundamenta o intervalo de confiança usual para \(\lambda\) no modelo Poisson: \(\hat\lambda \pm z_{1-\alpha/2} \sqrt{\hat\lambda/n}\).

Ilustração numérica e gráfica. Simulamos uma amostra de tamanho \(n=20\) de uma \(\text{Poisson}(\lambda=4)\), calculamos \(\hat\lambda\) e \(I(\hat\lambda)\), e comparamos a verossimilhança exata (que, como função de \(\lambda\), é o núcleo de uma densidade Gama (veja a questão (g) do exercício 2 ou a questão (b) do exercício 3) com a aproximação quadrática (normal) obtida pela série de Taylor.

set.seed(123)
lambdaT <- 4
n <- 20
y <- rpois(n, lambda = lambdaT)
##
sy <- sum(y)
lambdaMV <- mean(y)
IhatMV <- sy / lambdaMV^2   # = n / lambdaMV
##
cat("Amostra simulada:\n")
## Amostra simulada:
y
##  [1] 3 6 3 6 7 1 4 7 4 4 8 4 5 4 2 7 3 1 3 8
cat("\nn =", n, "  soma(y_i) =", sy,
    "  lambda_hat (EMV) =", round(lambdaMV, 4),
    "\nInformação observada I(lambda_hat) =", round(IhatMV, 4))
## 
## n = 20   soma(y_i) = 90   lambda_hat (EMV) = 4.5 
## Informação observada I(lambda_hat) = 4.4444
lambda.grid <- seq(1.5, 7.5, length.out = 200)

logL_exact <- -n * lambda.grid + sy * log(lambda.grid)
logL_exact_rel <- logL_exact - max(logL_exact)
L_exact_rel <- exp(logL_exact_rel)

logL_quad <- -0.5 * IhatMV * (lambda.grid - lambdaMV)^2
L_quad_rel <- exp(logL_quad)

op <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1), mgp = c(2, 1, 0))

plot(lambda.grid, L_exact_rel, type = "l", lwd = 2, col = "blue",
     xlab = expression(lambda), ylab = "Verossimilhança relativa",
     main = "Exata vs. aproximação normal")
lines(lambda.grid, L_quad_rel, lwd = 2, lty = 2, col = "darkorange")
abline(v = lambdaMV, lty = 3, col = "gray40")
legend("topright",
       legend = c("L exata (núcleo Gama)", "Aprox. quadrática (Normal)"),
       col = c("blue", "darkorange"), lty = c(1, 2), lwd = 2,
       bty = "n", cex = 0.8)

plot(lambda.grid, logL_exact_rel, type = "l", lwd = 2, col = "blue",
     xlab = expression(lambda), ylab = "Log-verossimilhança relativa",
     main = "Log-verossimilhança: exata vs. quadrática")
lines(lambda.grid, logL_quad, lwd = 2, lty = 2, col = "darkorange")
abline(v = lambdaMV, lty = 3, col = "gray40")

par(op)

Perto de \(\hat\lambda\), as duas curvas praticamente coincidem. É exatamente isto o que a série de Taylor garante localmente. Mas repare que a curva exata é assimétrica (mais pesada à direita, herdando o formato da densidade Gama), enquanto a aproximação normal é, por construção, simétrica em torno de \(\hat\lambda\). Para \(n\) pequeno, isso pode fazer a aproximação normal subestimar a plausibilidade de valores grandes de \(\lambda\) (um ponto importante para ter em mente ao usar essa aproximação para construir intervalos).

Exercícios:

  1. Escreva a expansão de Taylor de segunda ordem de uma função genérica \(f(\theta)\) em torno de um ponto \(\theta_0\), e explique por que, quando \(\theta_0=\hat\theta\) é um ponto de máximo de \(f(\theta)\), o termo linear da expansão desaparece.

  2. Sejam \(Y_1,\dots,Y_n \mid p \overset{iid}{\sim} \text{Bernoulli}(p)\). Escreva \(\ell(p)\), calcule \(\ell'(p)\) e \(\ell''(p)\), encontre o EMV \(\hat p\) e a informação observada \(I(\hat p)\), e escreva a aproximação normal resultante para \(p\).

  3. Sejam \(Y_1,\dots,Y_n \mid \theta \overset{iid}{\sim} \text{Exponencial}(\theta)\) (taxa \(\theta\)). Repita o procedimento do Exemplo 3 para obter o EMV \(\hat\theta\), a informação observada \(I(\hat\theta)\), e a aproximação normal resultante para \(\theta\).

  4. Sejam \(Y_1,\dots,Y_n \mid \mu \overset{iid}{\sim} N(\mu,\sigma_0^2)\), com \(\sigma_0^2\) conhecido. Mostre que, neste caso, \(\ell''(\mu)\) não depende de \(\mu\) e, portanto, a expansão de Taylor de segunda ordem de \(\ell(\mu)\) em torno de \(\hat\mu=\bar y\) é exata (sem termos de ordem superior sendo descartados). Confira que a “aproximação” normal resultante coincide com a distribuição amostral exata de \(\bar y\).