Contexto

Como proposta da disciplina CE087 (Estatística Computacional II), vamos comparar as propriedades de três estimadores da média de uma população, a média aritmética, a trimédia e o mid-range. Para isso, decidimos trabalhar com 6 diferentes cenários, buscando explorar o desempenho de cada estimador diante de diferentes situações.

Introdução

Podemos definir como um bom estimador aquele que possui menor viés, menor variância, e maior constância na estimação do parâmetro em questão. A seguir, examinaremos abordagens analíticas para avaliar esses critérios.

• Vício

Por definição, o vício (bias) de um estimador \(\theta\) pode ser calculado como:

\[B(\hat\theta)= E(\hat\theta)-\theta\]

Assim, conforme \(E(\hat{\theta})\) se aproxima do verdadeiro valor de \(\theta\), \(B(\hat{\theta})\) se aproxima de \(0\), configurando um estimador não viciado (esperança do estimador é igual ao valor do parâmetro).

• Variância

Escolher um estimador com menor variância proporciona estimativas mais precisas e consistentes do parâmetro de interesse, reduzindo a dispersão dos valores estimados. Além da variância, podemos calcular o erro padrão do estimador, importante para testes de hipóteses e intervalos de confiança, da seguinte forma:

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

• Erro Quadrático Médio (EQM)

O erro quadrático médio é uma medida que equilibra o vício e a variância de um estimador. Se o estimador for não viciado, seu EQM é igual à sua variância. Seu cálculo é feito da seguinte forma:

\[EQM = V(\hat\theta)+[B(\hat\theta)]^2\]

Assim, valores mais baixos indicam estimativas mais precisas.

• Eficiência Relativa

A eficiência relativa entre dois estimadores é definida como:

\[Efr(\hat\theta_1,\hat\theta_2) = \frac{EQM(\hat\theta_1)}{EQM(\hat\theta_2)}\]

Dessa forma, se \(Efr(\hat\theta_1,\hat\theta_2) < 1\), concluimos que \(\hat\theta_1\) fornece estimativas mais precisas para o parâmetro em consideração (e vice-versa).

Embora tenhamos listado formas analíticas de avaliar o desempenho de um estimador neste tópico, nosso trabalho priorizará técnicas computacionais, por meio de simulação, para uma análise mais prática e abrangente.

Cenários

Para comparar o desempenho dos estimadores, decidimos seguir o seguinte roteiro:

• Gerar 100 observações da distribuição desejada

• Calcular a média aritmética, a trimédia, e o mid-range

• Armazenar os resultados em um data frame

• Repetir o processo (1000 vezes)

• Analisar a distribuição das médias e comparar a média das médias (\(\hat{\mu}\)) com o verdadeiro valor de \(\mu\)

• Calcular o erro quadrático médio, o erro padrão e a eficiência relativa (a pares) dos estimadores

Considerar as linhas azuis dos gráficos de dispersão como o valor do parâmetro da distribuição e as linhas vermelhas como a média das médias calculadas.

set.seed(13)

n <- 1000
df <- data.frame(media_arit = numeric(n), trimedia = numeric(n), mid_range = numeric(n))

for (i in 1:n) {
  dados <- rnorm(100)
  df[i, "media_arit"] <- mean(dados)
  df[i, "trimedia"] <- (quantile(dados, 0.25) + 2 * quantile(dados, 0.5) + quantile(dados, 0.75)) / 4
  df[i, "mid_range"] <- (quantile(dados, 0.25) + quantile(dados, 0.75)) / 2
}

# Histograma
par(mfrow = c(1,3))
hist(df$media_arit, col = "salmon", main = "Média Aritmética", xlab = "mu.hat", xlim = c(-0.4, 0.4), ylim = c(0, 200))
hist(df$trimedia, col = "salmon", main = "Trimédia", xlab = "mu.hat", xlim = c(-0.4, 0.4), ylim = c(0, 200))
hist(df$mid_range, col = "salmon", main = "Mid-range", xlab = "mu.hat", xlim = c(-0.4, 0.4), ylim = c(0, 200))

# Scatterplot
par(mfrow = c(1,3))

plot(df$media_arit ~ c(1:1000), col = "darkgrey", 
     main = "Média Aritmética", xlab = "Amostra", ylab = "mu.hat", ylim = c(-0.4, 0.3))
abline(h = 0, col = "steelblue", lwd = 1.8)
abline(h = mean(df$media_arit), col = "#ff5454", lwd = 1.8)

plot(df$trimedia ~ c(1:1000), col = "darkgrey", 
     main = "Trimédia", xlab = "Amostra", ylab = "mu.hat", ylim = c(-0.4, 0.3))
abline(h = 0, col = "steelblue", lwd = 1.8)
abline(h = mean(df$trimedia), col = "#ff5454", lwd = 1.8)

plot(df$mid_range ~ c(1:1000), col = "darkgrey", 
     main = "Mid-range", xlab = "Amostra", ylab = "mu.hat", ylim = c(-0.4, 0.3))
abline(h = 0, col = "steelblue", lwd = 1.8)
abline(h = mean(df$mid_range), col = "#ff5454", lwd = 1.8)

dados2 <- c(mean(df$media_arit), mean(df$trimedia), mean(df$mid_range),sd(df$media_arit), sd(df$trimedia), sd(df$mid_range), (var(df$media_arit)+(mean(df$media_arit))^2), (var(df$trimedia)+(mean(df$trimedia))^2), (var(df$mid_range)+(mean(df$mid_range))^2))

tabela <- data.frame(t(matrix(dados2,3,3)))
coluna <- c("Media Aritmética", "Trimédia", "Mid-range")
linha <- c("Viés", "Erro Padrão", "EQM")
tabela <- `colnames<-`(tabela,coluna)
tabela <- `row.names<-`(tabela, linha)

kable(tabela)

# Média Aritmética x Trimédia
tabela[3,1]/tabela[3,2]

## [1] 0.7840854

# Média Aritmética x Mid-range
tabela[3,1]/tabela[3,3]

## [1] 0.8169422

# Trimédia x Mid Range
tabela[3,2]/tabela[3,3]

## [1] 1.041905

set.seed(13)

n <- 1000
df <- data.frame(media_arit = numeric(n), trimedia = numeric(n), mid_range = numeric(n))

for (i in 1:n) {
  dados <- rnorm(99) ; dados[100] <- 25
  df[i, "media_arit"] <- mean(dados)
  df[i, "trimedia"] <- (quantile(dados, 0.25) + 2 * quantile(dados, 0.5) + quantile(dados, 0.75)) / 4
  df[i, "mid_range"] <- (quantile(dados, 0.25) + quantile(dados, 0.75)) / 2
}

# Histograma
par(mfrow = c(1,3))
hist(df$media_arit, col = "salmon", main = "Média Aritmética", xlab = "mu.hat", xlim = c(-0.3, 0.6), ylim = c(0, 200))
hist(df$trimedia, col = "salmon", main = "Trimédia", xlab = "mu.hat", xlim = c(-0.3, 0.6), ylim = c(0, 200))
hist(df$mid_range, col = "salmon", main = "Mid-range", xlab = "mu.hat", xlim = c(-0.3, 0.6), ylim = c(0, 200))

# Scatterplot
par(mfrow = c(1,3))

plot(df$media_arit ~ c(1:1000), col = "darkgrey", 
     main = "Média Aritmética", xlab = "Amostra", ylab = "mu.hat", ylim = c(-0.4, 0.6))
abline(h = 0, col = "steelblue", lwd = 1.8)
abline(h = mean(df$media_arit), col = "#ff5454", lwd = 1.8)

plot(df$trimedia ~ c(1:1000), col = "darkgrey", 
     main = "Trimédia", xlab = "Amostra", ylab = "mu.hat", ylim = c(-0.4, 0.6))
abline(h = 0, col = "steelblue", lwd = 1.8)
abline(h = mean(df$trimedia), col = "#ff5454", lwd = 1.8)

plot(df$mid_range ~ c(1:1000), col = "darkgrey", 
     main = "Mid-range", xlab = "Amostra", ylab = "mu.hat", ylim = c(-0.4, 0.6))
abline(h = 0, col = "steelblue", lwd = 1.8)
abline(h = mean(df$mid_range), col = "#ff5454", lwd = 1.8)

dados2 <- c(mean(df$media_arit), mean(df$trimedia), mean(df$mid_range),sd(df$media_arit), sd(df$trimedia), sd(df$mid_range), (var(df$media_arit)+(mean(df$media_arit))^2), (var(df$trimedia)+(mean(df$trimedia))^2), (var(df$mid_range)+(mean(df$mid_range))^2))

tabela <- data.frame(t(matrix(dados2,3,3)))
coluna <- c("Media Aritmética", "Trimédia", "Mid-range")
linha <- c("Viés", "Erro Padrão", "EQM")
tabela <- `colnames<-`(tabela,coluna)
tabela <- `row.names<-`(tabela, linha)

kable(tabela)

# Média Aritmética x Trimédia
tabela[3,1]/tabela[3,2]

## [1] 5.883322

# Média Aritmética x Mid-range
tabela[3,1]/tabela[3,3]

## [1] 5.864696

# Trimédia x Mid Range
tabela[3,2]/tabela[3,3]

## [1] 0.9968342

set.seed(13)

n <- 1000
df <- data.frame(media_arit = numeric(n), trimedia = numeric(n), mid_range = numeric(n))

for (i in 1:n) {
  dados <- rexp(100)
  df[i, "media_arit"] <- mean(dados)
  df[i, "trimedia"] <- (quantile(dados, 0.25) + 2 * quantile(dados, 0.5) + quantile(dados, 0.75)) / 4
  df[i, "mid_range"] <- (quantile(dados, 0.25) + quantile(dados, 0.75)) / 2
}

# Histograma
par(mfrow = c(1,3))
hist(df$media_arit, col = "salmon", main = "Média Aritmética", xlab = "mu.hat", xlim = c(0.4, 1.4), ylim = c(0, 250))
hist(df$trimedia, col = "salmon", main = "Trimédia", xlab = "mu.hat", xlim = c(0.4, 1.4), ylim = c(0, 250))
hist(df$mid_range, col = "salmon", main = "Mid-range", xlab = "mu.hat", xlim = c(0.4, 1.4), ylim = c(0, 250))

# Scatterplot
par(mfrow = c(1,3))

plot(df$media_arit ~ c(1:1000), col = "darkgrey", 
     main = "Média Aritmética", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.5, 1.4))
abline(h = 1, col = "steelblue", lwd = 1.8)
abline(h = mean(df$media_arit), col = "#ff5454", lwd = 1.8)

plot(df$trimedia ~ c(1:1000), col = "darkgrey", 
     main = "Trimédia", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.5, 1.4))
abline(h = 1, col = "steelblue", lwd = 1.8)
abline(h = mean(df$trimedia), col = "#ff5454", lwd = 1.8)

plot(df$mid_range ~ c(1:1000), col = "darkgrey", 
     main = "Mid-range", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.5, 1.4))
abline(h = 1, col = "steelblue", lwd = 1.8)
abline(h = mean(df$mid_range), col = "#ff5454", lwd = 1.8)

dados2 <- c(mean(df$media_arit)-1, mean(df$trimedia)-1, mean(df$mid_range)-1,sd(df$media_arit), sd(df$trimedia), sd(df$mid_range), (var(df$media_arit)+(mean(df$media_arit)-1)^2), (var(df$trimedia)+(mean(df$trimedia)-1)^2), (var(df$mid_range)+(mean(df$mid_range)-1)^2))

tabela <- data.frame(t(matrix(dados2,3,3)))
coluna <- c("Media Aritmética", "Trimédia", "Mid-range")
linha <- c("Viés", "Erro Padrão", "EQM")
tabela <- `colnames<-`(tabela,coluna)
tabela <- `row.names<-`(tabela, linha)

kable(tabela)

# Média Aritmética x Trimédia
tabela[3,1]/tabela[3,2]

## [1] 0.1536921

# Média Aritmética x Mid-range
tabela[3,1]/tabela[3,3]

## [1] 0.2640429

# Trimédia x Mid Range
tabela[3,2]/tabela[3,3]

## [1] 1.717999

set.seed(13)

n <- 1000
df <- data.frame(media_arit = numeric(n), trimedia = numeric(n), mid_range = numeric(n))

for (i in 1:n) {
  dados <- runif(100)
  df[i, "media_arit"] <- mean(dados)
  df[i, "trimedia"] <- (quantile(dados, 0.25) + 2 * quantile(dados, 0.5) + quantile(dados, 0.75)) / 4
  df[i, "mid_range"] <- (quantile(dados, 0.25) + quantile(dados, 0.75)) / 2
}

# Histograma
par(mfrow = c(1,3))
hist(df$media_arit, col = "salmon", main = "Média Aritmética", xlab = "mu.hat", xlim = c(0.35,0.65), ylim = c(0,300))
hist(df$trimedia, col = "salmon", main = "Trimédia", xlab = "mu.hat", xlim = c(0.35,0.65), ylim = c(0,300))
hist(df$mid_range, col = "salmon", main = "Mid-range", xlab = "mu.hat", xlim = c(0.35,0.65), ylim = c(0,300))

# Scatterplot
par(mfrow = c(1,3))

plot(df$media_arit ~ c(1:1000), col = "darkgrey", 
     main = "Média Aritmética", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.35, 0.6))
abline(h = 0.5, col = "steelblue", lwd = 1.8)
abline(h = mean(df$media_arit), col = "#ff5454", lwd = 1.8)

plot(df$trimedia ~ c(1:1000), col = "darkgrey", 
     main = "Trimédia", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.35, 0.6))
abline(h = 0.5, col = "steelblue", lwd = 1.8)
abline(h = mean(df$trimedia), col = "#ff5454", lwd = 1.8)

plot(df$mid_range ~ c(1:1000), col = "darkgrey", 
     main = "Mid-range", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.35, 0.6))
abline(h = 0.5, col = "steelblue", lwd = 1.8)
abline(h = mean(df$mid_range), col = "#ff5454", lwd = 1.8)

dados2 <- c(mean(df$media_arit)-0.5, mean(df$trimedia)-0.5, mean(df$mid_range)-0.5,sd(df$media_arit), sd(df$trimedia), sd(df$mid_range), (var(df$media_arit)+(mean(df$media_arit)-0.5)^2), (var(df$trimedia)+(mean(df$trimedia)-0.5)^2), (var(df$mid_range)+(mean(df$mid_range)-0.5)^2))

tabela <- data.frame(t(matrix(dados2,3,3)))
coluna <- c("Media Aritmética", "Trimédia", "Mid-range")
linha <- c("Viés", "Erro Padrão", "EQM")
tabela <- `colnames<-`(tabela,coluna)
tabela <- `row.names<-`(tabela, linha)

kable(tabela)

# Média Aritmética x Trimédia
tabela[3,1]/tabela[3,2]

## [1] 0.5333539

# Média Aritmética x Mid-range
tabela[3,1]/tabela[3,3]

## [1] 0.654944

# Trimédia x Mid Range
tabela[3,2]/tabela[3,3]

## [1] 1.227973

set.seed(13)

n <- 1000
df <- data.frame(media_arit = numeric(n), trimedia = numeric(n), mid_range = numeric(n))

for (i in 1:n) {
  dados <- rpois(100, lambda = 5)
  df[i, "media_arit"] <- mean(dados)
  df[i, "trimedia"] <- (quantile(dados, 0.25) + 2 * quantile(dados, 0.5) + quantile(dados, 0.75)) / 4
  df[i, "mid_range"] <- (quantile(dados, 0.25) + quantile(dados, 0.75)) / 2
}

# Histograma
par(mfrow = c(1,3))
hist(df$media_arit, col = "salmon", main = "Média Aritmética", xlab = "mu.hat", xlim = c(3,6), ylim = c(0,400))
hist(df$trimedia, col = "salmon", main = "Trimédia", xlab = "mu.hat", xlim = c(3,6), ylim = c(0,400))
hist(df$mid_range, col = "salmon", main = "Mid-range", xlab = "mu.hat", xlim = c(3,6), ylim = c(0,400))

# Scatterplot
par(mfrow = c(1,3))

plot(df$media_arit ~ c(1:1000), col = "darkgrey", 
     main = "Média Aritmética", xlab = "Amostra", ylab = "mu.hat", ylim = c(4, 6))
abline(h = 5, col = "steelblue", lwd = 1.8)
abline(h = mean(df$media_arit), col = "#ff5454", lwd = 1.8)

plot(df$trimedia ~ c(1:1000), col = "darkgrey", 
     main = "Trimédia", xlab = "Amostra", ylab = "mu.hat", ylim = c(4, 6))
abline(h = 5, col = "steelblue", lwd = 1.8)
abline(h = mean(df$trimedia), col = "#ff5454", lwd = 1.8)

plot(df$mid_range ~ c(1:1000), col = "darkgrey", 
     main = "Mid-range", xlab = "Amostra", ylab = "mu.hat", ylim = c(4, 6))
abline(h = 5, col = "steelblue", lwd = 1.8)
abline(h = mean(df$mid_range), col = "#ff5454", lwd = 1.8)

dados2 <- c(mean(df$media_arit)-5, mean(df$trimedia)-5, mean(df$mid_range)-5,sd(df$media_arit), sd(df$trimedia), sd(df$mid_range), (var(df$media_arit)+(mean(df$media_arit)-5)^2), (var(df$trimedia)+(mean(df$trimedia)-5)^2), (var(df$mid_range)+(mean(df$mid_range)-5)^2))

tabela <- data.frame(t(matrix(dados2,3,3)))
coluna <- c("Media Aritmética", "Trimédia", "Mid-range")
linha <- c("Viés", "Erro Padrão", "EQM")
tabela <- `colnames<-`(tabela,coluna)
tabela <- `row.names<-`(tabela, linha)

kable(tabela)

# Média Aritmética x Trimédia
tabela[3,1]/tabela[3,2]

## [1] 0.5782791

# Média Aritmética x Mid-range
tabela[3,1]/tabela[3,3]

## [1] 0.3519523

# Trimédia x Mid Range
tabela[3,2]/tabela[3,3]

## [1] 0.6086201

set.seed(13)

n <- 1000
df <- data.frame(media_arit = numeric(n), trimedia = numeric(n), mid_range = numeric(n))

for (i in 1:n) {
  dados <- rpois(100, lambda = 1)
  df[i, "media_arit"] <- mean(dados)
  df[i, "trimedia"] <- (quantile(dados, 0.25) + 2 * quantile(dados, 0.5) + quantile(dados, 0.75)) / 4
  df[i, "mid_range"] <- (quantile(dados, 0.25) + quantile(dados, 0.75)) / 2
}

# Histograma
par(mfrow = c(1,3))
hist(df$media_arit, col = "salmon", main = "Média Aritmética", xlab = "mu.hat", xlim = c(0.2,1.6), ylim = c(0,600))
hist(df$trimedia, col = "salmon", main = "Trimédia", xlab = "mu.hat", xlim = c(0.2,1.6), ylim = c(0,600))
hist(df$mid_range, col = "salmon", main = "Mid-range", xlab = "mu.hat", xlim = c(0.2,1.6), ylim = c(0,600))

# Scatterplot
par(mfrow = c(1,3))

plot(df$media_arit ~ c(1:1000), col = "darkgrey", 
     main = "Média Aritmética", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.2, 1.4))
abline(h = 1, col = "steelblue", lwd = 1.8)
abline(h = mean(df$media_arit), col = "#ff5454", lwd = 1.8)

plot(df$trimedia ~ c(1:1000), col = "darkgrey", 
     main = "Trimédia", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.2, 1.4))
abline(h = 1, col = "steelblue", lwd = 1.8)
abline(h = mean(df$trimedia), col = "#ff5454", lwd = 1.8)

plot(df$mid_range ~ c(1:1000), col = "darkgrey", 
     main = "Mid-range", xlab = "Amostra", ylab = "mu.hat", ylim = c(0.2, 1.4))
abline(h = 1, col = "steelblue", lwd = 1.8)
abline(h = mean(df$mid_range), col = "#ff5454", lwd = 1.8)

dados2 <- c(mean(df$media_arit)-1, mean(df$trimedia)-1, mean(df$mid_range)-1,sd(df$media_arit), sd(df$trimedia), sd(df$mid_range), (var(df$media_arit)+(mean(df$media_arit)-1)^2), (var(df$trimedia)+(mean(df$trimedia)-1)^2), (var(df$mid_range)+(mean(df$mid_range)-1)^2))

tabela <- data.frame(t(matrix(dados2,3,3)))
coluna <- c("Media Aritmética", "Trimédia", "Mid-range")
linha <- c("Viés", "Erro Padrão", "EQM")
tabela <- `colnames<-`(tabela,coluna)
tabela <- `row.names<-`(tabela, linha)

kable(tabela)

# Média Aritmética x Trimédia
tabela[3,1]/tabela[3,2]

## [1] 0.4082183

# Média Aritmética x Mid-range
tabela[3,1]/tabela[3,3]

## [1] 0.1132121

# Trimédia x Mid Range
tabela[3,2]/tabela[3,3]

## [1] 0.2773322

Conclusão

A média aritmética apresentou um bom desempenho na maioria dos cenários, apesar de ser muito sensível a outliers. Os estimadores trimédia e mid-range tiveram performances similares e se destacaram na presença de outliers, sendo ótimos substitutos ao estimador convencional nesse sentido. Por fim, é interessante destacar que a análise gráfica foi importante, ao percebermos alguns padrões que passariam em branco nos cálculos.