Comandos R para análises estatísticas

Capítulo 10: Introdução à Inferência Estatística

10.5 Amostra Aleatória Simples

Para ilustrar como selecionar uma amostra aleatória no R, utilizemos os dados da Tabela 2.1 do livro. Temos, no total 36 indivíduos naquela tabela e vamos selecionar 15 deles ao acaso, de duas formas: com e sem reposição.

##### Amostragem sem reposição: Utilize a opção   'replace = FALSE'
AASs = sample(x = 1:dim(tab2_1)[1], size=15, replace = FALSE )
# Elementos selecionados na amostra
AASs
##  [1] 11 28 14 30 31  2 16 26 33 13 25 12 17 34  3
##### Amostragem com reposição: Utilize a opção   'replace = TRUE'
AASc = sample(x = 1:dim(tab2_1)[1], size=15, replace = TRUE )
# Elementos selecionados na amostra
AASc
##  [1] 33  9  2 12 35 33 25 24 36 24 26 20 22 11  6

O resultado destas duas amostragem é apresentado nas tabelas 10.A e 10.B a seguir:

Tabela 10.A: Amostra Aleatória Simples sem Reposição da Tabela 2.1
N estado_civil grau_instrucao n_filhos salario idade_anos idade_meses reg_procedencia
11 11 casado ensino médio 2 8.12 33 6 interior
28 28 casado ensino médio 0 14.69 29 8 interior
14 14 casado ensino fundamental 3 8.95 44 2 outra
30 30 casado ensino médio 2 15.99 35 10 capital
31 31 solteiro superior NA 16.22 31 5 outra
2 2 casado ensino fundamental 1 4.56 32 10 capital
16 16 solteiro ensino médio NA 9.35 38 8 outra
26 26 casado ensino médio 2 13.60 35 0 outra
33 33 casado superior 3 17.26 43 7 capital
13 13 solteiro ensino médio NA 8.74 37 5 outra
25 25 casado ensino médio 2 13.23 32 5 interior
12 12 solteiro ensino fundamental NA 8.46 27 11 capital
17 17 casado ensino médio 1 9.77 31 7 capital
34 34 solteiro superior NA 18.75 33 7 capital
3 3 casado ensino fundamental 2 5.25 36 5 capital
Tabela 10.B: Amostra Aleatória Simples com Reposição da Tabela 2.1
N estado_civil grau_instrucao n_filhos salario idade_anos idade_meses reg_procedencia
33 33 casado superior 3 17.26 43 7 capital
9 9 casado ensino médio 1 7.59 34 10 capital
2 2 casado ensino fundamental 1 4.56 32 10 capital
12 12 solteiro ensino fundamental NA 8.46 27 11 capital
35 35 casado ensino médio 2 19.40 48 11 capital
33.1 33 casado superior 3 17.26 43 7 capital
25 25 casado ensino médio 2 13.23 32 5 interior
24 24 casado superior 0 12.79 26 1 outra
36 36 casado superior 3 23.30 42 2 interior
24.1 24 casado superior 0 12.79 26 1 outra
26 26 casado ensino médio 2 13.60 35 0 outra
20 20 solteiro ensino médio NA 10.76 37 4 interior
22 22 solteiro ensino médio NA 11.59 34 2 capital
11 11 casado ensino médio 2 8.12 33 6 interior
6 6 casado ensino fundamental 0 6.66 28 0 interior

Exemplo 10.8

Já vimos no capítulo anterior como simular uma v.a. Normal:

x<-rnorm(5,mean=167, sd= 5)
x
## [1] 175.93 169.49 157.17 170.51 164.64

10.7 Distribuições Amostrais

Para ilustrar o comportamento das estatísticas média e variância, vamos simular 10.000 observações de uma distribuição normal com média 167 e variância 25, como no Exemplo 10.8.

Já vimos no capítulo anterior como simular uma v.a. Normal:

x_normal<-rnorm(10000,mean=167, sd= 5)
hist(x_normal, col="lightblue4", border="white",freq = FALSE)

Agora vamos extrair 200 amostras(com reposição) com k=15,30 e 100 observações deste conjunto e, para cada amostra, vamos calcular a média e variância. Para isto, faremos uma iteração com 200 ciclos(função for) e para cada ciclo selecionaremos uma amostra e calcularemos a mediana. Note que a função for começa com { e o ciclo é definido até o fechamento das chaves }. A variável de iteração, neste caso, será i e irá de 1 a 200:

#Inicializando as variaveis como vetores numericos
## Media e variancia para amostras de tam 15
xbar15<-numeric()
var_amostral15<-numeric()
## Media e variancia para amostras de tam 30
xbar30<-numeric()
var_amostral30<-numeric()
## Media e variancia para amostras de tam 100
xbar100<-numeric()
var_amostral100<-numeric()

for ( i in 1:200){
  # Extraindo amostras de tamanho 15 e calculando a média e variancia
  smp<-sample(x_normal,size = 15)
  xbar15[i]<-mean(x_normal[smp])
  var_amostral15[i]<-var(x_normal[smp])
  
  # Extraindo amostras de tamanho 30 e calculando a média e variancia
  smp<-sample(x_normal,size = 30)
  xbar30[i]<-mean(x_normal[smp])
  var_amostral30[i]<-var(x_normal[smp])
  
  # Extraindo amostras de tamanho 100 e calculando a média e variancia
  smp<-sample(x_normal,size = 100)
  xbar100[i]<-mean(x_normal[smp])
  var_amostral100[i]<-var(x_normal[smp])
}

A Figura a seguir apresenta os histogramas das 6 quantidades calculadas acima. Note que, apesar de \(\hat{\sigma}^2\sim \chi^2\), a medida que o número de observações da amostra aumenta, a distribuição de qui-quadrado se aproxima à distribuição normal.

par(mfrow=c(2,3))
hist(xbar15, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Média com n=15")
hist(xbar30, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Média com n=30")
hist(xbar100, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Média com n=100")
hist(var_amostral15, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Variância com n=15")
hist(var_amostral30, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Variância com n=30")
hist(var_amostral100, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Variância com n=100")

Figura 10.A: Histogramas das estatísticas das amostras

Como temos que \(\overline{X}\sim\mathcal{N}(\mu , \sigma^2/n)\), esperamos que, a medida que n aumenta, a distribuição de \(\overline{X}\) fique mais concentrada em torno da média. Veja a Figura 10.B a seguir, com a comparação entre as 4 curvas: a distribuição da Normal simulada e as três distribuições amostrais da média:

plot (density(x_normal), ylim=c(0,1), xlab="", main = "Funções densidades empíricas")
lines (density(xbar15),col="red")
lines (density(xbar30),col="blue")
lines (density(xbar100),col="orange")
legend("topright", c("X~N(167,25)", "Média com n=15", "Média com n=30", "Média com n=100"),
       col=c("black", "red", "blue", "orange"), lty=c(1,1,1,1))

Figura 10.B: Densidades empíricas das v.a. simulada e das médias para n=15, 30 e 100.

Exemplo 10.8

Para selecionar as 200 amostras de tamanho 5, faremos novamente uma iteração com 200 ciclos e para cada ciclo selecionaremos uma amostra e calcularemos a mediana:

## Mediana para amostras de tam 5
md5<-numeric()
for ( i in 1:200){
  # Extraindo amostras de tamanho 5 e calculando a mediana
  smp<-sample(x_normal,size = 5)
  md5[i]<-median(x_normal[smp])
}

Para imprimir as estatísticas descritivas das medianas calculadas, podemos utilizar os seguintes comandos:

print(summary(md5))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     162     165     166     167     167     173
cat("E(md)=",mean(md5), ", Var(md)=",var(md5),", dp(md)=",sd(md5),", min(md)=",min(md5),", max(md)=",max(md5))
## E(md)= 166.5 , Var(md)= 3.7276 , dp(md)= 1.9307 , min(md)= 161.68 , max(md)= 173.32

Finalmente, o gráfico da distribuição empírica (Figura 10.3) é produzido da seguinte forma:

hist(md5, col="lightblue3", border="white",freq = FALSE, breaks=seq(161,175,2), main="Mediana com n=5")

Figura 10.3: Distribuição amostral da mediana, obtida de 200 amostras de tamanho 5 de \(X\sim\mathcal{N}(167,25)\).

Para simular os dados da Figura 10.5, precisamos simular os dados para cada um dos gráficos. O código a seguir produz as simulações necessárias:

n=cbind(c(2, 5,25), # numero de amostras para os gráficos da primeira linha
        c(2,10,30), # numero de amostras para os gráficos da segunda linha
        c(3, 5,10), # numero de amostras para os gráficos da terceira linha
        c(5,15,25)) # numero de amostras para os gráficos da quarta linha
X<-list(
  x1=matrix(NA,nrow=200,ncol=3), # Matriz que armazenara as médias da primeira linha
  x2=matrix(NA,nrow=200,ncol=3), # Matriz que armazenara as médias da segunda linha
  x3=matrix(NA,nrow=200,ncol=3), # Matriz que armazenara as médias da terceira linha
  x4=matrix(NA,nrow=200,ncol=3)) # Matriz que armazenara as médias da quarta linha
for ( i in 1:3){  # iterando ao longo do número de amostras de cada linha
  for ( j in 1: 200){  # Construindo 200 observacoes para cada gráfico  
    X$x1[j,i] =  mean( runif( n=n[i,1])) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao uniforme
    X$x2[j,i] =  mean(  rexp( n=n[i,2], rate=2)) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao exponencial
    X$x3[j,i] =  mean(rchisq( n=n[i,3], df=1)) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao qui-quadrado
    X$x4[j,i] =  mean( rbeta( n=n[i,4], shape1 = .2, shape2 = .2)) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao beta
  }
}

Figura 10.5: Histogramas correspondentes às distribuições amostrais de \(\overline{X}\) para amostras extraídas de algumas populações.

Exemplo 10.11:

Para ilustrarmos o cálculo da probabilidade apresentada no Exemplo 10.5, simularemos 200 amostras de tamanho 100 da distribuição de \(X\sim \mathcal{N}(500,10^2)\) e para cada uma delas calcularemos a média para construir a distribuição empírica de \(\overline{X}\).

x_barr<-sapply(1:200, function(i){return(mean(rnorm(n=100,mean=500,sd=10)))})
hist(x_barr, col = "lightblue3", border = "white", freq=FALSE)

Finalmente, para estas amostras geradas, vejamos quantas delas encontram-se no intervalo de 498 a 502, que é uma estimativa (baseada nas amostras) da probabilidade:

sum( x_barr >= 498 & x_barr <= 502 )/200
## [1] 0.96

Exemplo 10.12

Novamente, simulemos 200 observações da distribuição desejada, isto é, \(Be(p=0.3)\) e calcularemos \(\hat{p}\)

x_barr <- numeric()
x_var <- numeric()
for (i in 1:200){
smp <- rbinom(n = 10, size = 1, prob = 0.3)
x_barr[i] <- mean(smp)
x_var[i] <- var(smp)
}

p_hat = sum( x_barr ) / 200
p_hat
## [1] 0.2905
hist(x_barr, col="lightblue3", border="white", breaks=seq(0,1,0.02))

Figura 10.C: Histograma das amostras geradas para $ p Be(0.3)$.

Note, na Figura 10.C, que a simetria da distribuição de \(\hat{p}\) é prejudicada pelo fato de \(p \in [0,1]\). Se fizéssemos o mesmo exercício com um n maior, teremos a variância de \(\hat{p}\) menor. Por exemplo, utilizando n=100, temos:

x_barr <- numeric()
x_var <- numeric()
for (i in 1:200){
smp <- rbinom(n = 100, size = 1, prob = 0.3)
x_barr[i] <- mean(smp)
x_var[i] <- var(smp)
}

p_hat = sum( x_barr ) / 200
p_hat
## [1] 0.3035
hist(x_barr, col="lightblue3", border="white", breaks=seq(0,1,0.02))

Exemplo 10.13

Para {1,3,5,5,7}, sabemos do Exemplo 10.7 que :

x 1 3 5 7
P(X=x) 1/5 1/5 2/5 1/5

Se simularmos a distribuição empírica de X a partir de amostras de tamanho n=2, teremos:

x<-c(1,3,5,5,7)
smp<-matrix(NA,nrow=200,ncol=2)
for ( i in 1:200){
  smp[i,] <- sample(x, size=2, replace=TRUE)
}

hist(smp, col="lightblue3", border="white")

Também sabemos (Exemplo 10.10) que \(\mu=4.2\), \(\sigma^2=4.16\) e \(E(\overline{X})=4.2\). Para \(Var(\overline{X})\), temos:

População \(\sigma^2=4.16\)
n = 1 \(Var(\overline{X})=4.16\)
n = 2 \(Var(\overline{X})=2.08\)
n = 3 \(Var(\overline{X})=1.39\)

Simularemos então distribuições amostrais para n=3 como anteriormente, para estimarmos os valores esperado e variância de \(S^2\), md e \(\hat(\sigma^2)\):

x<-c(1,3,5,5,7)
S2<-numeric()
md<-numeric()
var_hat<-numeric()

for ( i in 1:200){
  smp<-sample(x, size=3, replace=TRUE)
  S2[i] <- var(smp)
  var_hat[i] <- var(smp)*2/3
  md[i] <- median(smp)
}

Calculando então a média e variância destes estimadores ao longo das 200 amostras, teremos:

Estatística Média Variância
\(S^2\) 4.34 12.61
Mediana 4.58 2.4
\(\hat{\sigma}^2\) 2.89 5.61

Que são valores próximos dos valores apresentados nas Tabelas 10.7 a 10.9 do livro.

As Figuras 10.6 a 10.8 a seguir representam estimativas(distribuições empíricas) da distribuição dos estimadores \(S^2\), md e \(\hat{\sigma}^2\) por meio de simulações:

hist(S2, col="lightblue3", border="white")

Figura 10.6: Distribuição amostral de \(S^2\) para amostras de tamanho n=3 extraídas de {1,3,5,5,7}

hist(md, col="lightblue3", border="white")

Figura 10.7: Distribuição amostral de md para amostras de tamanho n=3 extraídas de {1,3,5,5,7}

hist(var_hat, col="lightblue3", border="white")

Figura 10.8: Distribuição amostral de \(\hat{\sigma}^2\) para amostras de tamanho n=3 extraídas de {1,3,5,5,7}

Exemplo 10.16

Como já vimos antes, R tem os mesmos comandos do S-plus para a construção de amostras:

x<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
sample(x, 7)
## [1] 14  4 15  9  6  5  8
sample(x, 7, replace=T)
## [1]  1  5  7 15 11 13  6

Capítulo Anterior | Introdução | Próximo Capítulo