ESTATÍSTICA COMPUTACIONAL I

PRÁTICA 2: ESTIMADORES PONTUAIS E PROPRIEDADES

Prof: Guilherme Augusto Veloso (guilhermev@id.uff.br)

\[\\[0.05in]\]

1 PACOTES NECESSÁRIOS

Para esta terceira aula prática de Estatística Computacional I, precisaremos dos seguintes pacotes:

require(ggplot2) # Abordagem gráfica
require(patchwork) # Múltiplos gráficos
require(fitdistrplus) # Método dos Momentos e da máxima verossimilhança
library(plotly) # Gráficos 3D

2 INTRODUÇÃO

A inferência estatística busca tirar conclusões sobre parâmetros populacionais a partir de uma amostra observada. A estimação pontual é uma das abordagens fundamentais da inferência, focando em fornecer um único valor (a “melhor aposta”) para o parâmetro desconhecido da população.

Imagine uma população real, como:

todos os eleitores de um país,
todos os carros fabricados por uma montadora,
todos os pacientes acometidos por determinada doença.

Geralmente, temos interesse em parâmetros que descrevem essa população:

média \(\mu\) da renda dos eleitores,
proporção \(p\) de carros com defeito,
desvio padrão \(\sigma\) da pressão alta dos pacientes.

Contudo, não temos acesso à população inteira. Em vez disso, trabalhamos com uma amostra aleatória e, a partir dela, estimamos os parâmetros desejados. Para além da teoria, a utilização de ferramentas computacionais como o R permite explorar os estimadores de maneira empírica. Com simulações, poderemos comparar propriedades dos estimadores como o viés e consistência, e experimentar diferentes métodos em cenários controlados, o que fortalece a compreensão dos conceitos e estimula a autonomia na análise de dados.

2.1 NOTAÇÃO

Antes de começarmos os nossos estudos computacionais envolvendo estimadores pontuais e proprieddes, é importante definir a notação Estatística utilizada nesta disciplina

2.1.1 VARIÁVEL ALEATÓRIA E PARÂMETRO

Seja \(X\) uma variável aleatória que descreve alguma característica sobre a população de interesse, com distribuição de probabilidade \(F\).

O parâmetro \(\theta\) é uma característica numérica de \(F\). Exemplos:

\(\theta=\mu=\mathbb{E}[X]\): média populacional
\(\theta=\sigma^2=\mathrm{Var}(X)\): variância populacional
\(\theta=p=\mathbb{P}(X = 1)\): proporção populacional

2.1.2 AMOSTRA ALEATÓRIA

Uma amostra aleatória simples de tamanho \(n\) é um vetor de variáveis aleatórias:

\[ (X_1, X_2, \ldots, X_n) \]

tal que \(X_1, X_2, \ldots, X_n\) são i.i.d. (independentes e identicamente distribuídas) com distribuição \(F\)

A realização observada da amostra é:

\[ (x_1, x_2, \ldots, x_n) \]

2.1.3 ESTIMADOR E ESTIMATIVA

Um estimador de um parâmetro \(\theta\) é uma função da amostra:

\[ \hat{\theta} = g(X_1, X_2, \ldots, X_n) \]

Ou seja, é uma variável aleatória que tenta aproximar o parâmetro \(\theta\).

Exemplos:

Estimador da média: \(\hat{\mu} = \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i\)
Estimador da variância: \(\hat{\sigma}^2 = \frac{1}{n - 1} \sum_{i=1}^{n} (X_i - \bar{X})^2\)

A estimativa é o valor numérico obtido ao aplicar o estimador a uma realização da amostra:

\[ \hat{\theta}_{obs} = g(x_1, x_2, \ldots, x_n) \]

Resumindo:

Estimador = função das variáveis aleatórias (tem variabilidade)
Estimativa = número fixo obtido com dados reais

3 ESTIMADORES PONTUAIS

Os estimadores pontuais fornecem, a partir de uma amostra, um único valor que represente a nossa melhor estimativa para um parâmetro. Para isso, utilizamos estimadores, que nada mais são do que funções dos dados amostrais. Basicamente, na inferência estatística, aprendemos três métodos de estimação pontual:

Método dos Momentos
Método da Máxima Verossimilhança
Método dos Mínimos Quadrados

Nas próximas seções veremos os detalhes computacionais dos dois primeiros métodos. O método dos mínimos quadrados será visto na disciplina de modelos lineares.

No R, podemos obter os estimadores de momentos e de máxima verossimilhança para diferentes distribuições de probabilidade a partir da função fitdist. Essa função tem os seguintes argumentos:

data: dados amostrais.
distr: distribuição de probabilidade.
method: mme para o método dos momentos e mle para o método da máxima verossimilhança.

Para mais detalhes:

?fitdist

3.1 MÉTODO DOS MOMENTOS

O método dos momentos é uma técnica clássica de estimação de parâmetros populacionais com base nas amostras. Ele foi proposto por Karl Pearson no final do século XIX e é uma das formas mais simples de construir estimadores. A ideia central do método dos momentos é igualar os momentos amostrais com os momentos populacionais teóricos e, a partir dessa equação, obter estimativas para os parâmetros da distribuição.

3.1.1 Definições

O \(r\)-ésimo momento populacional é definido por:

\[ \mu_r = \mathbb{E}[X^r] \]

O \(r\)-ésimo momento amostral é dado por:

\[ m_r = \frac{1}{n} \sum_{i=1}^{n} X_i^r \]

3.1.2 Procedimento Clássico

Suponha que desejamos estimar \(k\) parâmetros desconhecidos \(\theta_1, \dots, \theta_k\). O procedimento do método dos momentos é o seguinte:

Calcular os \(k\) primeiros momentos amostrais: \(m_1, m_2, \dots, m_k\).
Obter os \(k\) primeiros momentos teóricos (populacionais) como funções dos parâmetros: \(\mu_1, \mu_2, \dots, \mu_k\).
Igualar os momentos amostrais aos momentos teóricos: \[ m_r = \mu_r, \quad \text{para } r = 1, \dots, k \]
Resolver o sistema de equações para obter os estimadores dos parâmetros: \[ \hat{\theta}_1, \dots, \hat{\theta}_k \]

3.1.3 Procedimento Alternativo: Momentos Centrados

Em muitas distribuições, em vez de trabalhar diretamente com os momentos brutos, é mais conveniente utilizar os momentos centrados (média e variância). O procedimento é análogo ao clássico:

Calcular os momentos centrados amostrais:
- Média amostral: \[ \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i \]
- Segundo momento centrado amostral (variância não corrigida): \[ s_v^2 = \frac{1}{n} \sum_{i=1}^{n} (X_i - \bar{X})^2 \]
Obter os momentos centrados populacionais em função dos parâmetros:
- Média populacional: \[ \mu = \mathbb{E}[X] \]
- Variância populacional: \[ \sigma^2 = \mathbb{E}[(X - \mu)^2] \]
Igualar momentos centrados amostrais e populacionais: \[ \bar{X} = \mu, \quad s_v^2 = \sigma^2 \]
Resolver o sistema para encontrar os estimadores dos parâmetros.

3.1.4 Por que essa estratégia funciona?

A motivação é que média e variância são casos particulares de momentos centrados.

A igualdade \[ \bar{X} \approx \mathbb{E}[X] = \mu \] garante que o primeiro momento centrado da amostra aproxime o teórico.
A igualdade \[ s_v^2 \approx \mathbb{E}[(X - \mu)^2] = \sigma^2 \] garante que o segundo momento centrado da amostra aproxime o teórico.

Assim, ao igualar média e variância amostrais às populacionais, obtemos estimadores consistentes dos parâmetros em muitas distribuições usuais (Normal, Binomial, Poisson, Gama, etc.).

Esse procedimento é apenas uma reformulação do método dos momentos, mas geralmente mais intuitivo, pois média e variância já são medidas naturais de posição e dispersão.

3.1.5 Exemplo: Distribuição Gama

Seja \(X \sim \text{Gamma}(\alpha, \beta)\), com \(\alpha > 0\) (parâmetro de forma) e \(\beta > 0\) (parâmetro de taxa). Temos dois parâmetros para serem estimados: \(\alpha\) e \(\beta\).

Sabemos que:
\[ E(X) = \frac{\alpha}{\beta}, \quad Var(X) = \frac{\alpha}{\beta^2} \]

Considere que uma amostra de tamanho \(n\) \((X_1, X_2, \ldots, X_n)\) foi obtida dessa variável.

De acordo com o método dos momentos centrados, temos:

\[ \bar{X} = \frac{\alpha}{\beta}, \quad s_v^2 = \frac{\alpha}{\beta^2} \]

Logo, o sistema a resolver é:

\[ \begin{cases} \frac{\hat{\alpha}}{\hat{\beta}} = \bar{X} \\ \frac{\hat{\alpha}}{\hat{\beta}^2} = s_v^2 \end{cases} \]

A solução é:

\[ \hat{\alpha} = \frac{\bar{X}^2}{s_v^2}, \quad \hat{\beta} = \frac{\bar{X}}{s_v^2} \]

Considere a situação de uma clínica que quer estimar o tempo de atendimento por paciente, desde o momento em que a pessoa chega até a saída, para planejar a alocação de pessoal. Estudos anteriores sugerem que o tempo de atendimento (em horas) é assimétrico à direita — muitos atendimentos são rápidos, mas alguns demoram mais (por exemplo, em casos de crianças, idosos, pessoas com medo de agulha etc). Por isso, consideramos adequado modelar o tempo de atendimento com uma distribuição Gama. Foram coletados 100 valores de tempos de atendimento, em horas. Os dados estão apresentados abaixo:

Tempos = c(
  1.98, 0.60, 1.28, 1.57, 1.71, 1.34, 1.69, 1.14, 1.64, 0.43,
  1.11, 0.92, 2.40, 1.14, 1.32, 1.54, 2.36, 2.56, 3.76, 3.05,
  1.36, 1.88, 2.29, 0.33, 0.92, 1.17, 0.79, 1.49, 2.49, 0.73,
  1.78, 1.30, 1.01, 0.74, 0.88, 0.92, 0.85, 1.01, 1.16, 0.39,
  1.00, 1.90, 0.57, 2.31, 1.62, 1.97, 0.39, 1.67, 0.79, 1.98,
  1.72, 0.47, 1.87, 1.33, 0.90, 3.72, 3.46, 1.91, 2.85, 1.80,
  1.85, 0.55, 0.90, 1.51, 1.57, 2.37, 1.79, 0.48, 0.87, 1.47,
  1.85, 1.24, 0.32, 2.77, 1.69, 0.59, 0.84, 2.36, 2.03, 1.26,
  0.40, 0.63, 0.59, 1.50, 0.52, 0.20, 0.25, 3.09, 3.72, 3.16,
  1.93, 0.85, 1.87, 1.66, 1.31, 1.63, 1.94, 0.20, 2.72, 3.57
)

Vamos calcular as estimativas de \(\alpha\) e \(\beta\) segundo o método dos momentos:

# Tamanho da amostra
n = length(Tempos)
# Media amostral
xbarra = mean(Tempos)
# Segundo momento amostral centrado na média
s_2_v = var(Tempos)*(n-1)/n

# Estimativas pelo método dos momentos
Alfa_MM = xbarra^2/s_2_v
Beta_MM = xbarra/s_2_v

Alfa_MM

## [1] 3.125629

Beta_MM

## [1] 2.065712

Agora, podemos usar a função fitdist para conferir os nossos resultados:

# Função automática do R
fitdist(data=Tempos,distr="gamma",method="mme")$estimate

##    shape     rate 
## 3.125629 2.065712

Conforme esperado, os valores são iguais.

3.1.6 EXERCÍCIO: DISTRIBUIÇÃO UNIFORME

Uma confeitaria artesanal produz barras de chocolate e, devido a pequenas variações no processo manual de produção, o peso pode variar. A empresa garante que qualquer peso é igualmente provável. Ou seja, o peso da barra pode ser modelado por uma variável aleatória uniforme (a,b). Para fazer a inferência sobre os limites \(a\) e \(b\), foi coletada uma amostra de 50 barras de chocolate e foi verificado o peso de cada uma delas. Os dados estão descritos abaixo:

Peso = c(97.77, 101.74, 107.97, 102.23, 93.34, 105.71, 100.51,
  103.90, 93.34, 109.43, 103.24, 96.71, 93.80, 103.36,
  98.99, 93.22, 102.47, 109.30, 94.99, 105.54, 91.84,
  108.32, 101.80, 109.16, 99.82, 108.04, 90.53, 95.29,
  101.64, 100.59, 100.75, 107.65, 97.48, 90.91, 98.79,
  96.31, 96.82, 96.56, 107.75, 101.14, 105.16, 106.20,
  95.98, 105.90, 91.40, 102.19, 99.40, 93.33, 102.95,
  93.65)

Dados importantes:

Se \(X \sim \text{Uniforme}(a, b)\), então:

\[f(x) = \begin{cases} \dfrac{1}{b - a}, & \text{se } a \leq x \leq b \\ 0, & \text{caso contrário} \end{cases}\]

\[\mathbb{E}(X) = \frac{a + b}{2}\]

\[\text{Var}(X) = \frac{(b - a)^2}{12}\]

Encontre as estimativas de \(a\) e \(b\) pelo método dos momentos e confira o seu resultado com o da função fitdist.

3.2 MÉTODO DA MÁXIMA VEROSSIMILHANÇA

Imagine que você está tentando descobrir qual parâmetro de uma distribuição é o mais provável de ter gerado os dados que você tem em mãos. O método da máxima verossimilhança permite justamente isso: a partir dos dados observados, ele nos fornece um valor do parâmetro que torna a ocorrência dos dados mais “verossímil”. Essa ideia é fundamental em estatística inferencial e está por trás de muitas aplicações em ciência de dados, aprendizado de máquina, economia, biologia e diversas outras áreas.

Por exemplo, suponha que você está estudando o tempo de vida de peças de uma máquina, que você acredita seguir uma distribuição exponencial. Após observar uma amostra de tempos de vida, você quer descobrir qual valor de \(\lambda\) torna mais plausível a ocorrência dos dados observados. A resposta vem pela maximização de uma função: a função de verossimilhança.

Essa função mede o quanto os dados combinam com diferentes valores do parâmetro. O valor que maximiza essa função é então aquele que melhor explica os dados observados, no sentido probabilístico.

A função de verossimilhança é dada por:

\[ L(\theta) = \prod_{i=1}^{n} f(x_i; \theta) \]

em que \(f(x_i;\theta)\) é a função de probabilidade ou função densidadede probabilidade avaliada no valor \(x_i\). Trata-se de uma função de \(\theta\), com os dados \(x_1, \ldots, x_n\) fixos. Para facilitar os cálculos, utilizamos a log-verossimilhança:

\[ \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log f(x_i; \theta) \]

Para encontrar o estimador de máxima verossimilhança, precisamos maximizar \(L(\theta)\) ou \(\ell(\theta)\), com respeito a \(\theta\). É mais fácil maximizar a função de log-verossimilhança do que a função de verossimilhança tanto analiticamente quanto computacionalmente. Ambas tem o mesmo ponto de máximo, uma vez que a função logaritmica é uma função crescente e as funções de densidade e probabilidade são positivas. Assim, o estimador de máxima verossimilhança é dado por:

\[ \hat{\theta}_{\text{MV}} = \arg\max_{\theta \in \Theta} \ell(\theta) \]

Para encontrar \(\hat{\theta}_{\text{MV}}\), buscamos os pontos críticos:

\[ \frac{d}{d\theta} \ell(\theta) = 0 \]

No caso vetorial:

\[ \nabla_\theta \ell(\theta) = 0 \]

em que são necessários os cálculos de todas as derivadas parciais. É fundamental verificar se o ponto encontrado é um máximo e se está no interior do espaço paramétrico.

3.2.1 EXEMPLO: DISTRIBUIÇÃO EXPONENCIAL

Seja \(X_1, X_2, \ldots, X_n\) uma amostra aleatória simples de uma população com distribuição exponencial com parâmetro \(\lambda\), cuja função densidade de probabilidade é dada por:

\[ f(x; \lambda) = \lambda e^{-\lambda x}, \quad x \geq 0 \]

A função de verossimilhança é o produto das densidades avaliadas nos dados observados:

\[ L(\lambda) = \prod_{i=1}^{n} f(x_i; \lambda) = \prod_{i=1}^{n} \lambda e^{-\lambda x_i} = \lambda^n e^{-\lambda \sum_{i=1}^{n} x_i} \]

A função log-verossimilhança é dada por:

\[ \ell(\lambda) = \log L(\lambda) = \log \left( \lambda^n e^{-\lambda \sum x_i} \right) = n \log \lambda - \lambda \sum_{i=1}^{n} x_i \]

Sabemos que o estimador obtido pelo método da máxima verossimilhança para \(\lambda\) é dado por:

\[\hat{\lambda}_{MV}=\frac{1}{\bar{x}}\]

Podemos usar o R para provar computacionalmente que este estimador maximiza \(L(\lambda)\) e \(\ell(\lambda)\) de forma gráfica. Para isto, considere uma fábrica que possua uma linha de montagem completamente automatizada. Engenheiros perceberam que as falhas nessa linha ocorrem de forma esporádica e são independentes entre si. Após observar o sistema por muitos dias, eles verificam que o tempo entre duas falhas consecutivas pode ser modelado por uma variável aleatória exponencial com parâmetro \(\lambda\). Para estimar esse parâmetro a partir do método de máxima verossimilhança, eles coletaram 34 tempos, em horas, entre falhas:

# Amostra
Tempos = c(0.47, 0.27, 0.12, 0.13, 0.16, 0.16, 0.38, 0.59,
  0.04, 0.12, 0.25, 0.20, 0.21, 0.31, 0.34, 0.26,
  0.08, 0.20, 0.45, 0.18, 0.15, 0.29, 0.25, 0.25,
  0.06, 0.08, 0.07, 0.11, 0.02, 0.35, 0.17, 0.08,
  0.15, 0.21)

Sabemos que o estimador de máxima verossimilhança para \(\lambda\) é \(\hat{\lambda}_{MV}=\frac{1}{\bar{x}}\). Vamos calcular ele no R:

# Estimativa de máxima verossimilhança
rate_MV = 1/mean(Tempos)

Assim \(\hat{\lambda}_{MV}=4.74\) é a estimativa de máxima verossimilhança de \(\lambda\). Esse número deve maximizar as funções de verossimilhança ou log-verossimilhança, atreladas aos dados. Para iniciar, vamos implementar as funções de verossimilhança e log verossimilhança no R:

# Função de verossimilhança
Vero_Exp <- function(rate, x) {
  prod(dexp(x=x,rate=rate))
}

# Função de log-verossimilhança
Log_Vero_Exp <- function(rate,x){
  sum(dexp(x=x,rate=rate,log=T))
}

Note que a função de verossimilhança é o produtório de densidades e a função de log-verossimilhança é o somatório do logaritmo das densidades (por isso usamos o argumento log=T para usar a escala logaritmica). Vamos avaliar ambas as funções graficamente e verificar como a estimativa de máxima verossimilhança \(\hat{\lambda}_{MV}=4.74\) maximiza ambas as funções

# Taxas arbitrárias
rates = seq(0.01,10,length.out = 500)

# Data Frame com os valores das funções aplicado nas taxas
df = data.frame(rates = rates,
                verossim_vals = sapply(rates, Vero_Exp, x = Tempos),
                logverossim_vals = sapply(rates, Log_Vero_Exp, x = Tempos))

# Gráfico da função de verossimilhança
A = ggplot(df, aes(x = rates, y = verossim_vals)) +
  geom_line(color = "blue", size = 1) +
   geom_point(aes(x = rate_MV, y = Vero_Exp(rate_MV, Tempos)),
             color = "black", size = 3)+
  geom_vline(xintercept = rate_MV, linetype = "dashed", color = "blue") +
  labs(title = "Função de Verossimilhança",
       x = "Lambda (rate)",
       y = "Verossimilhança") +
  theme_minimal()

# Gráfico da função de log-verossimilhança
B = ggplot(df, aes(x = rates, y = logverossim_vals)) +
  geom_line(color = "darkred", size = 1) +
  geom_point(aes(x = rate_MV, y = Log_Vero_Exp(rate_MV, Tempos)),
             color = "black", size = 3)+
  geom_vline(xintercept = rate_MV, linetype = "dashed", color = "blue") +
  labs(title = "Função de Log-Verossimilhança",
       x = "Lambda (rate)",
       y = "Log-Verossimilhança") +
  theme_minimal()

A+B

Temos, portanto, o gráfico da função de verossimilhança (esquerda) e o gráfico da função de log-verossimilhança (direita). Para finalizar, vamos conferir com a função automática do R:

# Função automática do R
fitdist(data=Tempos,distr="exp",method="mle")$estimate

##     rate 
## 4.748603

O R utiliza a função optim dentro da função fitdist para maximizar a função de log-verossimilhança.

3.2.2 EXERCÍCIO: DISTRIBUIÇÃO POISSON

Suponha que estamos estudando a chegada de clientes em uma farmácia de plantão 24h em um bairro tranquilo. A gerência quer entender melhor como otimizar o número de atendentes durante a madrugada entre meia-noite e 1h da manhã. Um funcionário coletou dados de 36 madrugadas neste horário e anotou os seguintes resultados:

# Amostra
Clientes = c(9, 5, 6, 5, 2, 8, 3, 2, 4, 9, 8, 6, 6, 11, 6, 6, 5, 5, 4, 3, 9, 8, 6, 7, 1, 5, 6, 3, 4, 3, 3, 4, 4, 4, 3, 3)

Como a unidade de tempo é fixa (1 hora) e os eventos ocorrem de forma independente, considerou-se que o número de clientes que chegam nessa horário é uma variável aleatória Poisson com parâmetro \(\lambda\).

Encontre a estimativa de máxima verossimilhança de \(\lambda\).
Plote as funções de verossimilhança e log-verossimilhança, atrelada aos dados.
Marque nestas funções, a estimativa obtida no item 1-. Qual característica ele atende?
Confira seus resultados com a função fitdist do R.

3.2.3 EXEMPLO: DISTRIBUIÇÃO UNIFORME

Seja \(X_1, X_2, \ldots, X_n\) uma amostra aleatória simples de uma população com distribuição uniforme (0,\(\theta\)), cuja função densidade de probabilidade é dada por:

\[ f(x; \theta) = \begin{cases} \frac{1}{\theta}, & \text{se } 0 \le x \le \theta \\ 0, & \text{caso contrário} \end{cases} \]

A função de verossimilhança é o produto das densidades avaliadas nos dados observados:

\[ L(\theta) = \prod_{i=1}^{n}f(x_i; \theta) = \prod_{i=1}^{n} \frac{1}{\theta} \cdot \mathbb{I}_{[0, \theta]}(x_i) = \theta^{-n} \cdot \mathbb{I}_{[X_{(n)} \le \theta]} \]

Como \(\theta^{-n}\) é uma função decrescente e \(\theta \geq x_{(n)}\), o estimador de máxima verossimilhança de \(\theta\) é:

\[ \hat{\theta}_{\text{MV}} = X_{(n)} = \max\{X_1, X_2, \dots, X_n\} \]

Para ilustrar este resultado, suponha uma empresa de tecnologia que desenvolve um sistema de atendimento automático. A equipe de engenharia afirma que o tempo de resposta do sistema (em segundos), quando está operando corretamente, é distribuído de maneira uniforme entre 0 e um valor máximo \(\theta\), desconhecido. Por simplicidade, o tempo de resposta nunca é negativo (começa em 0) e nunca excede esse limite superior desconhecido \(\theta\). Foram registrados os tempos de resposta de 28 atendimentos e foram observados os seguintes valores:

# Amostra
Tempos = c(7.46, 14.91,  8.51, 27.45,  1.47, 14.15,
 25.57,  3.90, 17.95,  6.61,  4.08, 24.11,
 28.64, 11.98, 21.28,  3.03, 12.29,  8.78,
 26.07, 14.35, 25.92, 26.00, 25.42, 14.07,
 24.14, 20.14, 22.73,  0.02)

Sabemos que o estimador de máxima verossimilhança para \(\theta\) é o valor máximo \(X_{(28)}\):

# Estimativa de máxima verossimilhança
Theta_MV = max(Tempos)

Assim, \(\hat{\theta}_{\text{MV}}=28,64\) é a estimativa de máxima verossimilhança para \(\theta\). Agora, vamos implementar a função de máxima verossimilhança:

# Função de verossimilhança
Vero_Unif <- function(theta, x) {
  n <- length(x)
  if (theta < max(x)) 
    return(0)
  else
    return((1/theta)^n)
}

Agora, vamos avaliar esta função graficamente e marcar a estimativa de máxima verossimilhança \(\hat{\theta}_{MV}=28.64\):

# Thetas arbitrários
thetas = seq(25,40,length.out = 500)

# Data frame com o valor da função de verossimilhança
df = data.frame(thetas = thetas,
                verossim_vals = sapply(thetas, Vero_Unif, x = Tempos))

# Gráfico da função de verossimilhança
ggplot(df, aes(x = thetas, y = verossim_vals)) +
  geom_line(color = "blue", size = 1) +
   geom_point(aes(x = Theta_MV, y = Vero_Unif(Theta_MV, Tempos)),
             color = "black", size = 3)+
  geom_vline(xintercept = Theta_MV, linetype = "dashed", color = "blue") +
  labs(title = "Função de Verossimilhança",
       x = "Theta",
       y = "Verossimilhança") +
  theme_minimal()

Note como a estimativa de máxima verossimilhança maximiza esta função. Para finalizar, o comando a seguir certifica o resultado com a análise automática do R:

# Análise automática do R
fitdist(data=Tempos,distr="unif",method="mle",fix.arg=list(min = 0))$estimate

##   max 
## 28.64

O comando fix.arg=list(min = 0) fixa o primeiro parâmetro em 0.

3.2.4 EXEMPLO: DISTRIBUIÇÃO NORMAL

Seja \(X_1, X_2, \ldots, X_n\) uma amostra aleatória simples de uma população com distribuição Normal (\(\mu\),\(\sigma^2\)), cuja função densidade de probabilidade é dada por:

\[ f(x;\mu,\sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left( -\frac{(x - \mu)^2}{2\sigma^2} \right) \]

A função de verossimilhança é: \[ L(\mu, \sigma^2) = \prod_{i=1}^n f(x_i; \mu, \sigma^2) = \left( \frac{1}{\sqrt{2\pi \sigma^2}} \right)^n \exp\left( -\frac{1}{2\sigma^2} \sum_{i=1}^n (x_i - \mu)^2 \right) \]

A função de log-verossimilhança é dada por: \[ \ell(\mu, \sigma^2) = -\frac{n}{2} \log(2\pi) - \frac{n}{2} \log(\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (x_i - \mu)^2 \]

Após o processo de maximização da função de log-verossimilhança com respeito ao vetor (\(\mu,\sigma^2\)), os estimadores de máxima verossimilhança são:

\[ \hat{\mu} = \bar{x} = \frac{1}{n} \sum_{i=1}^n x_i, \quad \hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2 \]

Podemos usar o R para provar computacionalmente que este estimador maximiza \(L(\mu,\sigma^2)\) e \(\ell(\mu,\sigma^2)\) de forma gráfica. Um engenheiro agrônomo está analisando o peso das laranjas colhidas em uma fazenda para verificar se a produção está dentro dos padrões desejados. Sabe-se que o peso de cada laranja é influenciado por diversos fatores como irrigação, incidência solar, solo, nutrientes e genética — todos com efeitos pequenos e relativamente independentes. Por isso, é razoável assumir que o peso das laranjas segue uma distribuição normal com média \(\mu\) e variância \(\sigma^2\). Para estimar os parâmetros dessa distribuição, o engenheiro coleta aleatoriamente 38 laranjas da plantação e registra os seguintes pesos (em gramas):

# Amostra
Peso = c(129.7, 128.5, 128.1, 126.5, 129.0, 123.7, 140.8,
  136.0, 124.4, 128.0, 127.7, 133.9, 129.6, 131.3,
  129.9, 129.8, 136.8, 128.9, 137.6, 122.3, 132.9,
  130.6, 131.1, 131.9, 127.5, 128.3, 124.9, 124.6,
  131.5, 132.2, 130.3, 134.6, 140.3, 127.5, 118.5,
  135.0, 126.5, 126.6)

Sabemos que as estimativas de máxima verossimilhança para \(\mu\) e \(\sigma^2\)é \(\hat{\mu}_{MV}=\frac{1}{n} \sum_{i=1}^n x_i\) e \(\hat{\sigma^2}_{MV}=\frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2\). Vamos calculá-los no R:

n=length(Peso)

# Estimativa de máxima verossimilhança para mu
mean_MV = mean(Peso)

# Estimativa de máxima verossimilhança para sigma^2
sd_2_MV = var(Peso)*(n-1)/n

# Estimativa de máxima verossimilhança para sigma
sd_MV = sqrt(sd_2_MV)

Assim \(\hat{\mu}_{MV}=129,93\) e \(\hat{\sigma}_{MV}=4,68\) é a estimativa de máxima verossimilhança de \(\mu\) e \(\sigma\). Esse vetor deve maximizar as funções de verossimilhança ou log-verossimilhança, atreladas aos dados. Para iniciar, vamos implementar ambas as funções de verossimilhança e de log-verossimilhança no R:

# Função de verossimilhança
Vero_Norm <- function(mean,sd, x){
  prod(dnorm(x=x,mean=mean,sd=sd))
}

# Função de log-verossimilhança
Log_Vero_Norm <- function(mean,sd, x){
  sum(dnorm(x=x,mean=mean,sd=sd,log=T))
}

Agora, vamos fazer os gráficos das funções com o ponto (\(\hat{\mu}_{MV}=129,93\) , \(\hat{\sigma}_{MV}\)=4,68) maximizando ambas as funções (código omitido):

Para finalizar, vamos recorrer à função automática do R para conferir os resultados:

fitdist(data=Peso,distr="norm",method="mle")$estimate

##       mean         sd 
## 129.928947   4.682662

3.2.5 EXERCÍCIO: DISTRIBUIÇÃO GAMA

Seja \(x_1, x_2, \dots, x_n\) uma amostra aleatória de uma população com distribuição Gama de parâmetros \(\alpha > 0\) (forma) e \(\beta > 0\) (taxa), cuja função densidade de probabilidade é dada por:

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

A função de verossimilhança é o produto das densidades individuais:

\[ L(\alpha, \beta; x_1, \dots, x_n) = \prod_{i=1}^{n} f(x_i; \alpha, \beta) = \prod_{i=1}^{n} \left( \frac{\beta^\alpha}{\Gamma(\alpha)} x_i^{\alpha - 1} e^{-\beta x_i} \right). \]

Tomando o logaritmo da função de verossimilhança, obtemos a função log-verossimilhança:

\[ \ell(\alpha, \beta) = \log L(\alpha, \beta) = n \alpha \log \beta - n \log \Gamma(\alpha) + (\alpha - 1) \sum_{i=1}^{n} \log x_i - \beta \sum_{i=1}^{n} x_i. \]

Para encontrar os estimadores de máxima verossimilhança, não conseguimos fazer as contas na mão e é necessário recorrer a métodos numéricos como a função optim, presente na função fitdist do R para estimadores de máxima verossimilhança. Isso acontece também para muitas outras distribuições. O optim aplica métodos numéricos que maximizam funções. Para ilustrar, considere o mesmo exemplo realizado com o método dos momentos. Considere a situação de uma clínica que quer estimar o tempo de atendimento por paciente, desde o momento em que a pessoa chega até a saída, para planejar a alocação de pessoal. Estudos anteriores sugerem que o tempo de atendimento (em horas) é assimétrico à direita — muitos atendimentos são rápidos, mas alguns demoram mais (por exemplo, em casos de crianças, idosos, pessoas com medo de agulha etc). Por isso, consideramos adequado modelar o tempo de atendimento com uma distribuição Gama. Foram coletados 100 valores de tempos de atendimento, em horas. Os dados estão apresentados abaixo:

Tempos = c(
  1.98, 0.60, 1.28, 1.57, 1.71, 1.34, 1.69, 1.14, 1.64, 0.43,
  1.11, 0.92, 2.40, 1.14, 1.32, 1.54, 2.36, 2.56, 3.76, 3.05,
  1.36, 1.88, 2.29, 0.33, 0.92, 1.17, 0.79, 1.49, 2.49, 0.73,
  1.78, 1.30, 1.01, 0.74, 0.88, 0.92, 0.85, 1.01, 1.16, 0.39,
  1.00, 1.90, 0.57, 2.31, 1.62, 1.97, 0.39, 1.67, 0.79, 1.98,
  1.72, 0.47, 1.87, 1.33, 0.90, 3.72, 3.46, 1.91, 2.85, 1.80,
  1.85, 0.55, 0.90, 1.51, 1.57, 2.37, 1.79, 0.48, 0.87, 1.47,
  1.85, 1.24, 0.32, 2.77, 1.69, 0.59, 0.84, 2.36, 2.03, 1.26,
  0.40, 0.63, 0.59, 1.50, 0.52, 0.20, 0.25, 3.09, 3.72, 3.16,
  1.93, 0.85, 1.87, 1.66, 1.31, 1.63, 1.94, 0.20, 2.72, 3.57
)

# Estimadores de Momentos
Alfa_MM

## [1] 3.125629

Beta_MM

## [1] 2.065712

1- Encontre as estimativas de máxima verossimilhança a partir da função automática do R, a fitdist. Guarde estas estimativas nos objetos Alfa_MV e Beta_MV. Compare com as estimativas obtidas pelo método dos momentos.

A título de ilustração, seguem os gráficos das funções de verossimilhança e log-verossimilhança para o caso da distribuição gama (códigos omitidos):

4 DISTRIBUIÇÃO DOS ESTIMADORES

Imagine que você acabou de coletar uma amostra — talvez alturas de pessoas, tempos de espera, medidas de pressão arterial — e precisa estimar algum parâmetro populacional, como uma média, uma proporção, uma variância. Você aplica o método dos momentos ou o método da máxima verossimilhança, e ele te devolve um número, uma estimativa. Ótimo. Mas… e agora? Será que esse número está próximo do verdadeiro valor do parâmetro? Quão confiável ele é? Será que, se eu coletasse outra amostra, obteria um valor muito diferente? E como eu posso expressar matematicamente essa incerteza? É exatamente nesse ponto que entra a distribuição dos estimadores. Ela é a chave que abre a porta da inferência estatística para as propriedades dos estimadores (viés, eficiência e consistência).

Mesmo que a gente tenda a pensar nos estimadores como simples fórmulas ou funções — média amostral, proporção amostral, etc. —, é fundamental lembrar que eles são funções da amostra. E a amostra é aleatória. Portanto, os estimadores também são aleatórios.

Isso quer dizer que:

A cada nova amostra, o valor do estimador muda.
Existe uma distribuição de probabilidades associada a esse estimador.
E essa distribuição tem uma média, uma variância, uma forma.

Ou seja, entender a distribuição do estimador é entender como ele se comporta ao longo de muitas amostras possíveis. É como imaginar um universo paralelo onde você pudesse repetir o experimento milhares de vezes e observar os valores que o estimador toma. Isso é essencial para responder à pergunta: “quão confiável é o meu estimador?”. Dentro dessa distribuição, duas características se destacam: o valor esperado e a variância do estimador.

4.1 ESPERANÇA DO ESTIMADOR

O valor esperado de um estimador é, de forma simplificada, a média de longo prazo dos valores que ele assume se a gente repetir o experimento muitas vezes.

Formalmente, se \(\hat{\theta}\) é o estimador do parâmetro \(\theta\), então o valor esperado de \(\hat{\theta}\) é \(E(\hat{\theta})\). Esse valor esperado nos diz qual é a tendência central da distribuição do estimador e está diretamente relacionada ao seu viés.

EXEMPLO: Considere a coleta da amostra \(X_1, X_2, ...,X_n\) de uma população cuja média é \(\mu\) e a variância é \(\sigma^2\). No caso da média amostral \(\bar{X}\), por exemplo, sabemos que:

\[E(\bar{X})=\mu\]

4.2 VARIÂNCIA DO ESTIMADOR

Já a variância de um estimador mede o quanto ele oscila de uma amostra para outra. É uma medida de precisão: quanto menor a variância, mais estável o estimador, e mais próximos uns dos outros tendem a ser os valores estimados em diferentes amostras. Formalmente, se \(\hat{\theta}\) é o estimador do parâmetro \(\theta\), então a variância de \(\hat{\theta}\) é

\[Var(\hat{\theta})\]

E consequentemente, o erro padrão do estimador \(\hat{\theta}\) é \[EP(\hat{\theta})=\sqrt{Var(\hat{\theta})}\]

Considere a coleta da amostra \(X_1, X_2, ...,X_n\) de uma *população cuja média é \(\mu\) e a variância é \(\sigma^2\). No caso da média amostral \(\bar{X}\), por exemplo, sabemos que:

\[Var(\bar{X})=\frac{\sigma^2}{n}\]

\[EP(\bar{X})=\sqrt{Var(\bar{X})}=\frac{\sigma}{\sqrt{n}}\]

A função fitdist além de soltar as estimativas dos parâmetros, também solta uma estimativa do erro padrão do estimador, tanto para o método dos momentos quanto para o método da máxima verossimilhança. Para o caso da Poisson, o estimador de máxima verossimilhança é a média amostral, assim

# Amostra de uma Poisson
Clientes = c(9, 5, 6, 5, 2, 8, 3, 2, 4, 9, 8, 6, 6, 11, 6, 6, 5, 5, 4, 3, 9, 8, 6, 7, 1, 5, 6, 3, 4, 3, 3, 4, 4, 4, 3, 3)

fitdist(data=Clientes,distr="pois",method="mle")

## Fitting of the distribution ' pois ' by maximum likelihood 
## Parameters:
##        estimate Std. Error
## lambda 5.166667  0.3788384

No caso da distribuição Poisson \(E(\bar{X})=\lambda\) e \(\hat{\lambda}=5,16\) é uma estimativa para \(\lambda\) e, consequentemente, é uma estimativa da esperança do estimador \(\bar{X}\).

Analogamente, no caso da distribuição Poisson, \(Var(\bar{X})=\lambda/n\) e \(EP(\bar{X})=\sqrt{\lambda/n}\). Assim, uma estimativa do erro padrão de \(\bar{X}\) é \(\sqrt{\hat{\lambda}/n}\) = \(\sqrt{5,16/36}\)= 0.378, o valor encontrado na função acima.

De um modo geral, existem fórmulas matemáticas para obter os valores dos erros padrões, algumas são analíticas, ou seja, podem ser obtidas diretamente a partir da teoria das distribuições, e outras são estimativas numéricas, obtidas via métodos computacionais como o bootstrap ou a aproximação da informação de Fisher.

No caso da função fitdist, ela nos fornece uma estimativa do erro padrão do estimador junto com os próprios parâmetros estimados, tanto pelo método dos momentos quanto pelo método da máxima verossimilhança.

5 PROPRIEDADES DOS ESTIMADORES

Tradicionalmente, as propriedades dos estimadores são estudadas de forma teórica: deduzimos a esperança, a variância, o viés, a distribuição assintótica… Mas no contexto da Estatística Computacional, abrimos espaço para outra abordagem: avaliar as propriedades dos estimadores por meio de simulações.

Essa abordagem é valiosa porque:

Permite observar na prática o comportamento de um estimador.
Ajuda a comparar diferentes estimadores mesmo quando as soluções teóricas são complicadas ou não existem.
Serve para verificar empiricamente propriedades assintóticas (por exemplo, como o estimador se comporta quando o tamanho da amostra aumenta).
Cria um caminho para ensinar inferência de maneira mais visual e intuitiva.

A ideia central é simular muitas amostras de uma distribuição conhecida (com parâmetros conhecidos), calcular o estimador em cada amostra e analisar os resultados. Por exemplo, podemos gerar 1.000 amostras de tamanho 30 de uma distribuição normal com média 5, calcular a média amostral em cada uma e, com isso, estimar empiricamente o comportamento desse estimador. Em muitos casos, utilizaremos a função replicate do R para agilizar as replicações da amostra. Essa função tem a seguinte estrutura:

replicate("B réplicas","Função a ser replicada")

As principais propriedades dos estimadores são:

Viés
Consistência
Eficiência
Erro quadrático médio

Mas antes de vermos cada uma dessas propriedades é fundamental que saibamos como estimar a esperança e a variância dos estimadores.

5.1 ESTIMANDO A ESPERANÇA E A VARIÂNCIA DE UM ESTIMADOR

Antes de verificar, com detalhes, cada uma destas propriedades, é importante destacar como estimamos a esperança, variância e erro padrão de um estimador, computacionalmente. Essas estimativas são fundamentais para o bom entendimento das propriedades. Abaixo está o passo a passo computacional para estimar a esperança, variância e erro padrão de um estimador usando simulação:

1- Defina o modelo populacional: escolha uma distribuição populacional conhecida (por exemplo, normal, binomial, exponencial) com parâmetros fixos. Esse será o modelo gerador dos dados.

2- Escolha o estimador \(\hat\theta\): defina claramente qual quantidade populacional \(\theta\) você deseja estimar e qual estimador \(\hat{\theta}\) será utilizado.

3- Fixe o tamanho amostral \(n\): determine quantas observações cada amostra terá.

4- Simule \(B\) amostras independentes: gere \(B\) amostras de tamanho \(n\) da população especificada (em geral, B=1000 ou 10000).

5- Calcule o estimador em cada amostra: para a \(b\)-ésima amostra, calcule o o estimador \(\hat{\theta}^{(b)}\), com \(b=1,2,3,...,B\).

6- Armazene os valores obtidos: guarde os \(B\) valores do estimador em um vetor ou lista: \(\{\hat\theta^{(1)}, \hat\theta^{(2)}, \dots, \hat\theta^{(B)}\}\). Este vetor representa a distribuição empírica do estimador \(\hat{\theta}\).

Com essa distribuição, \(\hat{\theta}\), podemos obter:

A estimativa da esperança de \(\hat{\theta}\):

\[ \widehat{\operatorname{E}}(\hat\theta) = \frac{1}{B} \sum_{b=1}^B \hat\theta^{(b)} \]
A estimativa da variância de \(\hat{\theta}\):

\[ \widehat{\operatorname{Var}}(\hat\theta) = \frac{1}{B - 1} \sum_{b=1}^B (\hat\theta^{(b)} - \widehat{\operatorname{E}}(\hat\theta) )^2 \]
A estimativa do erro padrão de \(\hat{\theta}\):

\[ \widehat{\operatorname{EP}}(\hat\theta) = \sqrt{ \widehat{\operatorname{Var}}(\hat\theta)} \]

5.1.1 EXEMPLO: MÉDIA AMOSTRAL

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população com média \(\mu\) e variância \(\sigma^2\). A média amostral é dada por:

\[ \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i \]

Sabemos que:

\[E(\bar{X})=\mu\] \[Var(\bar{X})=\frac{\sigma^2}{n}\] \[EP(\bar{X})=\frac{\sigma}{\sqrt{n}}\]

Vamos conduzir um estudo assumindo que \(X_1, X_2, \dots, X_n\) é uma amostra aleatória de uma população normalmente distribuída com média \(\mu\)=170 e variância \(\sigma^2=100\). Nosso objetivo será estimar empiricamente a esperança, variância e erro padrão de \(\bar{X}\). Para isso, vamos extrair \(B=1000\) amostras de tamanho \(n=25\) dessa população. Seguem os comandos:

# Parâmetros populacionais
media_real <- 170
sigma_real <- 10

# Quantidade de réplicas
B <- 1000

# Tamanho de amostra
n = 25

# Guarda as médias das B réplicas (distribuição do estimador)
set.seed(2025)
medias <- replicate(B, mean(rnorm(n, mean = media_real, sd = sigma_real)))

# Estimativa da Esperança
Est_Esperanca = mean(medias)

# Estimativa da variância
Est_Variância <- var(medias)

# Estimativa do erro padrão
Est_ErroPadrão <- sd(medias)

# Data frame para o ggplot
df <- data.frame(media = medias)

# Gráfico com ggplot2
ggplot(df, aes(x = media)) +
  # Histograma
  geom_histogram(fill = "lightblue", color = "black") +
  labs(title = "Histograma das medias amostrais",
       x = "Media amostral", y = "Frequência") +
  theme_minimal()

Assim, temos as seguintes estimativas: \(\widehat{\operatorname{E}}(\hat\theta)=169,98\), \(\widehat{\operatorname{Var}}(\hat\theta)=3,94\) e \(\widehat{\operatorname{EP}}(\hat\theta)=1,98\). Os valores reais são \(E(\hat\theta)=\mu=170\), \(Var(\hat\theta)=\frac{\sigma^2}{n}=\frac{100}{25}=4\) e \(EP(\hat\theta)=\frac{\sigma}{\sqrt{n}}=\frac{10}{5}=2\). Tem-se que os resultados teóricos são validados computacionalmente.

5.1.2 EXERCÍCIO: PROPORÇÃO AMOSTRAL

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população com proporção populacional \(p\), onde cada \(X_i\) segue uma distribuição Bernoulli com parâmetro \(p\) (isto é, \(X_i \sim \text{Bernoulli}(p)\)).

A proporção amostral é dada por:

\[ \hat{p} = \frac{1}{n} \sum_{i=1}^{n} X_i \]

Sabemos que:

\[ E(\hat{p}) = p \]

\[ Var(\hat{p}) = \frac{p(1-p)}{n} \]

\[ EP(\hat{p}) = \sqrt{\frac{p(1-p)}{n}} \]

Conduza um estudo assumindo que \(X_1, X_2, \dots, X_n\) é uma amostra aleatória de uma população onde \(p = 0{,}6\) e coletou-se uma amostra de tamanho \(n=25\). Nosso objetivo será estimar empiricamente a esperança, variância e erro padrão de \(\hat{p}\) dados por:

Proporção populacional:
\[ p = 0,6 \]
Esperança de \(\hat{p}\):
\[ E(\hat{p}) = p = 0,6 \]
Variância de \(\hat{p}\):
\[ Var(\hat{p}) = \frac{p(1-p)}{n} = \frac{0,6 \cdot 0,4}{25} = \frac{0,24}{25} = 0,0096 \]
Erro padrão de \(\hat{p}\):
\[ EP(\hat{p}) = \sqrt{0,0096} \approx 0,098 \]

Para isso, extraia \(B = 1000\) amostras de tamanho \(n = 25\) dessa população.
Faça o histograma da distribuição deste estimador.

5.2 VIÉS

O viés (bias) de um estimador \(\hat{\theta}\) é definido como:

\[ \text{Viés}(\hat{\theta}) = \mathbb{E}[\hat{\theta}] - \theta \]

onde \(\mathbb{E}[\hat{\theta}]\) representa o valor esperado do estimador \(\hat{\theta}\).

Se \(\text{Viés}(\hat{\theta}) = 0\), dizemos que o estimador é não viesado
Se \(\text{Viés}(\hat{\theta}) \ne 0\), o estimador é viesado.

No R, podemos estimar o viés por meio de simulação computacional. Assim, a estimativa do viés do estimador é dada por:

\[\widehat{\text{Viés}}(\hat{\theta}) = \widehat{\operatorname{E}}(\hat\theta) - \theta\]

5.2.1 EXEMPLO: VIÉS DA MÉDIA AMOSTRAL

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população com média \(\mu\) e variância \(\sigma^2\). A média amostral é dada por:

\[ \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i \]

Queremos verificar se a média amostral \(\bar{X}\) é viesada como estimador de \(\mu\).

Sabemos que, teoricamente:

\[ \mathbb{E}[\bar{X}] = \mu \quad \Rightarrow \quad \text{Viés}(\bar{X}) = \mathbb{E}[\bar{X}] - \mu = 0 \]

Ou seja, a média amostral é um estimador não-viesado da média populacional.

Vamos conduzir um estudo computacional assumindo que a população é descrita por uma variável Normal(\(\mu=170\),\(\sigma^2=100\)) e estaremos coletando amostras de tamanho \(n=30\):

# Parâmetros populacionais
media_real <- 170
sigma_real <- 10

# Quantidade de réplicas
B <- 1000

# Tamanho de amostra
n = 30

# Guarda as médias das B réplicas (distribuição do estimador)
set.seed(2025)
medias <- replicate(B, mean(rnorm(n, mean = media_real, sd = sigma_real)))

# Estimativa da Esperança
Est_Esperanca = mean(medias)

# Estimativa do viés
Est_Vies <- Est_Esperanca-media_real

# Data frame para o ggplot
df <- data.frame(media = medias)

# Gráfico com ggplot2
ggplot(df, aes(x = media)) +
  # Histograma
  geom_histogram(fill = "lightblue", color = "black") +
  # Marcando de vermelho o verdadeiro valor do parâmetro.
  geom_vline(xintercept = media_real, color = "red", linewidth = 1) +
  # Marcando de azul a estimativa da esperança do estimador
  geom_vline(xintercept = Est_Esperanca, color = "blue", linewidth = 1) +
  labs(title = "Histogramas das medias amostrais",
       x = "Media amostral", y = "Frequência") +
  theme_minimal()

Veja que o viés estimado é muito próximo de 0. Comprova-se, portanto, a falta de viés da média amostral.

5.2.2 EXEMPLO: VIÉS DE ESTIMADORES DE \(\theta\) DA UNIFORME(0,\(\theta\))

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população distribuição Uniforme (0,\(\theta\)). Para estimar \(\theta\), temos dois estimadores, por exemplo:

Estimador de momentos: \(\hat{\theta}_{MM}=2\bar{X}\)
Estimador de máxima verossimilhança: \(\hat{\theta}_{MV} = X_{(n)}\)

Pode-se provar que \(E(\hat{\theta}_{MM})=\theta\) e \(E(\hat{\theta}_{MV})=\frac{n}{n+1}\theta\). Assim, \(\hat{\theta}_{MM}\) é um estimador não viesado pra \(\theta\) e \(\hat{\theta}_{MV}\) é um estimador viesado pra \(\theta\) que, particularmente, subestima o parâmetro.

Vamos conduzir um estudo computacional coletando B=1000 amostras de tamanho \(n=30\) de uma população uniformemente distribuida com \(\theta=100\). Nesse contexto:

Para o estimador dos momentos: \[ E(\hat{\theta}_{MM}) = \theta = 100 \] \[ Viés(\hat{\theta}_{MM}) = E(\hat{\theta}_{MM}) - \theta = 0 \]
Para o estimador de máxima verossimilhança: \[ E(\hat{\theta}_{MV}) = \frac{n}{n+1}\theta = \frac{30}{31} \cdot 100 \approx 96,77 \] \[ Viés(\hat{\theta}_{MV}) = E(\hat{\theta}_{MV}) - \theta \approx 96,77 - 100 = -3,23 \]

Portanto:

\(\hat{\theta}_{MM}\) é não viesado.
\(\hat{\theta}_{MV}\) apresenta viés negativo, subestimando \(\theta\) em média por cerca de 3,23 unidades quando \(n=30\).

Segue o script onde vamos estimar estes resultados teóricos:

# Parâmetro populacional
theta_real <- 100

# Quantidade de réplicas
B <- 10000

# Tamanho de amostra
n = 30

# Guarda as estimativas de momentos das B réplicas 
set.seed(2025)
theta_mm <- replicate(B, 2*mean(runif(n, max=theta_real)))

# Guarda as estimativas de máxima verossimilhança das B réplicas 
set.seed(2025)
theta_mv <- replicate(B, max(runif(n, max=theta_real)))

# Estimativa da esperança dos estimadores
Est_Esperanca_mm = mean(theta_mm)
Est_Esperanca_mv = mean(theta_mv)

# Estimativa do viés
Est_Vies_mm <- Est_Esperanca_mm-theta_real
Est_Vies_mv <- Est_Esperanca_mv-theta_real

# Data frames para o ggplot
df_mm <- data.frame(theta_mm = theta_mm)
df_mv <- data.frame(theta_mv = theta_mv)

# Gráficos com ggplot2
A = ggplot(df_mm, aes(x = theta_mm)) +
      geom_histogram(fill = "lightblue", color = "black") +
      geom_vline(xintercept = theta_real, color = "red", linewidth = 1) +
      geom_vline(xintercept = Est_Esperanca_mm, color = "blue", linewidth = 1) +
      labs(title = "Histograma do estimador MM",
               x = "Estimador MM", y = "Frequência") +
      theme_minimal()+
      xlim(60,130)

B = ggplot(df_mv, aes(x = theta_mv)) +
      geom_histogram(fill = "lightblue", color = "black") +
      geom_vline(xintercept = theta_real, color = "red", linewidth = 1) +
      geom_vline(xintercept = Est_Esperanca_mv, color = "blue", linewidth = 1) +
      labs(title = "Histograma do estimador MV",
               x = "Estimador MV", y = "Frequência") +
      theme_minimal()+
      xlim(60,130)
A/B

Assim, comprovamos que o estimador de MV é viesado e tende a subestimar o parâmetro \(\theta\), especialmente, quando \(n\) não é suficientemente grande.

5.2.3 EXERCÍCIO: VARIÂNCIAS AMOSTRAIS

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população com média \(\mu\) e variância \(\sigma^2\). Temos dois estimadores para a variância \(\sigma^2\), dados por:

\[S^2_{\text{V}} = \frac{1}{n} \sum_{i=1}^n (X_i - \bar{X})^2\]

\[S^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2\]

Pode-se provar que \(E(S^2_V) = \sigma^2 \cdot \frac{n-1}{n}\) e \(E(S^2) = \sigma^2\). Portanto \(S^2\) é um estimador não viesado para \(\sigma^2\) (Viés = 0) e o estimador \(S^2_v\) é viesado para \(\sigma^2\) (Viés = \(-\frac{\sigma^2}{n}\)).

Considere uma população normalmente distribuída com média \(\mu=170\) e variância \(\sigma^2=100\). Para a coleta de amostras de tamanho \(n=30\), faça, em sequencia:

Calcule o viés teórico de ambos os estimadores.
Gere 1000 amostras de tamanho n=30 dessa população. Calcule, para cada uma, os estimadores \(S^2_{\text{V}}\) e \(S^2\). Calcule as estimativas das esperanças dos dois estimadores e os respectivos vieses estimados. Faça um gráfico para cada estimador com os resultados.

5.3 CONSISTÊNCIA

A consistência de um estimador é uma propriedade fundamental que garante que, à medida que o tamanho da amostra aumenta, o estimador converge para o valor verdadeiro do parâmetro que está sendo estimado. Utilizaremos agora a notação \(\hat{\theta}_n\) para o estimador apenas para reforçar que ele depende do tamanho amostral \(n\). Formalmente, um estimador \(\hat{\theta}_n\) é consistente para o parâmetro \(\theta\) se:

\[ \hat{\theta}_n \xrightarrow{P} \theta \quad \text{quando} \quad n \to \infty, \]

onde \(\xrightarrow{P}\) indica convergência em probabilidade, ou seja, para qualquer \(\epsilon > 0\),

\[ \lim_{n \to \infty} P(|\hat{\theta}_n - \theta| \geq \epsilon) = 0. \]

Isso significa que, com o aumento do tamanho da amostra \(n\), a probabilidade de o estimador \(\hat{\theta}_n\) se desviar de \(\theta\) por mais de \(\epsilon\) vai para zero.

Para verificar computacionalmente se um estimador é consistente, podemos usar simulações para observar seu comportamento à medida que o tamanho da amostra aumenta. O processo é o seguinte:

Para verificar computacionalmente se um estimador é consistente, podemos usar simulações para observar seu comportamento à medida que o tamanho da amostra aumenta. O processo é o seguinte:

Escolha um parâmetro \(\theta\)
Defina o parâmetro para o qual você deseja verificar a consistência do estimador.
Gere amostras
Para diferentes tamanhos de amostras \(n\) (por exemplo, \(n = 10, 30, 50, 100, 500, \dots\)), gere \(B\) amostras da população indexada por \(\theta\).
Aplique o estimador
Para cada amostra, calcule o estimador \(\hat{\theta}_n\).
Estime média e variância do estimador
Para cada tamanho de amostra \(n\), estime:
- A esperança de \(\hat{\theta}_n\).
- A variância de \(\hat{\theta}_n\).
Analise a convergência
- Verifique se, à medida que \(n\) aumenta, a esperança estimada de \(\hat{\theta}_n\) se aproxima de \(\theta\).
- Verifique se a variância estimada de \(\hat{\theta}_n\) diminui com \(n\), tendendo a 0.

Se essas duas condições forem atendidas, temos evidências de que o estimador é consistente.

5.3.1 EXEMPLO: UNIFORME (0,\(\theta\))

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população descrita por uma variável Uniforme \((0,\theta)\).

O valor máximo da amostra, dado por

\[ X_{(n)} = \max(X_1,\dots,X_n), \]

é o estimador de máxima verossimilhança de \(\theta\).

Sabe-se que:

Esperança do estimador:
\[ \mathbb{E}[X_{(n)}] = \frac{n}{n+1}\theta \]
Variância do estimador:
\[ \mathrm{Var}(X_{(n)}) = \frac{n}{(n+1)^2(n+2)}\theta^2 \]

📌 Isso mostra que o estimador é viesado para baixo, mas:
- \(\mathbb{E}[X_{(n)}] \to \theta\) quando \(n \to \infty\)
- \(\mathrm{Var}(X_{(n)}) \to 0\) quando \(n \to \infty\)

Ou seja, o estimador é consistente. Segue o script abaixo que ilustra isso:

# Parâmetro populacional
theta_real = 50

# Quantidade de réplicas
B <- 1000

# Função para gerar o gráfico
gerar_grafico <- function(n) {

  set.seed(2025)
  thetas <- replicate(B, max(runif(n, min=0, max=theta_real)))
  
  # Estimativa da esperança e da variância
  Est_Esperanca <- mean(thetas)
  
  # Data frame para o ggplot
  df <- data.frame(thetas = thetas)
  
  # Gráfico com ggplot2
  ggplot(df, aes(x = thetas)) +
    geom_histogram(fill = "lightblue", color = "black") +
    geom_vline(xintercept = theta_real, color = "red", linewidth = 1) +
    geom_vline(xintercept = Est_Esperanca, color = "blue", linewidth = 1) +
    labs(title = paste("n=", n, sep = ""),
         x = "Máximo Amostral", y = "Frequência absoluta") +
    theme_minimal() +
    xlim(25,50)
}

n <- c(10, 30, 50, 100)
graficos <- lapply(n, gerar_grafico)

# Combinar os gráficos com patchwork
graficos[[1]] + graficos[[2]] + graficos[[3]] + graficos[[4]] + plot_layout(ncol = 2)

Note como a estimativa da esperança do estimador vai se aproximando do verdadeiro valor do parâmetro e a variância diminui à medida que o tamanho da amostra aumenta refletindo na consistência deste estimador.

5.3.2 EXERCÍCIO: VARIÂNCIA AMOSTRAL VIESADA

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população descrita por uma variável Normal com parâmetros \((\mu,\sigma^2)\).

Nosso interesse está na estimação do parâmetro \(\sigma^2\). Considere o seguinte estimador:

\[ S_v^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2 \]

Este é o chamado estimador da variância amostral viesada, pois:

\[ \mathbb{E}[S_v^2] = \frac{n-1}{n}\sigma^2 \neq \sigma^2 \]

Ou seja, subestima a variância verdadeira.

No entanto, apesar de viesado, pode-se provar que ele é consistente. Verifique isso computacionalmente:

Crie uma função chamada ConsistenciaVarViesada que retornará as estimativas da esperança e da variância do estimador \(S_v^2\), com base em \(B=1000\) réplicas. Essa função deverá ter como argumentos \(\mu\), \(\sigma\) e \(n\).
Suponha que a população possui distribuição \(N(\mu=170,\sigma^2=100)\). Avalie a consistência do estimador com a função criada observando o que ocorre quando \(n=10, 50, 500, 1000\).

5.4 EFICIÊNCIA

Na inferência estatística, frequentemente nos deparamos com mais de um estimador não viesado para o mesmo parâmetro populacional. Para decidir qual estimador utilizar, comparamos a variância desses estimadores. Surge então o conceito de eficiência.

Sejam \(\hat{\theta}_1\) e \(\hat{\theta}_2\) dois estimadores não viesados de \(\theta\). Dizemos que \(\hat{\theta}_1\) é mais eficiente que \(\hat{\theta}_2\) se:

\[ \operatorname{Var}(\hat{\theta}_1) < \operatorname{Var}(\hat{\theta}_2) \]

Ou seja, entre dois estimadores não viesados, o mais eficiente é aquele com menor variância. A eficiência de \(\hat{\theta}_1\) em relação a \(\hat{\theta}_2\) é definida como:

\[ \text{Eficiência}_{var}(\hat{\theta}_1,\hat{\theta}_2) = \frac{\operatorname{Var}(\hat{\theta}_1)}{\operatorname{Var}(\hat{\theta}_2)} \]

Se a eficiência for maior que 1, então \(\hat{\theta}_2\) é mais eficiente.
Se for menor que 1, então \(\hat{\theta}_1\) é mais eficiente.
Se for igual a 1, os dois estimadores têm a mesma eficiência.

No R, para estimar a eficiência de dois estimadores não viesados \(\hat{\theta_1}\) e \(\hat{\theta_2}\), basta fazer:

\[ \widehat{\text{Eficiência}_{var}}(\hat{\theta}_1,\hat{\theta}_2) = \frac{\widehat{\operatorname{Var}}(\hat{\theta}_1)}{\widehat{\operatorname{Var}}(\hat{\theta}_2)} \]

5.4.1 EXEMPLO: MÉDIA E MEDIANA AMOSTRAIS

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população com distribuição normal de média \(\mu = 170\) e variância \(\sigma^2 = 100\).

Nosso objetivo é comparar dois estimadores que podem ser usadas para estimar \(\mu\):

A média amostral: \(\bar{X} = \frac{1}{n} \sum_{i=1}^n X_i\).
A mediana amostral: \(\tilde{X}\) = valor central da amostra ordenada.

Sabemos que:

A média amostral é um estimador não viesado. Isto é, \(E(\bar{X})=\mu\).
A mediana amostral é um estimador não viesado. Isto é, \(E(\tilde{X})=\mu\).

Vamos comparar esses dois estimadores quanto à sua eficiência, definida por:

\[ \text{Eficiência}_{Var}(\bar{X},\tilde{X}) = \frac{\text{Var}(\bar{X})}{\text{Var}(\tilde{X})} \]

Segue o código que estima a eficiência para \(n=25\):

# Parâmetros populacionais
media_real <- 170
sigma_real <- 10

# Quantidade de réplicas
B <- 1000
# Tamanho de amostra
n = 25

# Guarda as médias das B réplicas (distribuição do estimador)
set.seed(2025)
medias <- replicate(B, mean(rnorm(n, mean = media_real, sd = sigma_real)))

# Guarda as medianas das B réplicas (distribuição do estimador)
set.seed(2025)
medianas <- replicate(B, median(rnorm(n, mean = media_real, sd = sigma_real)))

# Estimativa da Variância do estimador média amostral
Est_Variancia_media = var(medias)
# Estimativa da Variância do estimador mediana amostral
Est_Variancia_mediana = var(medianas)

# Estimativa da eficiência:
Est_Eficiencia = Est_Variancia_media/Est_Variancia_mediana

# Resultados:
Est_Variancia_mediana

## [1] 6.294808

Est_Variancia_media

## [1] 3.946686

Est_Eficiencia

## [1] 0.6269747

Agora vamos fazer um estudo variando o tamanho da amostra \(n\):

# Parâmetros populacionais
media_real <- 170
sigma_real <- 10

# Quantidade de réplicas
B <- 1000

# Tamanhos de amostra
n_vals <- c(10, 20, 30, 50, 100, 200, 300, 400, 500)

# Vetor para armazenar as eficiências
eficiencias <- numeric(length(n_vals))

# Loop para calcular a eficiência para cada tamanho de amostra
for (i in 1:length(n_vals)) {
  n <- n_vals[i]
  
  # Guarda as médias das B réplicas (distribuição do estimador)
  set.seed(2025)
  medias <- replicate(B, mean(rnorm(n, mean = media_real, sd = sigma_real)))
  
  # Guarda as medianas das B réplicas (distribuição do estimador)
  set.seed(2025)
  medianas <- replicate(B, median(rnorm(n, mean = media_real, sd = sigma_real)))
  
  # Estimativa da Variância do estimador média amostral
  Est_Variancia_media = var(medias)
  # Estimativa da Variância do estimador mediana amostral
  Est_Variancia_mediana = var(medianas)
  
  # Estimativa da eficiência:
  Est_Eficiencia = Est_Variancia_media / Est_Variancia_mediana
  
  # Armazenando a eficiência
  eficiencias[i] <- Est_Eficiencia
}

# Criar um gráfico da eficiência em função de n
dados <- data.frame(
  n = n_vals,
  eficiencia = eficiencias
)

# Gráfico
ggplot(dados, aes(x = n, y = eficiencia)) +
  geom_line(color = "blue") +
  geom_point(color = "blue") +
  labs(
    title = "Eficiência da Media em relação à Mediana para Diferentes Tamanhos de Amostra",
    x = "Tamanho da Amostra (n)",
    y = "Eficiência"
  ) +
  theme_minimal()

Como esperado, a média amostral é mais eficiente.

5.4.2 EXERCÍCIO: DISTRIBUIÇÃO EXPONENCIAL

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória simples de uma distribuição exponencial \(\text{Exp}(\lambda)\).

Dois estimadores não viesados são propostos para estimar a média \(\mu = \frac{1}{\lambda}\):

Estimador 1 (média amostral):

\[ \hat{\mu}_1 = \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i \]

A média amostral é um estimador não tendencioso para \(\mu\). A variância da média amostral \(\bar{X}\) é dada por:

\[ \text{Var}(\hat{\mu}_1) = \frac{\sigma^2}{n} \]

Para a distribuição exponencial \(\text{Exp}(\lambda)\), temos que \(\sigma^2 = \frac{1}{\lambda^2}\), logo:

\[ \text{Var}(\hat{\mu}_1) = \frac{1}{n\lambda^2} \]

Estimador 2 (mínimo ajustado)

\[ \hat{\mu}_2 = n \times \min(X_1, \dots, X_n) \]

Pode-se provar que:

\[ \mathbb{E}[\hat{\mu}_2] = n \times \mathbb{E}[\min(X_1, \dots, X_n)] = n \times \frac{1}{n\lambda} = \frac{1}{\lambda} = \mu \]

Logo, \(\hat{\mu}_2\) também é um estimador não tendencioso para \(\mu\).

A variância de \(\min(X_1, \dots, X_n)\) é:

\[ \text{Var}(\min(X_1, \dots, X_n)) = \frac{1}{(n\lambda)^2} \]

Portanto:

\[ \text{Var}(\hat{\mu}_2) = n^2 \times \text{Var}(\min(X_1, \dots, X_n)) = n^2 \times \frac{1}{n^2\lambda^2} = \frac{1}{\lambda^2} \]

A eficiência teórica do estimador 2 em relação ao estimador 1 é dada pela razão das variâncias:

\[ \text{Eficiência}_{var}(\hat{\mu}_1,\hat{\mu}_2) = \frac{\text{Var}(\hat{\mu}_1)}{\text{Var}(\hat{\mu}_2)} = \frac{\frac{1}{n\lambda^2}}{\frac{1}{\lambda^2}} = \frac{1}{n} \]

Portanto, a eficiência teórica do estimador 2 em relação ao estimador 1 é:

\[ \boxed{\text{Eficiência}_{var} = \frac{1}{n}} \] 1- Considere \(\lambda=1/50\), \(n=25\) e \(B=1000\). Estime, empiricamente, a eficiência entre estes dois estimadores.

2- Faça um estudo variando o tamanho da amostra \(n\) com a eficiência, considerando n = 10, 20, 30, 50, 100, 200, 300, 400 e 500.

5.5 ERRO QUADRÁTICO MÉDIO

Na prática, um estimador pode ser não viesado, mas ainda assim apresentar alta variância, tornando-o indesejável. Para avaliar a qualidade global de um estimador, consideramos o seu Erro Quadrático Médio (EQM).

Dado um estimador \(\hat{\theta}\) de um parâmetro \(\theta\), o EQM é definido como:

\[ \operatorname{EQM}(\hat{\theta}) = \mathbb{E}[(\hat{\theta} - \theta)^2] \]

O EQM pode ser decomposto em duas partes fundamentais:

\[ \operatorname{EQM}(\hat{\theta}) = \operatorname{Var}(\hat{\theta}) + \left(\operatorname{Viés}(\hat{\theta})\right)^2 \]

Ou seja, o erro quadrático médio incorpora tanto a variância quanto o viés do estimador. Isso significa que mesmo um estimador viesado pode ter baixo EQM, caso seu viés seja pequeno e sua variância também.

Dado dois estimadores \(\hat{\theta}_1\) e \(\hat{\theta}_2\), podemos dizer que:

\(\hat{\theta}_1\) é melhor que \(\hat{\theta}_2\) se \(\operatorname{EQM}(\hat{\theta}_1) < \operatorname{EQM}(\hat{\theta}_2)\).

Essa comparação leva em conta tanto a precisão (variância) quanto a exatidão (viés) dos estimadores. Podemos definir, ainda, a eficiência de dois estimadores \(\hat\theta_1\) e \(\hat\theta_2\) que é a razão entre seus erros quadráticos médios (EQMs). Suponha que ambos sejam estimadores de um mesmo parâmetro \(\theta\). A eficiência de \(\hat\theta_1\) em relação a \(\hat\theta_2\) é dada por:

\[ \text{Eficiência}_{EQM}(\hat{\theta_1},\hat{\theta_2}) = \frac{\text{EQM}(\hat\theta_1)}{\text{EQM}(\hat\theta_2)} \]

Uma eficiência menor que 1 indica que \(\hat\theta_1\) é mais eficiente que \(\hat\theta_2\), ou seja, possui menor EQM.

No R, para estimar o erro quadrático médio de um estimador \(\hat{\theta}\), basta fazer:

\[\widehat{EQM}(\hat{\theta})=\widehat{\operatorname{Var}}(\hat{\theta}) + \left(\widehat{\operatorname{Viés}}(\hat{\theta})\right)^2\]

Consequentemente, para estimar a eficiência entre dois estimadores \(\hat\theta_1\) e \(\hat\theta_2\), tem-se:

\[\widehat{\text{Eficiência}_{EQM}}(\hat{\theta_1},\hat{\theta_2}) = \frac{\widehat{EQM}(\hat{\theta_1})}{\widehat{EQM}(\hat{\theta_2})}\]

5.5.1 EXEMPLO: DISTRIBUIÇÃO UNIFORME

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população distribuição Uniforme (0,\(\theta\)). Para estimar \(\theta\), temos dois estimadores, por exemplo:

Estimador de momentos: \(\hat{\theta}_{MM}=2\bar{X}\)
Estimador de máxima verossimilhança: \(\hat{\theta}_{MV} = X_{(n)}\)

Pode-se provar que, analiticamente:

\[EQM(\hat{\theta}_{MM}) = \frac{\theta^2}{3n}\]

\[EQM(\hat{\theta}_{MV}) = \frac{2\theta^2}{(n+1)(n+2)}\]

A eficiência do estimador de máxima verossimilhança em relação ao estimador de momentos é dada por:

\[ \text{Eficiência}_{EQM}(\hat{\theta}_{MM},\hat{\theta}_{MV}) = \frac{EQM(\hat{\theta}_{MM})}{EQM(\hat{\theta}_{MV})} = \frac{\dfrac{\theta^2}{3n}}{\dfrac{2\theta^2}{(n+1)(n+2)}} = \frac{(n+1)(n+2)}{6n} \]

Assim, quanto maior o tamanho amostral \(n\), mais eficiente se torna o estimador de máxima verossimilhança em comparação com o de momentos.

Vamos avaliar o erro quadrático médio de ambos os estimadores e a eficiencia relativa entre eles coletando amostras de tamanho \(n=30\) de uma população uniformemente distribuida com \(\theta=100\). Para isso, os resultados teóricos são:

\[EQM(\hat{\theta}_{MM}) = \frac{100^2}{3*30}=111,11\] \[EQM(\hat{\theta}_{MV}) = \frac{2*100^2}{(30+1)(30+2)}=20,16\] \[\text{Eficiência}_{EQM} = \frac{111,11}{20,16} = 5,51\]

Seguem os comandos que estimam, empiricamente, estes erros quadráticos médios e eficiência.

# Parâmetro populacional
theta_real <- 100

# Quantidade de réplicas
B <- 1000

# Tamanho de amostra
n = 30

# Guarda as estimativas de momentos das B réplicas 
set.seed(2025)
theta_mm <- replicate(B, 2*mean(runif(n, max=theta_real)))

# Guarda as estimativas de máxima verossimilhança das B réplicas 
set.seed(2025)
theta_mv <- replicate(B, max(runif(n, max=theta_real)))

# Estimativa da esperança dos estimadores
Est_Esperanca_mm = mean(theta_mm)
Est_Esperanca_mv = mean(theta_mv)

# Estimativa do viés dos estimadores
Est_Vies_mm <- Est_Esperanca_mm-theta_real
Est_Vies_mv <- Est_Esperanca_mv-theta_real

# Estimativa da variância dos estimadores
Est_Var_mm <- var(theta_mm)
Est_Var_mv <- var(theta_mv)

# Estimativa do erro quadrático médio dos estimadores
Est_EQM_mm <- Est_Var_mm + Est_Vies_mm^2
Est_EQM_mv <- Est_Var_mv + Est_Vies_mv^2

# Estimativa da eficiência relativa entre os dois estimadores
Est_Eficiencia <- Est_EQM_mm/Est_EQM_mv

Est_EQM_mm

## [1] 106.9369

Est_EQM_mv

## [1] 19.05424

Est_Eficiencia

## [1] 5.612236

Conforme esperado, as estimativas são próximas dos valores reais. Dessa forma, apesar do estimador de momentos ser não viesado, ele apresenta um erro quadrático médio bem inferior ao do estimador de máxima verossimilhança e, consequentemente, é menos eficiente. Agora vamos verificar como essa eficiência varia de acordo com o tamanho da amostra \(n\):

# Parâmetro populacional
theta_real <- 100

# Quantidade de réplicas
B <- 1000

# Tamanhos de amostra
n_vals <- c(10, 20, 30, 50, 100, 200, 300, 400, 500)

# Vetor de eficiências
eficiencias <- numeric(length(n_vals))

# Loop para calcular a eficiência para cada tamanho de amostra
for (i in 1:length(n_vals)) {
  n <- n_vals[i]
  # Guarda as estimativas de momentos das B réplicas 
  set.seed(2025)
  theta_mm <- replicate(B, 2*mean(runif(n, max=theta_real)))
  # Guarda as estimativas de máxima verossimilhança das B réplicas 
  set.seed(2025)
  theta_mv <- replicate(B, max(runif(n, max=theta_real)))
  # Estimativa da esperança dos estimadores
  Est_Esperanca_mm = mean(theta_mm)
  Est_Esperanca_mv = mean(theta_mv)

  # Estimativa do viés dos estimadores
  Est_Vies_mm <- Est_Esperanca_mm-theta_real
  Est_Vies_mv <- Est_Esperanca_mv-theta_real

  # Estimativa da variância dos estimadores
  Est_Var_mm <- var(theta_mm)
  Est_Var_mv <- var(theta_mv)

  # Estimativa do erro quadrático médio dos estimadores
  Est_EQM_mm <- Est_Var_mm + Est_Vies_mm^2
  Est_EQM_mv <- Est_Var_mv + Est_Vies_mv^2

  # Estimativa da eficiência relativa entre os dois estimadores
  Est_Eficiencia <- Est_EQM_mm/Est_EQM_mv
  
  # Armazenando a eficiência
  eficiencias[i] <- Est_Eficiencia
}

  # Criar um gráfico da eficiência em função de n
  dados <- data.frame(
  n = n_vals,
  eficiencia = eficiencias
 )

  # Gráfico
  ggplot(dados, aes(x = n, y = eficiencia)) +
   geom_line(color = "blue") +
   geom_point(color = "blue") +
   labs(
    title = "Eficiência do estimador mm em relação ao estimador mv para Diferentes Tamanhos de Amostra",
    x = "Tamanho da Amostra (n)",
    y = "Eficiência"
  ) +
  theme_minimal()

5.5.2 EXERCÍCIO: VARIÂNCIAS AMOSTRAIS

Seja \(X_1, X_2, \dots, X_n\) uma amostra aleatória de uma população com distribuição Normal(\(\mu\), \(\sigma^2\)). Para estimar \(\sigma^2\), temos dois estimadores, por exemplo:

Estimador não viesado: \(S^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2\)
Estimador viesado: \(S^2_v = \frac{1}{n} \sum_{i=1}^n (X_i - \bar{X})^2\)

Pode-se provar que, analiticamente:

\[ EQM(S^2) = \frac{2}{n-1} \sigma^4 \]

\[ EQM(S^2_v) = \frac{2n-1}{n^2} \sigma^4 \]

A eficiência do estimador não viesado em relação ao viesado é dada por:

\[ \text{Eficiência}_{EQM}(S^2,S^2_v) = \frac{EQM(S^2)}{EQM(S^2_v)} = \frac{\dfrac{2}{n-1}\sigma^4}{\dfrac{2n-1}{n^2}\sigma^4} = \frac{2n^2}{(n-1)(2n-1)} \]

Assim, para grandes amostras (\(n\) grande), a eficiência de ambos os estimadores se torna muito próxima, em termos do MSE.

1- Avalie o erro quadrático médio de ambos os estimadores e a eficiência entre eles coletando B=1000 amostras de tamanho \(n=30\) de uma população normalmente distribuída com média \(\mu=170\) e variância \(\sigma^2=100\). Para isso, os resultados teóricos são:

\[ EQM(S^2) = \frac{2}{30-1} \cdot 100^2 = \frac{2}{29} \cdot 10000 = 689.66 \]

\[ EQM(S^2_v) = \frac{2*30-1}{30^2} \cdot 100^2 = \frac{59}{900} \cdot 10000 = 655.56 \]

\[ \text{Eficiência}_{EQM} = \frac{689.66}{655.56} \approx 1,05 \]

Estime estes resultados a partir de simulação.

2- Faça um estudo variando o tamanho da amostra \(n\) com a eficiência, considerando n = 10, 20, 30, 50, 100, 200, 300, 400 e 500.