Última atualização: 13 de maio, 2025 às 20:08.

1 Distribuição amostral da média

A figura a seguir copiada do livro Estatística Básica dos Profs Wilton Bussab e Pedro Morettin (Bussab and Morettin 2017) é uma ilustração do Teorema Central do Limite. Faça um código computacional para produzir um resultado semelhante.

Figura 10.5: Histogramas correspondentes à distribuição amostral da média amostral para amostras extraídas de algumas populações.

Seguem algumas ideias e códigos iniciais para o Caso 1 da figura.

## estimativa para 1 amostra simulada
(am <- runif(2))
## [1] 0.1813371 0.3105324
(m1 <- mean(am))
## [1] 0.2459347
## Repetindo para várias (N) amostras de tamanho n
N <- 100000
n <- 10

ams <- matrix(runif(n*N), ncol=n)
estimativas <- rowMeans(ams)
length(estimativas)
## [1] 100000
hist(estimativas, freq=FALSE)
lines(density(estimativas))

Para facilitar produzir resultados e gráficos para diferentes tamanhos de amostra \(n\) fazemos um função que combina os comandos para simulação e gráfico da distribuição amostral (representada pelo historiograma e gráfico de densidade) e a aproximação normal (representada pela linha vermelha).

## Montando uma função
damUnif <- function(n, N=1e5, plot=TRUE){
    ams <- matrix(runif(n*N), ncol=n)
    estimativas <- rowMeans(ams)
    if(plot){
        hist(estimativas, freq=FALSE, col = "transparent",
             main = bquote(n == .(n)), nclass = 25, xlab = expression(bar(y)))
        lines(density(estimativas))
        EX <- (0+1)/2
        VX <- ((1-0)^2)/12
        curve(dnorm(x, m=EX, sd=sqrt(VX/n)),from=0, to=1, col=2, add=TRUE)
    }
    return(invisible())
}
par(mfrow = c(1,4), mar = c(3.5, 3.5, 1, 1), mgp = c(2,1,0))
POP1 <- runif(1e5)
hist(POP1, freq = FALSE, main = "população", xlab = "y")
curve(dunif(x), add = TRUE)
damUnif(2)
damUnif(5)
damUnif(25)

Uma outra opção é sortear uma população (finita) e dela extrair amostras aleatórias.

## Pop finita
POP <- runif(1536)
(am <- sample(POP, 2))
mean(am)

Vejamos agora o Caso 4 da figura que representa uma mistura de (duas) distribuições.

Sejam \(f_1(y), f_2(y), \ldots, f_k(y)\) e pesos \(w_1, w_2, \ldots, w_k\) com \(w_i > 0\) e \(\sum_i w_i = 1\). A densidade de probabilidades \[ \sum_i w_i f_i(y) \] é uma mistura finita com \[\begin{align*} \mu = \text{E}[Y] &= \sum_i w_i \mu_i \\ \sigma^2 = \text{Var}[Y] &= \sum_i w_i (\sigma^2_i + \mu_i^2) - \mu^2 \end{align*}\]

Para este caso vamos escolher \(k = 2\), mistura de duas densidades normais sendo a primeira \(\text{N}(10, 1)\) com peso \(w_1 = w = 0.3\) e a segunda \(\text{N}(15, 1)\) com \(w_2 = (1 - w) = 0.7\). A densidade e suas duas componentes são mostradas na figura a seguir.

m1 <- 10; s1 <- 1
m2 <- 15; s2 <- 1
w1 <- 0.3
LIMS <- c(m1 - 3*s1, m2 + 3*s2)
dmixnorm <- function(x, mu1 = m1, sd1 = s1, mu2 = m2, sd2 = s2, w=w1){
    w*dnorm(x, mu1, sd1) + (1-w)*dnorm(x, mu2, sd2)
}
curve(dmixnorm(x), from = LIMS[1], to = LIMS[2])
curve(w1*dnorm(x, m1, s1), from = LIMS[1], to = LIMS[2], col = "blue", add = TRUE)
curve((1-w1)*dnorm(x, m2, s2), from = LIMS[1], to = LIMS[2], col = "red", add = TRUE)

Amostras da distribuição de mistura podem ser obtidas de diversas formas e usamos uma a seguir.

POP4a <- rnorm(1e5, m1, s1)    
POP4b <- rnorm(1e5, m2, s2)
IND <- sample(c(0,1), 1e5, prob = c(w1, 1-w1), replace = TRUE)
POP4 <- ifelse(IND == 0 , POP4a, POP4b)
rm(POP4a, POP4b, IND)
hist(POP4, freq = FALSE, nclass = 25)
curve(dmixnorm(x), from = LIMS[1], to = LIMS[2], add = TRUE)

A seguir são obtidas as distribuições amostrais para diferentes tamanho de amostra e a aproximação normal é indicada na linha vermelha.

damMix <- function(n, N=1e5, plot=TRUE,
                   mu1 = m1, sd1 = s1, mu2 = m2, sd2 = s2, w=w1){
    Nn <- N*n
    Nn1 <- round(w * Nn)
    Nn2 <- Nn - Nn1
    ##
    POPa <- rnorm(Nn1, mu1, sd1)
    POPb <- rnorm(Nn2, mu2, sd1)
    ams <- matrix(sample(c(POPa, POPb)), ncol=n)
    estimativas <- rowMeans(ams)
    if(plot){
        h <- hist(estimativas, freq=FALSE, col = "transparent",
                  main = bquote(n == .(n)), nclass = 25, xlab = expression(bar(y)))
        lines(density(estimativas))
        EX <- w * mu1 + (1-w) * mu2
        VX <- sum(c(w,1-w)*(c(sd1,sd2)^2 + c(mu1, mu2)^2)) - EX^2
        curve(dnorm(x, m=EX, sd=sqrt(VX/n)),from=h$beaks[1], to=h$beaks[1], col=2, add=TRUE)
    }
    return(invisible())
}

par(mfrow = c(1,4), mar = c(3.5, 3.5, 1, 1), mgp = c(2,1,0))
hist(POP4, freq = FALSE, nclass = 25, main = "população", xlab = "y")
curve(dmixnorm(x), from = LIMS[1], to = LIMS[2], add = TRUE, lwd = 2)
damMix(3)
damMix(5)
damMix(10)

Existem pacotes para cálculos estatísticos envolvendo misturas que contém funções diversas que podem ser utilizadas.

2 Distribuição amostral da variância

2.1 Verificando propriedades do estimador

Implemente uma ilustração computacional para avaliar propriedades do estimador da variância de uma população normal.
\[ S^2 = \frac{\sum_{i=1}^n (y_i - \overline{y})^2}{n-1} \]

Lembre-se que, neste caso, temos uma distribuição amostral conhecida para \(S^2\) que pode ser usada como referência. \[ \frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1}. \]

Avalie as propriedades de (não) tendenciosidade (vício), variância/erro padrão, erro quadrático médio, eficiência (relativa) e consistência. Veja no final desta página uma rápida revisão das propriedades dos estimadores.
Considere diferentes tamanhos de amostra, por exemplo, \(n=10, 20, 30\).

Algumas ideias para começar:

N <- 100000
n <- 10
mu <- 70; sigma <- 10

(am <- rnorm(n, m=mu, sd=sigma))
var(am)                     ## S^2
(n-1)*var(am)/n             ## \hat{\sigma}^2
(mean(abs(am-mean(am))))^2  ## DM^2

ams <- matrix(rnorm(n*N, m=mu, sd=sigma), ncol=n)
dim(ams)
S2 <- apply(ams, 1, var)
S2sc <- (n-1)*S2/sigma^2   ## escalonando
## ou poderia ter calculado diretamente por:
S2sc <- apply(ams, 1,function(y) (n-1)*var(y)/(sigma^2))
par(mfrow=c(1,2), mgp = c(2,1,0), mar=c(3.5,3.5,2,0))
hist(S2, freq=FALSE)
hist(S2sc, freq=FALSE)
lines(density(S2sc))
curve(dchisq(x, df=n-1), from=0, to=30, col=2, add=T)

Vício: \(E[S^2] = \sigma^2\) ?

mean(S2)
## [1] 99.93377
sigma^2
## [1] 100
mean(S2) - sigma^2  ## viés
## [1] -0.06623342
(mean(S2) - sigma^2)/sigma^2  ## viés relativo
## [1] -0.0006623342

## variância e erro padrão do estimador
var(S2)
## [1] 2211.732
sd(S2)
## [1] 47.02905

E agora faça para uma sequencia de valores de n e verifique o efeito do tamanho da amostra.

2.2 Comparando estimadores

Considere outros estimadores a seguir (e/ou outros que queira propor/adotar). \[ \hat{\sigma}^2 = \frac{\sum_{i=1}^n (y_i - \overline{y})^2}{n} \;\;;\;\; {\rm DM}^2 = \left(\frac{\sum_{i=1}^n |y_i - \overline{y}|}{n}\right)^2 \]

2.3 Distribuição amostral da variância: população não normal.

Repita o exercício anterior para observações provenientes de \(Y \sim {\rm Gama}(2,4)\).
Compare e discuta os resultados.

Apêndice: Propriedades de estimadores

Seja \(T\) um estimador de um parâmetro \(\theta\).

  1. O estimador é dito não-viciado se \({\rm E}[T] = \theta\).
  2. A variância do estimador é \({\rm E}[(T - \text{E}[T])^2]\) e o erro padrão é a raiz quadrada deste valor.
  3. O erro quadrático médio (MSE de sigla em inglês) é \({\rm MSE}(T) = {\rm E}[(T - \theta)^2] = {\rm Var[T]} + (\text{E}(T) - \theta)^2\).
  4. Um estimador \(T_2\) é dito mais eficiente que \(T_1\) se \(\text{Var}(T_1) > \text{Var}(T_2)\) para todos valores de \(\theta\).
  5. \(T\) é dito consistente se \(\underset{n \to \infty}{\lim} \text{Pr}(|T_n - \theta| > \epsilon) = 0\).

References

Bussab, Wilton, and Pedro Morettin. 2017. Estatística Básica. 9th ed. Saraiva.