Esse artigo é uma explicação do código descrito nesse link referente a uma palestra de valor esperado do curso de Inferência Estatística da especialização em Ciências de Dados do Coursera.
Primeiro, definiremos as bibliotecas que vamos usar, o número de simulações nosim e os tamanhos das amostras para calcular as médias em amostra, e a semente da simulação para que seja reprodutível.
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
nosim <- 10000
amostra <- c(10,20,30)
set.seed(123)
Os dataframes serão criados com duas colunas, uma com a média dos lançamentos e outra com o tamanho da amostra. A coluna com a média será determinada pelos os resultados dos lançamentos simples e dos lançamentos baseadas nos tamanhos dados em amostra. Por exemplo, se o primeiro elemento de amostra é 10, cada elemento relativa a 10 terá a média de 10 lançamentos simulados. Para fazer esses cálculos mais rapidamente, criaremos uma matriz com um número de linhas igual a nosim e o número de colunas igual ao valor correspondente na amostra, com cada entrada referente a um lançamento com reposição, determinado pela função sample. Após isso, utilizaremos a função apply em cada linha de cada matriz para calcular a média. Faremos um exemplo abaixo.
Vamos supor que 0 represente “cara” e 1 “coroa” num lançamento de moedas. Vamos criar uma matriz com 10 linhas, onde cada linha representa 10 lançamentos aleatórios de moeda, e armazenar na variável exemplo.
exemplo <- matrix(sample(0 : 1, 10 * amostra[1], replace = TRUE),
10)
exemplo
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 1 0 0 0 1 0 1 1 1
## [2,] 0 1 1 1 1 0 1 0 1 1
## [3,] 0 1 0 0 1 0 1 0 0 0
## [4,] 1 0 0 1 0 0 0 0 1 1
## [5,] 0 1 0 1 0 0 0 0 1 1
## [6,] 1 0 0 0 0 1 1 1 1 0
## [7,] 1 1 1 0 0 1 0 1 0 0
## [8,] 1 0 1 0 1 0 0 0 0 1
## [9,] 0 0 0 0 0 1 0 1 1 0
## [10,] 0 0 1 1 0 0 0 1 0 0
Agora, vamos aplicar a função apply nessa matriz, aplicando a função mean na linhas (número 1 na fórmula).
apply(exemplo,1,mean)
## [1] 0.5 0.7 0.3 0.4 0.4 0.5 0.5 0.4 0.3 0.3
Utilizaremos esse procedimento para cada valor em amostra, com cada matriz tendo o tamanho de nosim. Esses valores serão juntados na variável x do dataframe. Para terminar o dataframe, criaremos uma coluna do tipo factor para dizer qual o tamanho da amostra de cada lançamento. Assim, para nosso exemplo, o dataframe terá 400000 linhas e duas colunas.
Com isso, o dataframe para as simulações de lançamentos de moedas vai ser dado mostrado abaixo. Adicionaremos uma amostra de linhas desse dataframe e vamos criar o dataframe mu que será a média total das amostras.
dat <- data.frame(
#Lançamentos simples
x = c(sample(0 : 1, nosim, replace = TRUE),
#Médias de lançamentos para primeiro elemento de amostra
apply(matrix(sample(0 : 1, nosim * amostra[1], replace = TRUE),
nosim), 1, mean),
#Médias de lançamentos para segundo elemento de amostra
apply(matrix(sample(0 : 1, nosim * amostra[2], replace = TRUE),
nosim), 1, mean),
#Médias de lançamentos para terceiro elemento de amostra
apply(matrix(sample(0 : 1, nosim * amostra[3], replace = TRUE),
nosim), 1, mean)
),
#fatores
size = factor(rep(c(1, amostra[1], amostra[2], amostra[3]),
rep(nosim, 4)))
)
dat[sample(1:40000,10),]
## x size
## 16229 0.6000000 10
## 16742 0.5000000 10
## 33643 0.4666667 30
## 1488 1.0000000 1
## 5228 1.0000000 1
## 14061 0.6000000 10
## 8860 1.0000000 1
## 14447 0.6000000 10
## 35467 0.4333333 30
## 21383 0.3500000 20
#Dataframe com a média por grupo
mu <- dat %>% group_by(size) %>% summarise(grp.mean=mean(x),.groups="drop")
mu
## # A tibble: 4 x 2
## size grp.mean
## <fct> <dbl>
## 1 1 0.499
## 2 10 0.498
## 3 20 0.500
## 4 30 0.499
O valor da média grp.mean corresponde ao valor esperado calculado:
\[E[X] = 0.5 \times 0 + 0.5 \times 1 = 0.5 \]
Agora faremos um gráfico com a biblioteca ggplot2, contendo um histograma das médias para cada fator, e uma linha tracejada em cada histograma com a média total de cada fator.
g <- ggplot(dat, aes(x = x, fill = size)) +
geom_histogram( binwidth = 1/12, colour = "black")+
geom_density(alpha=0.2)+ facet_grid(. ~ size)+
geom_vline(data=mu, aes(xintercept=grp.mean),
linetype="dashed")+
labs( x="Média das amostras",y = "Número de lançamentos",fill="Tamanho \n amostra")
g
Distribuição de médias de amostras de lançamentos de moedas
Faremos o mesmo do caso anterior, só trocando a amostra, que antes era de entre 0 e 1 no caso das moedas e agora faremos entre 1 e 6 representando cada face do dado.
dat <- data.frame(
#Lançamentos simples
x = c(sample(1 : 6, nosim, replace = TRUE),
#Médias de lançamentos para primeiro elemento de amostra
apply(matrix(sample(1 : 6, nosim * amostra[1], replace = TRUE),
nosim), 1, mean),
#Médias de lançamentos para segundo elemento de amostra
apply(matrix(sample(1 : 6, nosim * amostra[2], replace = TRUE),
nosim), 1, mean),
#Médias de lançamentos para terceiro elemento de amostra
apply(matrix(sample(1 : 6, nosim * amostra[3], replace = TRUE),
nosim), 1, mean)
),
#fatores
size = factor(rep(c(1, amostra[1], amostra[2], amostra[3]),
rep(nosim, 4)))
)
dat[sample(1:40000,10),]
## x size
## 36487 3.433333 30
## 12403 3.400000 10
## 17617 3.700000 10
## 17386 4.600000 10
## 355 3.000000 1
## 17040 1.700000 10
## 13165 3.600000 10
## 10373 3.200000 10
## 498 5.000000 1
## 1438 2.000000 1
#Dataframe com a média por grupo
mu <- dat %>% group_by(size) %>% summarise(grp.mean=mean(x),.groups="drop")
mu
## # A tibble: 4 x 2
## size grp.mean
## <fct> <dbl>
## 1 1 3.48
## 2 10 3.49
## 3 20 3.50
## 4 30 3.50
O valor da média grp.mean corresponde ao valor esperado calculado:
\[ E[X] = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} + 3 \times \frac{1}{6} + 4 \times \frac{1}{6} + 5 \times \frac{1}{6} + 6 \times \frac{1}{6} = 3.5 \]
E o gráfico é dado por
g <- ggplot(dat, aes(x = x, fill = size)) +
geom_histogram(alpha=0.4, binwidth = .25, colour = "black")+
geom_density(alpha=0.2)+ facet_grid(. ~ size)+
geom_vline(data=mu, aes(xintercept=grp.mean),
linetype="dashed")+
labs(x="Média das amostras",y = "Número de lançamentos",fill="Tamanho \n amostra")
g
Distribuição de médias de amostras de lançamentos de dados