Lq1 <- function(x, p) {
if (!is.numeric(x) | !is.numeric(p))
stop("Os argumentos 'x' e 'p' precisam ser numéricos.")
n = length(x) # Captando o tamanho da amostra.
f1 = ((1 - p)^(sum(x))) * (p^n) # Definindo L(p|X), a função de verossimilhança.
return(f1) # Retornando o valor da função calculado com base no vetor passado em 'x'.
}GET00130 — Resolução da Lista 03
Questão 01
Relembremo-nos, primeiro, o que é uma função de verossimilhança.
def. Dada uma amostra aleatória \(X_{1}, X_{2},...,X_{n}\) de uma população \(X\), com função de probabilidade \(P(X=x|p)\), a função de verossimilhança de \(p\) será definida como:
\[ \begin{align*} \mathcal{L}(p|x_i) & = \mathcal{L}(p|x_1, x_2,...,x_n) \\ & = P(X=x_1|p)\times P(X=x_2|p) \times ... \times P(X=x_n|p) \\ & = \prod_{i = 1}^{n} P(X=x_i|p) \end{align*} \]
Vamos, agora, aplicar isso à função de distribuição dada na questão.
\[ \begin{align*} \mathcal{L}(p|x_i) & = \prod_{i=1}^{n} P(X=x_i|p) \\ & = (1-p)^{x_{1}}p \times (1-p)^{x_{2}}p \times ... \times (1-p)^{x_{n}}p \\ & = (1-p)^{x_1 + x_2 +...+ x_n} \times p^{n} \\ & = (1-p)^{\sum_{i=1}^{n}{x_i}} \times p^{n} \end{align*} \] A estimativa para o estimador de máxima verossimilhança para nossa variável \(X\) nada mais é do que o ponto que maximiza a função de verosimilhança dessa variável. Em termos matemáticos,
\[ \begin{align*} \hat{p} = {\operatorname {arg max}} \, \mathcal{L}(p|\textbf{X}),\ \textbf{X} = X_1, X_2,..., X_n \end{align*} \]
No R, é possível obter esse valor através de passos simples.
Primeiro, é necessário criar uma função que calcule o valor da verossimilhança e isso será feito da seguinte forma…
Vamos avaliar o comportamento dessa função, utilizando a amostra dada no enunciado e um intervalo de plotagem condizente com o parâmetro — que é uma probabilidade, ou seja, limitado no intervalo \([0, 1]\).
options(scipen = 999) # Comando para prevenir a notação científica em números muito grandes ou muito pequenos.
aa_q1 <- c(3, 0, 1, 3, 1, 1, 0, 0, 0, 0, 0, 0, 0, 3, 1) # Vetor de amostras dadas no enunciado da questão.
(gLq1 <- ggplot2::ggplot(data = dplyr::tibble(x = c(0, 1)),
mapping = ggplot2::aes(x = x)) +
ggplot2::stat_function(fun = Lq1,
args = list(x = aa_q1)) +
ggplot2::labs(y = expression(L(p)),
x = expression(p))) +
ggplot2::theme_minimal()
Por último, é necessário encontrar de fato o argumento que maximiza a função de verossimilhança, isto é, encontrar seu ponto de máximo. Essa tarefa será realizada com os comandos abaixo.
(estMVq1 <- optimize(f = Lq1,
interval = c(0, 1),
maximum = T,
x = aa_q1))$maximum
[1] 0.5357114
$objective
[1] 0.00000000400136Portanto, a estimativa para o estimador de máxima verossimilhança de \(p\) é, aproximadamente, \(0.54\). Isto é:
\[ \begin{align*} \hat{p} & = {\operatorname {arg max}} \, \mathcal{L}(p|\textbf{X}),\ \textbf{X} = X_1, X_2,..., X_n \\ & = 0.5357114 \end{align*} \] Vejamos onde se encontra essa estimativa no gráfico de \(\mathcal{L}(p|\textbf{X})\).
(gLq1_2 <- gLq1 +
ggplot2::geom_vline(xintercept = estMVq1$maximum,
color = "red2")) +
ggplot2::theme_minimal()
Questão 03
Item a
Queremos averiguar, de maneira numérica, se \(E(T_1) = E(\bar{X}) = 20\) e se \(E(T_2) = E(2\bar{X}) = 20\).
Vejam que, sobre a população, só sabemos de sua distribuição, mas não temos necessariamente uma amostra. Então, façamos um código que gere \(N\) amostras de tamanho \(n = 5\) a partir de uma \(U(0, 20)\). Aqui, vamos considerar \(N = 50000\), mas isso fica a critério pessoal.
# Criando uma matriz nXN de NAs para posterior alocação das amostras.
amostrasQ3a <- matrix(data = NA,
ncol = 50000,
nrow = 5)
# Sorteando e alocando as amostras na matriz criada acima.
for (i in 1:50000) {
amostrasQ3a[,i] = runif(n = 5, min = 0, max = 20)
}O próximo e último passo consiste em calcular as médias das amostras e a média dessas médias, para, então, comparar com o valor assumido pelo parâmetro, que é 20. Façamos isso:
# Calculando a média de cada amostra.
X_barra_Q3aT1 <- apply(X = amostrasQ3a,
MARGIN = 2,
FUN = mean)
X_barra_Q3aT2 <- 2*X_barra_Q3aT1
# Calculando a média das médias de cada estimador acima.
media_X_barra_Q3aT1 <- mean(X_barra_Q3aT1)
media_X_barra_Q3aT2 <- mean(X_barra_Q3aT2)
# Visualizando os resultados.
dplyr::tibble("T1" = media_X_barra_Q3aT1,
"T2" = media_X_barra_Q3aT2)# A tibble: 1 × 2
T1 T2
<dbl> <dbl>
1 10.0 20.0Os resultados indicam que \(T_1 = \bar{X}\) é um estimador viesado para \(\theta\), enquanto \(T_2 = 2\bar{X}\) é não-viesado, devido a alta precisão da aproximação de \(E(T_2)\) com o valor de \(\theta\).
A título de curiosidade, \(T_1\) é viesado apenas para \(\theta\). Vejam que, no contexto de estimação da média populacional, \(\mu = \frac{0+20}{2} = 10\), \(E(T_1) = E(\bar{X}) = 10 = \mu\) e, portanto, \(T_1\) é não viesado para a estimação da média populacional.
Item b
Agora, queremos averiguar a consistência desses dois estimadores. O processo é similar ao do item anterior, acrescido de alguns passos para poder calcular a variância das médias amostrais. Portanto, concentremo-nos no código.
# Definindo o número de amostras a serem sorteadas.
N = 50000
# Definindo os diferentes tamanhos dessas amostras
n = c(2, 5, 10, 50, 100, 500, 1000, 10000)
# Definindo os objetos que armazenarão a variância das médias dessas amostras.
varT1 = NULL
varT2 = NULL
# Criando a matriz que armazenará as amostras.
for (i in 1:length(n)) {
amostrasQ3b = matrix(data = NA,
ncol = N,
nrow = n[i])
for (j in 1:N) {
amostrasQ3b[,j] = runif(n = n[i], min = 0, max = 20)
}
# Calculando a média para cada amostra obtida
X_barra_Q3bT1 = apply(X = amostrasQ3b,
MARGIN = 2,
FUN = mean)
X_barra_Q3bT2 = 2*X_barra_Q3bT1
# Calculando a variância dessas médias
varT1[i] = var(X_barra_Q3bT1)
varT2[i] = var(X_barra_Q3bT2)
}
# Visualizando os resultados obtidos.
dplyr::tibble("Tamanho da amostra (n)" = n,
"Variância de T1" = varT1,
"Variância de T2" = varT2)# A tibble: 8 × 3
`Tamanho da amostra (n)` `Variância de T1` `Variância de T2`
<dbl> <dbl> <dbl>
1 2 16.8 67.0
2 5 6.66 26.6
3 10 3.32 13.3
4 50 0.671 2.68
5 100 0.336 1.35
6 500 0.0668 0.267
7 1000 0.0337 0.135
8 10000 0.00331 0.0132A variância dos dois estimadores diminui, tendendo a zero, conforme o tamanho de \(n\) aumenta. Entretanto, dado o vício de \(T_1\) comprovado no item anterior, apenas \(T_2\) é um estimador consistente. Notem que uma amostra de tamanho \(n = 5\) reduz drasticamente a variância desse estimador, quando comparado a uma amostra de tamanho \(n = 2\). Veja no gráfico logo abaixo.
g1 <- dplyr::tibble(
"Tamanho da amostra (n)" = n,
"Variância de T1" = varT1,
"Variância de T2" = varT2
) |>
ggplot2::ggplot(mapping = ggplot2::aes(x = factor(n),
y = `Variância de T1`)) +
ggplot2::geom_bar(stat = "identity") +
ggplot2::labs(x = "Tamanho da amostra (n)",
y = "Variância de T1") +
ggplot2::theme_minimal()
g2 <- dplyr::tibble(
"Tamanho da amostra (n)" = n,
"Variância de T1" = varT1,
"Variância de T2" = varT2
) |>
ggplot2::ggplot(mapping = ggplot2::aes(x = factor(n),
y = `Variância de T2`)) +
ggplot2::geom_bar(stat = "identity") +
ggplot2::labs(x = "Tamanho da amostra (n)",
y = "Variância de T2") +
ggplot2::theme_minimal()
library(patchwork)
g1 + g2
Questão 5
Temos \(\textbf{X} = \{71, 66, 81, 69, 68, 72, 77, 76, 80, 75, 75, 81, 78, 69, 79, 75, 78, 72, 84, 77\}\) uma a.a. de \(X \sim N(\mu,\sigma^2)\).
Item a
Neste item, queremos o \(IC\) para a média de população normal que não possui variância conhecida, considerando \(97\%\) de confiança, e tendo como base a amostra observada. Ou seja, queremos:
\[ IC_{(\mu,1-\alpha)} = \biggl(\bar{X}-t_{(n-1;\alpha/2)} \frac{S}{\sqrt{n}}; \bar{X}+t_{(n-1;\alpha/2)} \frac{S}{\sqrt{n}} \biggr) \]
Para isso, façamos:
# Criando o vetor da amostra de X
amostraQ5 <- c(71, 66, 81, 69, 68, 72, 77, 76, 80, 75, 75, 81, 78, 69, 79, 75, 78, 72, 84, 77)
# Calculando o IC
Publish::ci.mean(x = amostraQ5,
alpha = .03) mean CI-97%
75.15 [72.58;77.72]A interpretação padrão é que, dada a construção de diversos intervalos para diversas amostras de mesmo tamanho, oriundas de uma população normal, \(97\%\) desses intervalos vão conter o peso médio desses 20 indivíduos.
Item c
Vamos definir uma variável binária \(Y = Classificação\ do\ peso\ do\ indivíduo\). Se o indivíduo possui um peso maior que \(75 \ kg\), \(Y = 1\). Caso contrário, \(Y = 0\). No R, isso fica da seguinte forma:
tibbleQ5 <- dplyr::tibble("Amostra" = amostraQ5) |>
dplyr::mutate("Classificação" = dplyr::case_when(Amostra > 75 ~ 1,
TRUE ~ 0))Na construção acima, optei por utilizar a função case_when() em vez da cut(), no intuito de mostrar uma alternativa para criação de colunas que dependem de outras. Seu uso não é tão intuitivo, mas após um pouco de prática e leitura do manual, ele se torna.
O primeiro argumento, se não tivessemos utilizado o pipe, seria o tibble de interesse. Os próximos argumentos se baseiam em testes de condições.
Especificamente no caso acima, eu testo, primeiro, se os valores da coluna Amostra são maiores que 75 e, quando forem, atribuo 1 a linha referente a esse valor na coluna nova, Classificação.
Por último, o segundo argumento avalia se os valores da coluna Amostra são menores ou iguais a 75, mas ele o faz de uma maneira um pouco diferente. Vejam que o valor 1 só é atribuído na coluna Classificação quando Amostra > 75 == TRUE. Dessa forma, quando esse condicional retorna FALSE, é necessário que haja algo para atribuir. O formato do segundo argumento foi moldado para que, quando (Amostra > 75 == FALSE) == TRUE, o valor 0 seja inputado na coluna Classificação. Em resumo, o segundo argumento atribuiu um valor com base no resultado do teste lógico feito sobre o teste lógico do argumento anterior.
Feito isso, podemos dizer que o que queremos é, na verdade o \(IC_{(p,1-\alpha)}\). Isto é, visto a partir da perspectiva otimista, considerando o nível de confiança de \(95\%\):
\[ \begin{align*} IC_{(p,95\%)} = \Biggr(\hat{p} - z_{0,05/2} \sqrt{\frac{p(1-p)}{n}}; \hat{p} + z_{0,05/2} \sqrt{\frac{p(1-p)}{n}} \Biggl) \end{align*} \] Façamos isso no R, então.
binom::binom.asymp(x = sum(tibbleQ5$`Classificação`),
n = length(tibbleQ5$`Classificação`),
conf.level = .95) |>
dplyr::as_tibble()# A tibble: 1 × 6
method x n mean lower upper
<chr> <dbl> <int> <dbl> <dbl> <dbl>
1 asymptotic 10 20 0.5 0.281 0.719Fazendo uso da interpretação padrão, dada a construção de diversos intervalos para diversas amostras de mesmo tamanho, \(95\%\) desses intervalos vão conter a proporção de indivíduos com peso acima de \(75 \ kg\).