4.1 População e Amostra.

Primeiro, gostaria de lembrar dos conceitos de população e amostra. Enquanto a população representa completamente um conjunto, a amostra representa parte desse conjunto. Por exemplo, em uma sala de aula há 30 alunos. A população são os alunos nessa sala de aula. Caso eu selecione aleatoriamente 15 alunos nessa sala de aula, esses alunos que eu selecionei representam uma amostra da população total. Mas por que isso é relevante?

No caso de populações pequenas, como o caso do exemplo acima, uma amostra não faz muito sentido, pois o custo associado de selecionar todos os 30 alunos é muiro baixo. E no caso de populações muito grandes? Imagine que eu faça uma pesquisa na POPULAÇÃO brasileira inteira. Qual seria o custo de perguntar para cada um dos 211 milhões de brasileiros perguntas relativas a essa minha pesquisa? Enorme! Por isso, uma amostra dessa população reduz gritantemente os custos associados a essa pesquisa, mantendo certo nível de confiabilidade a ela.

Observação: existem diversos métodos de amostragem e não entrarei no mérito delas.

Falei de confiabilidade agora, o que é isso? Bom, quando faço uma amostra, passo a desconhecer certos aspectos de uma população, e isso introduz incerteza na minha pesquisa. A estatística tem absolutamente tudo a ver com a incerteza!

Sabendo muito bem a diferença entre amostra e população, quero introduzir alguns conceitos (não tão) abstratos:

A esperança matemática de uma variável \(E(X_i)\) é uma média dessa variável considerando a população toda, enquanto que a média aritmética \(\mu\) é uma estimativa do valor populacional (verdadeiro) em dada amostra. Exemplo: a média de idade para aquela turma de 30 alunos é de 20 anos, assim, \(E(X_i) = 20\), enquanto que a média para a amostra de 15 alunos é 22, sendo \(\mu = 22\).

A variância de uma variável é uma representação de como cada ponto observado se distancia de sua média. Novamente, existem diferenças entre os valores populacionais e amostrais, mas na prática, representam fenômenos similares. Tomamos \(Var(X_i)\) como uma notação populacional para a variância e \(\sigma^2\) como a notação da variância amostral. Ufa! Eu vou mostrar isso abaixo no exemplo que dei para sala de aula.

# Gero os dados para sala de aula
dataset <- data.table(aluno = 1:30, idade = sample(18:24, 10, replace = TRUE))

dataset %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

Calculando a média das idades nessa sala, encontro a esperança da idade dessa população de alunos, que é 21 anos, enquanto que a variância populacional é de 3.1034483! Se eu selecionar sortear 15 alunos nessa população (amostra simples), teria essa seguinte tabela:

dataset_sorteado <- slice_sample(dataset, n = 15)

dataset_sorteado %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

Calculando a média das idades nessa amostra, encontro a estimativa da média populacional dos alunos, que é 20.7333333 anos, enquanto que a variância da amostra é de 2.7809524!

Note que os valores de média são diferentes, apesar de se tratar de uma mesma população.

Para finalizar esses conceitos, vou apresentar um novo conceito, a média das médias. Matematicamente eu me refiro a esse princípio:

\[ E(X_i) = \frac{1}{n}\sum_{i \rightarrow n}^n \mu_i \] A lei dos grandes números é um princípio estatístico muito interessante que gostaria de relembrar. Quanto mais amostras são realizadas em uma população, a média das médias dessas amostras tende “assintoticamente” para o valor populacional (caso as amostras tenham boa representatividade da população). Vou fazer esse mesmo sorteio de 15 alunos e calcular a média. Na prática, posso fazer o seguinte:

dataset_sorteado <- data.table(replicate(15, sample(dataset$idade, 15)))
for (i in 1:15) {
  setnames(dataset_sorteado, paste0("V", i), paste0("Sorteio ", i))
}

dataset_sorteado %>%
  kbl(align = "cc") %>%
  kable_classic() %>%
  scroll_box(width = "800px", height = "350px")

A média de cada uma dessas colunas é:

medias <- transpose(data.table(apply(dataset_sorteado, 1, mean)))
for (i in 1:15) {
  setnames(medias, paste0("V", i), paste0("Média da amostra ", i))
}

medias %>%
  kbl(align = "cc", digits = 3) %>%
  kable_classic() %>%
  scroll_box(width = "800px", height = "300px")

Se eu somar essas médias amostrais, encontro o valor de 21.0133333 que, comparado com a média populacional de 21 é extremamente próximo.

Voilà! Definimos os conceitos de media das médias, esperança populacional e média da amostra. Esses conceitos são relevantes pois podemos realizar inferências na população com base na nossa amostra. A variância populacional pode ser definida com base na esperança das observações.

Por curiosidade, a \(Var(X_i)\) é entendida como: \[ Var(X_i) = E[(X_i - E(X_i))^2] \]

4.2 As distribuições: o que são e para que servem?

dataset <- data.table(idade = rnorm(20000, 21, 1))

ggplot(dataset, aes(idade)) +
  geom_histogram(fill = "steelblue") +
  theme_classic()

ggplot(dataset, aes(idade)) +
  geom_density(color = "steelblue") +
  theme_classic()

dataset <- data.table(idade = sn::rsn(n = 20000, 18, 3.5, 6))

ggplot(dataset, aes(idade)) +
  geom_density(color = "steelblue") +
  theme_classic() +
  scale_x_continuous(labels = scales::number_format())

dataset <- data.table(face_tirada = sample(1:6, 10, replace = TRUE))

ggplot(dataset, aes(face_tirada)) +
  geom_bar(fill = "steelblue") +
  theme_classic()

dataset <- data.table(face_tirada = sample(1:6, 10000, replace = TRUE))

ggplot(dataset, aes(face_tirada)) +
  geom_bar(fill = "steelblue") +
  theme_classic()

dataset <- data.table(face_tirada = sample(1:6, 10000, replace = TRUE, prob = c(0.5/6, 0.5/6, 0.5/6, 1/6, 1/6, 1.5/6)))

ggplot(dataset, aes(face_tirada)) +
  geom_bar(fill = "steelblue") +
  theme_classic()

dataset <- data.table(x = rnorm(100000))

ggplot(dataset, aes(x)) +
  geom_density(color = "steelblue") +
  theme_classic()

dataset[, fda := pnorm(x)]

ggplot(dataset, aes(x, fda)) +
  geom_line(color = "steelblue") +
  theme_classic()

4.3 Testes de hipótese.

Antes de começar, recomendo assistir esse vídeo para uma compreensão sobre como funcionam os testes de uma forma descontraída. Por mais que os testes realizados não se concentram no campo dos testes bayesianos, a lógica de aplicação é a mesma (e o conteúdo vale a pena).

Mas então, o objetivo central dos testes de hipótese advém da área inferencial da estatística. Isto é, da área estatística em que se busca inferir sob uma população quando conhecemos apenas uma amostra dela. No processo de teste, utilizamos níveis de significância para determinar se os resultados deste teste específico são significantes ou não. Mas o que é significância?

Imagine que realizamos um experimento qualquer em que testamos a ocorrência de um evento. A cada 95 de 100 vezes o evento esperado ocorre, e apenas 5 vezes esse evento não ocorre. Essas 5 vezes que o evento não ocorreu é fruto do acaso, sendo a significância do acaso de 5%. Nesse sentido, significância é o quanto um experimento atribui a não ocorrência do evento esperado ao acaso. Por isso, optar por um nível maior de significância significa reduzir a confiança do experimento, visto que o acaso passa a ter maior impacto.

O que a significância representa no caso das distribuições? Lembrando que a área abaixo de uma fdp representa a probabilidade de ocorrência de um evento, os valores extremos dessa distribuição representam o acaso! Por exemplo, em que intervalo se concentram 95% das informações mais frequentes na distribuição normal abaixo?

dataset <- data.table(x = rnorm(100000))

ggplot(dataset, aes(x)) +
  geom_density(color = "steelblue") +
  theme_classic()

Pela fdp fica difícil de ver isso né? E pela fda? (Vou facilitar ainda mais colocando linhas!)

dataset[, fda := pnorm(x)]

ggplot(dataset, aes(x, fda)) +
  geom_line(color = "steelblue") +
  geom_hline(yintercept = 0.025, color = "red") +
  geom_hline(yintercept = 0.975, color = "red") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic() +
  labs(title = "Extremos de um teste bicaudal de uma distribuição normal clássica.")

Se você entendeu bem o conceito de fda, pode ter compreendido que 95% dos dados estão contidos entre as linhas vermelhas. Mas por que essas linhas estão saindo de 0.025 e 0.975 você deve estar se perguntando. A resposta é simples, os valores mais extremos possíveis dessa distribuição são suas caudas. Testar a significância de um valor implica em testar se esse valor está dentro do intervalo do acaso!

Com base no que vimos da fda, sabemos então que 95% dos valores estão contidos entre aproximadamente -2 e 2. De volta para a fdp:

ggplot(dataset, aes(x)) +
  geom_density(color = "steelblue") +
  geom_vline(xintercept = -2, color = "red") +
  geom_vline(xintercept = 2, color = "red") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic()

Talvez agora eu esteja refrescando sua memória (caso você tenha aprendido bem sobre testes de hipótese!). De \(-\infty\) a aproximadamente -2 e de 2 até \(\infty\), tenho intervalos atribuídos ao acaso, considerando um nível de significância (normalmente referido a \(\alpha\)) igual a 5%.

A partir dessa consideração, eu já te garanto que não tem muita abstração! Testes de hipótese vão considerar exatamente a distribuição de um evento observado e sua relação com o acaso. Vamos para um exemplo numérico!

Retornando ao exemplo da sala de aula: foram realizados 2 testes de conhecimentos gerais em uma turma de 100 alunos, um no início do semestre letivo, e outro no fim do semestre letivo. Considerando que cada aluno foi avaliado de 0 a 100, gostaria de testar se os conhecimentos gerais dos alunos é o mesmo que no início do ano. Na prática, isso é um teste de diferença de médias! A tabela abaixo representam as notas de cada aluno:

dataset <- data.table(aluno = 1:100, primeira_prova = rnorm(100, 60, 10), segunda_prova = rnorm(100, 65, 10))

dataset %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

A distribuição das notas dos alunos pode ser representada como segue:

dataset_long <- melt(dataset, id.vars = "aluno", variable.name = "prova", value.name = "nota")

ggplot(dataset_long, aes(x = nota, color = prova)) +
  geom_density() +
  theme_classic()

Nota, o desvio padrão das duas distribuições é o mesmo e igual a 10!

Visualmente, a média das notas parece maior na segunda aplicação dessa prova. Posso dar validade a essa observação do experimento por meio do teste Z! O teste é o seguinte:

\[ Z = \frac{\bar{X}_1 - \bar{X}_2}{\frac{\sigma}{\sqrt{n}}} \] Abaixo apresento o maior princípio dos testes de hipótese: as hipóteses! Os testes de hipótese testam uma hipótese nula contra uma hipótese alternativa, sendo o critério de aceitação atrelado à significância! No caso do teste Z, temos as seguintes hipóteses:

\[ H_0: \bar{X_1} = \bar{X_2} \] \[ H_a: \bar{X}_1 \neq \bar{X}_2 \] Sendo \(\bar{X}_1\) e \(\bar{X}_2\) as médias de nota da primeira e segunda aplicação, respectivamente.

Substituo abaixo os valores na equação para encontrar o valor de Z.

\[ Z = \frac{60 - 65}{\frac{10}{\sqrt{100}}} = \frac{-5}{\frac{10}{10}} = -5 \] O que esse valor quer me dizer? Bom, esse valor de Z que calculamos segue a distribuição normal clássica que vimos lá em cima! Isso significa dizer que nosso Z calculado, em módulo, se situa em 2 pontos distintos de nossa distribuição, como segue:

dataset <- data.table(x = rnorm(100000))

ggplot(dataset, aes(x)) +
  geom_density(color = "steelblue") +
  geom_vline(xintercept = -2, color = "red") +
  geom_vline(xintercept = 2, color = "red") +
    geom_vline(xintercept = -5, color = "gray23") +
    geom_vline(xintercept = 5, color = "gray23") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic()

O intervalo de -2 e 2 foram definidos como sendo os pontos que definem o intervalo de 95% das observações, lembra? O que significa que os outros pontos são resultantes do acaso. Talvez você já tenha visto isso em algum lugar, mas 2 é nosso valor Z tabelado! Como nosso valor z calculado é superior ao tabelado, podemos afirmar que nosso teste se encontra no intervalo mais extremo ainda que o acaso, tornando nossos resultados significantes. Em outras palavras, pelo Z calculado ser maior que o Z tabelado, rejeitamos a hipótese nula que as médias são iguais!

Essa interpretação por meio do Z calculado pode ser cansativa por um motivo simples: preciso checar qual o valor de Z tabelado para poder comparar! Seria tão bom se eu tivesse uma forma de olhar somente para um valor e concluir se eu rejeito ou não a hipótese nula…

P-valor vem ao resgate! E se a gente calculasse a probabilidade no intervalo entre meu Z calculado em módulo? Voltando para a fda da distribuição normal:

dataset[, fda := pnorm(x)]

ggplot(dataset, aes(x, fda)) +
  geom_line(color = "steelblue") +
  geom_vline(xintercept = -5, color = "red") +
  geom_vline(xintercept = 5, color = "red") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic() +
  labs(title = "Pvalor de um valor Z em um teste bicaudal.")

Qual a probabilidade de encontrarmos valores menores ou iguais a estatística Z calculada? Essa é a definição do pvalor do nosso teste! Como é um teste bicaudal (na próxima seção explicarei isso melhor), devemos considerar a fda de Z nos dois extremos.

Lembre-se que o objetivo da fda é representar \(P(X \leq x_i)\), ou seja, a probabilidade de encontramos valores de X maiores que \(x_i\). A compreensão da parte esquerda do gráfico faz sentido quando queremos encontrar valores mais extremos do que o nosso Z calculado, e a parte direita?

Não faz! O pvalor considera a chance de encontrarmos valores mais extremos que nossa estatística de teste. Portanto, conceitualmente, para um teste bicaudal, buscamos \(P(X \leq -Z)\) e \(P(X \geq Z)\)! A segunda expressão seria alcançável matematicamente fazendo 1 - fda.

dataset[, fda := pnorm(x)]

ggplot(dataset, aes(x, 1-fda)) +
  geom_line(color = "steelblue") +
  geom_vline(xintercept = -5, color = "red") +
  geom_vline(xintercept = 5, color = "red") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic() +
  labs(title = "P(X >= xi)")

Pela teoria de probabilidades, a probabilidade conjunta de dois valores independentes é a multiplicação dos dois. Nesse caso, o pvalor do teste seria:

\[ pvalor = \int_{-\infty}^{-Z_{calculado}} \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} dx \int_{Z_{calculado}}^{\infty} \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} dx \]

Credo (2)! Por mais que essa equação seja cabeluda, não há motivo para assustar. Veja nos dois gráficos acima. O ponto de x = -5 no primeiro gráfico é 0. O ponto de x = 5 no gráfico de baixo é 0. A probabilidade conjunta dos dois é 0! Então nosso pvalor dá um valor próximo a zero!

Mas então qual a relação do pvalor com o nível de significância? Bom, como a probabilidade de encontrarmos valores mais extremos que nosso valor Z calculado é muito próximo de zero, rejeitamos a hipótese nula por ele ser mais raro que o próprio acaso!

dataset <- data.table(x = rnorm(100000))

ggplot(dataset, aes(x)) +
  geom_density(color = "steelblue") +
  theme_classic()

dataset[, fda := pnorm(x)]

ggplot(dataset, aes(x, fda)) +
  geom_line(color = "steelblue") +
  geom_hline(yintercept = 0.95, color = "red") +
  geom_vline(xintercept = 1.64, color = "gray23") +
  scale_x_continuous(n.breaks = 20) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic() +
  labs(title = "Extremo de um teste unicaudal de uma distribuição normal clássica.")

ggplot(dataset, aes(x)) +
  geom_density(color = "steelblue") +
  geom_vline(xintercept = 1.64, color = "red") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic()

ggplot(dataset, aes(x)) +
  geom_density(color = "steelblue") +
  geom_vline(xintercept = 1.64, color = "red") +
  geom_vline(xintercept = 5, color = "black") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic()

ggplot(dataset, aes(x, 1-fda)) +
  geom_line(color = "steelblue") +
  geom_vline(xintercept = 5, color = "black") +
  geom_hline(yintercept = 0.05, color = "red") +
  scale_x_continuous(n.breaks = 20) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic() +
  labs(title = "P(X >= xi)")

4.4 Graus de Liberdade.

Quando existem poucas observações em uma amostra, considerar a distribuição normal é um erro! Por que? Lembre da lei dos grandes números! Quanto mais observações, mais confiança poderemos ter nas nossas inferências simplesmente porque podemos captar mais características de um evento.

E o que são graus de liberdade? Esse conceito está diretamente relacionado com o número de informações independentes menos o número de parâmetros estatísticos a serem avaliados.

Confuso? Vamos de exemplo! Eu tenho a média de uma amostra (média é um parâmetro estatístico). Essa amostra é composta de 5 observações: 2, 5, 3, 6, x e a média dessa amostra é 3,6. Qual o valor de x? Caso você tenha percebido, o valor de x é conhecido a partir do momento que conhecemos a média da amostra, sendo x necessáriamente igual a 2. Te pergunto, x poderia assumir outro valor?

Como x não pode assumir outro valor, ele não é independente!. Logo, existe um valor que nunca será livre nessa amostra. Nesse caso, a amostra apresenta 4 graus de liberdade.

Mas e se eu desconhecer mais um dos valores dessa mesma amostra de 5? Teríamos 2, 5, 3, z, x com média de 3.6. Conseguimos encontrar valores para x e z? Não! Um ficaria em função do outro! Na prática, em uma amostra simples, sempre podemos calcular a média, então os graus de liberdade são reduzidos a simples equação de \(GL = n - 1\)!

Uma distribuição muito utilizada para penalizar poucas observações (menores graus de liberdade) é a t de student. Sua fdp é dada abaixo:

\[ f(x) = \frac{\Gamma(\frac{gl + 1}{2})}{\sqrt{\pi gl \Gamma (\frac{gl}{2})}} (1 + \frac{x^2}{gl})^{-\frac{(gl + 1)}{2}} \] Credo ao infinito! Por mais que seja uma função bem da esquisita, não há motivo para preocupação. Só mostrei ela aí pra deixar claro de onde estão saindo os desenhos do gráfico abaixo! Vamos analisar a fdp de uma distribuição t de student com 5 e 30 graus de liberdade:

dataset <- data.table(referencia = 1:10000, cinco_gl = rt(10000, 5), trinta_gl = rt(10000, 30))

ggplot(dataset) +
  geom_density(aes(cinco_gl), color = "steelblue") +
  geom_density(aes(trinta_gl), color = "indianred1") +
  labs(x = "x") +
  theme_classic()

Em vermelho represento a distribuição com 30 graus de liberdade e, em azul a distribuição com 5. Qual a diferença das duas?

A “saia” da distribuição em azul é maior. Em termos de práticos isso significa que mais valores estarão nas extremidades na distribuição azul do que na distribuição vermelha. Essa diferença é visível pela fda!

dataset[
  ,
  `:=` (
    fda_cinco = pt(cinco_gl, 5),
    fda_trinta = pt(trinta_gl, 30)
  )
]

ggplot(dataset) +
  geom_line(aes(cinco_gl, fda_cinco), color = "steelblue") +
  geom_line(aes(trinta_gl, fda_trinta), color = "indianred1") +
  geom_hline(yintercept = 0.025, color = "gray23") +
  geom_hline(yintercept = 0.975, color = "gray23") +
  scale_x_continuous(n.breaks = 10) +
  scale_y_continuous(n.breaks = 20) +
  theme_classic() +
  theme_classic() +
  labs(x = "x")

Caso você tenha entendido bem a seção de testes de hipótese, sabe o que isso significa! Pela distribuição azul, ao penalizar mais os graus de liberdade, os valores t (pois nesse caso são distribuidos pela fdp de student), deverão ser mais extremos, implicando que serão precisos t calculados maiores para validarmos o nosso teste de hipótese!

Uma curiosidade: quantos mais graus de liberdades temos em uma amostra, mais nossa distribuição tende a distribuição normal. Fiz uma simulação sobre esse efeito aqui.

Com isso, encerramos a revisão teórica de estatística! Espero que tenha sido esclarecedora e tenha facilitado sua compreensão sobre essa ferramenta científica maravilhosa!

5.1 Estimador de Máxima Verossimilhança

Quando falamos de estimação por máxima verossimilhança, na prática, falamos de uma técnica que busca máximizar função de verossimilhança, de forma a encontrar os valores mais “verossímeis” possíveis. Mas o que é função de verossimilhança? O que é verossimilhança? O que estou maximizando? Todas essas são questões relevantes! Para respondê-las, inicialmente utilizarei um exemplo!

Analise o gráfico abaixo:

rm(list = ls())

dataset <- data.table(x = rnorm(10000, 30, 1))
set.seed(30000) #Semente para garantir a reprodutibilidade
dataset_plot <- copy(slice_sample(dataset, n = 100))
dataset_plot[, count := .N, by = .(x = round(x))]

ggplot(dataset_plot, aes(x, count)) +
  geom_point(color = "gray23") +
  theme_classic()

Com vase nesse gráfico, eu poderia supor que a variável x é normalmente distribuído. Mas qual o valor de sua média \(\mu\) e desvio padrão \(\sigma\)? Uma solução seria de gerar várias fdp, com diferentes valores para média e desvio padrão, até que a linha desenhada pela fdp represente esses pontos.

Para simplificar, vou apenas estimar o valor ideal de média, sendo o desvio padrão conhecido e igual a 1.

dataset_plot[
  ,
  `:=` (
    fdp1 = dnorm(x, 29, 1),
    fdp2 = dnorm(x, 30, 1),
    fdp3 = dnorm(x, 31, 1)
  )
]

ggplot(dataset_plot) +
  geom_point(aes(x, count), color = "gray23") +
  geom_line(aes(x, fdp1 * 100), color = "steelblue") +
  geom_line(aes(x, fdp2 * 100), color = "indianred1") +
  geom_line(aes(x, fdp3 * 100), color = "orange") +
  theme_classic()

Com base nas 3 fdps desenhadas no gráfico acima, qual você diria estimar o melhor valor de média possível para nossa distribuição? Se você falou que é a vermelha você está mais do que certo! Na prática, o que você acabou de fazer foi uma estimação por máxima verossimilhança!

É claro que na maioria das vezes estimar os parâmetros que maximizam a verossimilhança não é uma tarefa tão simples quanto chutar (na mão) valores de alguma coisa até acertar. Nesse caso, já conhecíamos o parâmetro de desvio padrão, e se não fosse o caso? Para isso, vou introduzir a notação matemática do problema abaixo:

\[ L(\mu, \sigma^2;x_1,x_2,...,x_n) = \prod_{j=1}^n f(x_j;\mu,\sigma^2) \]

Tá, mas o que eu quero dizer com a equação acima? Ao lado esquerdo da equação temos os parâmetros da fdp normal, dados \(x_1,x_2,...,x_n\) constantes (como já conhecemos os valores, eles são constantes. \(x_1\) será sempre \(x_1\), \(x_2\) será sempre \(x_2\) e por aí em diante). Do lado direito, temos a introdução de um produtório, que é a multiplicação dos valores de densidade de x, dados os parâmentros média \(\mu\) e variância \(\sigma^2\).

Mas por quê um produtório? Tem um vídeo excelente no canal statquest que vai explicar melhor que eu, mas vou viabilizar uma explicação. Veja o gráfico abaixo:

ggplot(dataset_plot) +
  geom_point(aes(x, count), color = "gray23") +
  geom_line(aes(x, fdp2 * 100), color = "indianred1") +
  theme_classic()

Na função de densidade de probabilidade, a densidade (a linha vermelha nesse caso) equivale a verossimilhança de encontrar x na distribuição. A multiplicação de verossimilhanças em cada valor de x, visto que x são observações com verossimilhanças independentes. Em outras palavras, a multiplicação de cada valor de x da nossa fdp é nossa verossimilhança da estimativa dos parâmetros média \(\mu\) e variância \(\sigma^2\). Por isso um produtório.

Então o objetivo é maximizar a função de verossimilhança considerando que podemos mudar os valores de média e variância! No caso que fizemos a maximização no olho, deixa eu provar que a verossimilhança é máxima somente na curva do meio:

dataset_plot_long <- melt(dataset_plot, id.vars = c("x", "count"), variable.name = "linha", value.name = "fdp")

dataset_plot_long[, verossimilhança := prod(fdp), linha]
dataset_plot_long[
  , 
  `:=` (
    fdp = fdp * 100,
    verossimilhança = as.factor(format(verossimilhança, digits = 2))
)]

ggplot(dataset_plot_long, aes(x, fdp, color = verossimilhança)) +
  geom_line() +
  geom_point(aes(x, count), color = "gray23") +
  theme_classic() +
  labs(y = "count")

Note: notação científica de computador tem essa aparência. Quando você vir e-63 é a mesma coisa que \(\cdot10^{-63}\)!

A maior verossimilhança é realmente nessa estimação da linha do meio, visto que tem um valor maior do que as outras!

Agora que vimos a maximização pela forma visual, podemos estimar por meio da matemática! Como um bom aluno de economia, você já deve ter derivado e igualado a zero milhares de vezes ao longo do curso. Derivar e igualar a zero é uma condição de primeira ordem para máximo local de uma função! Ou seja, podemos derivar aquela função de verossimilhança em meio a encontrar o ponto em que a contribuição para o “incremento” da média seja zero.

\[ \frac{\partial L(\mu, \sigma^2;x_1,x_2,...,x_n)}{\partial \mu} = \frac{\prod_{j=1}^n f(x_j;\mu,\sigma^2)}{\partial \mu} \]

Me dá até arrepios ver uma função dessas! Essa é realmente uma derivada de assustar, por isso não faremos ela! O logarítmo tem sido uma ferramenta da humanidade há anos, e a contribuição dele para a matemática, na área de otimização é ainda mais gracioso!

Por que estou dizendo isso? Os pontos de máximo e mínimo de uma função continuam sendo os mesmos, independente se aplicarmos log na função, ou não! Isso quer dizer que eu posso derivar em cima de um somatório, por causa de uma propriedade do logaritmo, e facilita em muito essa conta aí!

A partir de agora, não tenho interesse em maximizar minha função de verossimilhança, e sim, a função do logarítmo da verossimilhança! (note que agora a função passa a se chamar LL)

\[ LL(\mu, \sigma^2;x_1,x_2,...,x_n) = ln(\prod_{j=1}^n f(x_j;\mu,\sigma^2)) = \sum_{j=1}^n ln(f(x_j;\mu,\sigma^2)) \]

Agora a derivada fica mais fácil, veja:

\[ \frac{\partial LL(\mu, \sigma^2;x_1,x_2,...,x_n)}{\partial \mu} = \frac{\sum_{j=1}^n ln(f(x_j;\mu,\sigma^2))}{\partial \mu} \]

E se eu te contar que essa derivada, substituindo aquela parte \(f(x_j;\mu,\sigma^2))\) pela fdp normal, resulta numa condição de maximização da verossimilhança nessa expressão:

\[ \hat{\mu_n} = \frac{1}{n} \sum_{j=1}^n x_j \] Que é simplesmente a média aritmética clássica! Uma prova matemática mais completa pode ser conferida aqui!

O mesmo pode ser feito caso a variância não fosse conhecida:

\[ \frac{\partial LL(\mu, \sigma^2;x_1,x_2,...,x_n)}{\partial \sigma^2} = \frac{\sum_{j=1}^n ln(f(x_j;\mu,\sigma^2))}{\partial \sigma^2} \]

Que também apresenta nosso estimador amostral de variância clássico:

\[ \hat{\sigma}_n^2 = \frac{1}{n} \sum_{j=1}^n (x_j - \hat{\mu_n}) ^ 2 \]

Pode parecer desanimadora a quantidade de conta que fizemos para chegar em conclusões óbvias, mas não desanime! Utilizei a distribuição normal por ser uma introdução ao assunto, mas, no geral, encontrar o maximizador da função de verossimilhança não é tão fácil assim!

Além disso, muitas vezes a estimação por verossimilhança vai incorrer em muita conta na mão e derivadas cansativas. Uma solução para chegar na estimativa ótima é fazer na força bruta computacional! Podemos maximizar essa função por meio de algoritmos que otimizam o logaritmo da verossimilhança, com base nos parâmetros estimados!

Mais para frente, trabalharemos com estimadores que utilizam da estimação por máxima verossimilhança, e, consequentemente, retornaremos a falar sobre esse assunto. Até aqui, achei melhor dar uma introdução a ideia da estimação e as diferentes formas que ela pode ser feita.

5.2 Método dos Momentos Generalizados

O método de momentos generalizados é um método de estimação de parâmetros populacionais. Na prática, isso quer dizer que utilizamos um pressuposto para gerar estimadores consistentes que obedecem esse pressuposto.

Comparado com o método de Máxima Verossimilhança (MLE), o método de Momentos é extremamente simples e chega em resultados iguais ou próximos aos de MLE, com uma única consideração, de que os estimadores apresentam viés, precisando de correção.

Inicialmente, devem ser considerados os momentos populacionais (os valores esperados de algum valor em um evento) como funções do parâmetro de interesse. Assim, a solução para os parâmetro de interesse é uma estimativa desse parâmetro de interesse!

Por ser uma abordagem mais matemática, e menos visual, como foi o MLE, não vou dar muitos exemplos sobre esse método. Na próxima subseção estimei os estimadores de MQO pelo método de momentos e pela máxima verossimilhança. Quebrei ao máximo as contas pra não deixar nada abstrato, então eu sugiro dar uma olhada! (veja a diferença de estimar na mão pelos dois métodos!)

5.3 Minimos Quadrados Ordinários

Atenção! Essa subseção é recheada de conta na mão para chegar nos estimadores de MQO. Assim, meu objetivo aqui é provar que os estimadores do modelo apresentado em GUJARATI podem ser estimados pelo método de Máxima Verossimilhança e pelo método de Momentos Generalizados.

O modelo linear tem a forma funcional: \[ Y_i = \hat{\beta_0} + \hat{\beta_1}X_i + u_i \] Ou seja, queremos estimadores para \(\hat{\beta_0}\) e \(\hat{\beta_1}\) que melhor se ajustem a X e Y. Note que as fórmulas utilizadas servem somente para essa forma funcional, com um intercepto e um parâmetro de inclinação. Outra forma funcional teria outra fórmula para os estimadores. (Por isso a fórmula matricial do MQO é melhor de se utilizar, pois generaliza as contas e facilita a nossa vida.)

O estimador para o intercepto é: \[ \hat{\beta_0} = \frac{1}{n}\sum_{i=1}^n Y_i - \hat{\beta_1}\frac{1}{n}\sum_{i=1}^nX_i = \bar{Y} - \hat{\beta_1}\bar{X} \] E para o parâmetro de inclinação: \[ \hat{\beta_1}= \frac{n\sum_{i=1}^n Y_iX_i - \sum_{i=1}^nX_i\sum_{i=1}^n Y_i}{n\sum_{i=1}^nX_i^2-(\sum_{i=1}^nX_i)^2} \] Mãos a obra!

6.1 Monte Carlo

Monte Carlo é uma pequena cidade situada em Mônaco, muito conhecida por sua cultura dos jogos de azar. Nesse sentido, o método de Monte Carlo, assim como boa parte da estatística, advém de tentativas de melhorar as chances de ganhar (ou entender) a chance em jogos de azar!

No caso, o método de Monte Carlo contribui para as ciências como um todo, por propor uma solução ao problema de inferência em grandes amostras. Isto é, quando existe uma população muito grande, e não existe possibilidade de coletar grandes amostras para inferir precisamente sobre características dessa população.

O método pode ser descrito da seguinte forma:

Suponha que você quer estimar \(\theta\), o valor esperado de uma variável aleatória X: \[ \theta = E(X) \] Esse valor pode representar uma estimativa qualquer de um valor populacional que se deseja estimar. Em meio a exemplificar o método, desejo estimar a proporção de homens com relação ao número de habitantes em uma cidade. Nesse caso, meu \(\theta\) seria a proporção estimada pelo Método de Monte Carlo.

O método de Monte Carlo permite a compreensão de fenômenos estatísticos de uma forma que não precisamos de um pós doutorado em estatística para chegar em uma resposta! Uns dias atrás eu resolvi um problema de cartas de baralho usando simulações de Monte Carlo, que (provavelmente) poderia usar teoria estatística para responder. Não vou entrar no mérito desse problema, mas se quiser saber me procura que eu explico!

Abaixo, utilizarei o método de Monte Carlo de forma bem genérica. Ele não vai exaltar a beleza do método, porque nesse documento busco a simplicidade dos exemplos, mesmo que eles não sejam os melhores possíveis para demonstrar uma técnica.

Tomemos uma cidade de 5000 habitantes como exemplo. Você realizou uma pesquisa e coletou dados de 500 moradores por meio de amostragem simples:

rm(list = ls())

populacao <- data.table(sexo = sample(c("Mulher", "Homem"), 5000, prob = c(0.7, 0.3), replace = TRUE))
set.seed(42)
amostra <- copy(slice_sample(populacao, n = 500))
amostra[, pessoa := 1:500]
setcolorder(amostra, c("pessoa", "sexo"))

amostra %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

A população média de homens é de 0.301% dos moradores. Sua amostra concluiu que existem 0.294% de homens! Por meio de simulações de Monte Carlo, podemos tornar mais precisa essa estimativa, sem adicionar custos de pesquisa adicional!

A simulação dessa estimativa é a quantidade de vezes que selecionamos um homem na amostra sobre a quantidade de vezes que rodarmos a simulação.

monte_carlo <- mean(replicate(100000, sample(amostra$sexo, size = 1)) == "Homem")

No exemplo acima, eu fiz 100.000 simulações de Monte Carlo em valores aleatórios por simulação, e encontrei uma estimativa mais próxima da amostral, 0.29434% de homens. Esse valor é mais próximo do valor populacional de 0.301% do que da amostra inicial de 0.294%.

Mas você deve estar se questionando: pra que tanto esforço para um ganho de precisão tão pequeno? Bom, no caso do problema utilizado, a utilização de métodos mais complexos é um tanto quanto exagerada. Como estou utilizando como exemplo a estimação de um parâmetro populacional “simples” o método de Monte Carlo não tem muito destaque. Porém, em caso de estimações mais complexas, esse método se destaca muito! (o mesmo vale para o próximo método.)

Para ver uma aplicação mais engenhosa do método de Monte Carlo, clique aqui.

6.2 Bootstrap

A estimação pelo método de reamostragem de bootstrap é uma solução criativa e eficiente para gerar uma estimativa dos valores populacionais utilizando poder computacional.

Antes de mais nada, gostaria de relembrar o que é amostragem simples sem reposição.

Imagine que você está realizando um sorteio com vários nomes em um papel. Existem 3 possíveis ganhadores e 10 pessoas estão competindo. Uma mesma pessoa não pode ganhar mais de um prêmio. Nesse caso, você, como sorteador, tiraria o nome da pessoa que ganhou dos nomes no papel, e fim de jogo. Certo? Pois bem, o que você fez foi uma amostragem simples sem reposição. É uma técnica de amostragem simples pois você não atribui nenhum critério na hora de selecionar um nome no papel, e todos os nomes podem ser sorteados. É uma amostragem sem reposição, porque você não voltou com o nome do primeiro ganhador para o montinho de papéis a serem sorteados!

No caso da amostragem simples com reposição, você colocaria de novo o nome no papel, balançaria para não saber qual papel é qual, e selecionaria um papel.

E o que é reamostragem? É gerar uma nova amostra em cima de uma amostra!

Na prática, um exemplo de reamostragem é dado como segue:

Em uma sala de aula com 30 alunos (população), seleciono aleatoriamente 15 (amostra). Os quinze sorteados foram numerados de 1 a 15, como segue abaixo:

rm(list = ls())
dataset <- data.table(aluno = 1:15)

dataset %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

Uma forma de reamostrar essa amostra é sortear novamente esses 15 alunos, 15 vezes, mas, dessa vez, eu posso sortear um aluno mais de uma vez. Abaixo, a nova amostra pelo processo de reamostragem:

set.seed(42)
dataset_sorteado <- slice_sample(dataset, prop = 1, replace = TRUE)[order(aluno)]

dataset_sorteado %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

Agora um aluno pode aparecer na nova amostra mais de uma só vez! Isso introduz em nossa amostra uma espécie de “ruido” aleatório, importante no processo de bootstrap.

Por que isso tudo é relevante? Bom, a técnica de bootstrap abusa de características estatísticas da reamostragem até não querer mais! Para de mostrar a eficiência do método de bootstrap, vou aplicar um caso numérico!

Você quer conhecer a média de idade da cidade Blob, que conta com 5000 moradores. Você não dispõe de recursos para sair perguntando para todos os moradores quantos anos o cidade Blob tem em média. Por conta disso, você decide construir uma amostra de 50 moradores, de forma a inferir sobre a média de idade da cidade. A tabela abaixo são os dados coletados:

populacao <- data.table(idade = sample(0:99, 500, prob = c(rep(.02, 20), rep(.01, 40), rep(.005, 40)), replace = TRUE))
set.seed(412)
amostra <- copy(slice_sample(populacao, n = 50))
amostra[, pessoa := 1:50]
setcolorder(amostra, c("pessoa", "idade"))

amostra %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

A média da sua amostra é de 33.92 anos! A média de idade da população é de 34.624 anos. A sua amostra erra a média de idade por muito pouco. E caso você tenha interesse de estimar a média de idade da população pelo método de bootstrap?

O método funciona da seguinte forma:

1. Reamostragem.
2. Cálculo da média dessa reamostragem.
3. Salvar a média em uma lista.
4. Repetir esse processo n vezes.
5. Por fim, calcular a média das médias do método de bootstrap.

Quanto mais vezes se repete o método de bootstrap, mais precisa fica a estimativa de média populacional. Realizarei uma estimativa utilizando 10.000 iterações desse processo. Em outras palavras, farei o processo de reamostragem 10.000 vezes. A título de demonstração, realizarei uma iteração do processo de bootstrap abaixo:

reamostragem <- copy(slice_sample(amostra, prop = 1, replace = TRUE)[order(pessoa)])

reamostragem %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")

A média de idade dessa reamostragem é de 30.8! Salvo esse valor, e repito muitas vezes.

reamostragens <- data.table(medias = replicate(10000, mean(sample(amostra$idade, length(amostra$idade), replace = TRUE))))

reamostragens %>%
  kbl(align = "cc") %>%
  kable_material() %>%
  scroll_box(width = "800px", height = "350px")