Teorema Central Limite em R

Amostras aleatórias

Teorema Central do Limite

Teorema: Seja uma amostra aleatória \((X_1,X_2,\ldots,X_n)\) retiradas de uma população com média \(\mu\) e variância \(\sigma\). A distribuição amostral de \(\bar{X}\) aproxima-se, para n grande, de uma distribuição normal com média \(E[\bar{X}]=\mu\) e variância \(\sigma^2/n\).

# Definindo os parâmetros
mu <- 0
sigma <- 1
n <- 3
m <- 4

# Gerando as amostras
for (i in 1:n) {
  X <- rnorm(m, mean = mu, sd = sigma)
  cat("Amostra", i, ":", X, "\n")
}

Amostra 1 : -1.08823 -0.7172377 -1.129133 0.06722497 
Amostra 2 : 0.06446549 0.4328004 0.365717 -0.7146501 
Amostra 3 : 0.273727 -0.1629615 -0.3819411 -0.5278877 

Teorema: Seja \(X\) uma variável aleatória com esperança \(E[X] = \mu\) e variância \(V(X) = \sigma^2\). Seja \(\bar{X}\) a média amostral de uma amostra de tamanho \(n\). Então:

\[ E[\bar{X}]= \mu, \quad V(\bar{X})=\frac{\sigma^2}{n}. \]

Vamos verificar esse teorema através de simulações:

library(ggplot2)
library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

# Definindo os parâmetros
vn <- c()
vmean <- c()
mu <- 2
std <- 1

# Loop para gerar amostras de uma distribuição exponencial
for (n in seq(1, 10000, by = 10)) {
  X <- rexp(n, rate = 1/mu)  # Gerando dados exponenciais com média mu
  vmean <- c(vmean, mean(X))  # Calculando a média amostral
  vn <- c(vn, n)
}

# Plotando os resultados
plot(vn, vmean, type = "l", col = "blue", lwd = 2, xlab = "n", ylab = expression(bar(X)), main = "Média amostral vs. Tamanho da amostra")
abline(h = mu, col = "red", lwd = 2, lty = 2)  # Linha da média populacional
legend("topright", legend = c("Média amostral", "Média populacional"), col = c("blue", "red"), lty = c(1, 2), lwd = 2)

# Criando o gráfico com ggplot2
p <- ggplot(data.frame(vn, vmean), aes(x = vn, y = vmean)) +
  geom_line(color = "blue", size = 1.2) +
  geom_hline(yintercept = mu, color = "red", linetype = "dashed", size = 1.2) +
  labs(title = "Média amostral vs. Tamanho da amostra",
       x = "n", y = "Média Amostral") +
  theme_minimal()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

# Convertendo o gráfico para plotly
ggplotly(p)

Notamos que a média amostral converge para a média populacional \(\mu\) quando aumentamos o tamanho da amostra.

Teorema Central do Limite

Teorema: Seja uma amostra aleatória \((X_1,X_2,\ldots,X_n)\) retiradas de uma população com média \(\mu\) e variância \(\sigma\). A distribuição amostral de \(\bar{X}\) aproxima-se, para n grande, de uma distribuição normal com média \(E[\bar{X}]=\mu\) e variância \(\sigma^2/n\).

Vamos verificar esse teorema através de simulações:

library(ggplot2)
library(fitdistrplus)

Loading required package: MASS

Attaching package: 'MASS'

The following object is masked from 'package:plotly':

    select

Loading required package: survival

library(gganimate)

vS <- c(1, 2, 4, 8, 50, 100, 1000)  # Tamanhos das amostras
S <- 500  # Número de amostras
mu <- 2

# Preparando o dataframe para renderização
df_combined <- data.frame()

for (n in vS) {
  vmean <- numeric(S)
  
  # Gerando amostras e calculando médias
  for (s in 1:S) {
    X <- runif(n, min = 0, max = 2 * mu)  # Amostras de distribuição uniforme
    vmean[s] <- mean(X)
  }
  
  # Ajuste da curva normal
  fit <- fitdist(vmean, "norm")
  
  # Preparando a curva teórica
  xmin <- min(vmean)
  xmax <- max(vmean)
  lnspc <- seq(xmin, xmax, length.out = 100)
  pdf_g <- dnorm(lnspc, mean = fit$estimate[1], sd = fit$estimate[2])
  
  # Armazenando os dados para a renderização
  df <- data.frame(vmean = vmean, lnspc = lnspc, pdf_g = pdf_g, sample_size = n)
  df_combined <- rbind(df_combined, df)
}

# Gerando o gráfico com renderização dinâmica
p <- ggplot(df_combined, aes(x = vmean)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "#0504aa", alpha = 0.7) +  # Calcula a densidade correta
  geom_line(aes(x = lnspc, y = pdf_g), color = "blue", size = 1) +
  labs(x = expression(bar(X)), y = expression(Prob(bar(X)))) +  # Ajustando os nomes dos eixos
  theme_minimal(base_size = 15) +
  theme(axis.title.x = element_text(size = 25, face = "bold", family = "serif"),
        axis.title.y = element_text(size = 25, face = "bold", family = "serif")) +
  transition_states(sample_size, transition_length = 2, state_length = 1) +
  ease_aes('linear')

# Renderizando a animação
#animate(p, duration = 10, fps = 20)
# Renderizando e salvando a animação como GIF

# Salvar e renderizar a animação como GIF
library(gifski)

anim <- animate(p, duration = 10, fps = 20, renderer = gifski_renderer())
anim_save("animacao.gif", anim)

Aqui está a animação gerada:

Notamos que a distribuição da média amostral tende à distribuição Normal quando aumentamos o tamanho da amostra. Notem que a média da distribuição converge para a média da população, isto é, \(E[\bar{X}]=\mu\).

Teorema Central do Limite: Exemplo

Exemplo: Seja a variável aleatória com distribuição de probabilidade: P(X=3)=0,4; P(X=6)=0,3; P(X=8)=0,3. Uma amostra com 40 observações é sorteada. Qual é a probabilidade de que a média amostral ser maior do que 5?

A média e o desvio padrão dessa população:

# Definindo os valores de X e P
X <- c(3, 6, 8)
P <- c(0.4, 0.3, 0.3)

# Inicializando as variáveis de Esperança e Esperança do quadrado
E <- 0
E2 <- 0

# Calculando a Esperança (E) e Esperança do quadrado (E2)
for (i in 1:length(X)) {
  E <- E + X[i] * P[i]
  E2 <- E2 + (X[i]^2) * P[i]
}

# Calculando a Variância (V)
V <- E2 - E^2

# Exibindo os resultados
cat("Esperança:", E, "\nVariância:", V, "\n")

Esperança: 5.4 
Variância: 4.44 

Vamos sortear várias amostras de tamanho n=40 e verificar qual a probabilidade da média dessa amostra ser maior do que 5.

library(ggplot2)

# Definindo os valores de X, P e variáveis
X <- c(3, 6, 8)
P <- c(0.4, 0.3, 0.3)
n <- 40
ns <- 1000  # número de simulações

# Armazenar a média amostral
vx <- numeric(ns)

# Simulações
for (s in 1:ns) {
  A <- sample(X, n, replace = TRUE, prob = P)  # Amostragem com reposição usando as probabilidades
  vx[s] <- mean(A)  # Calcula a média amostral
}

# Plotar o histograma
g_1 <- ggplot(data.frame(vx), aes(x = vx)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "#0504aa", alpha = 0.7) +
  labs(x = "Mean of X", y = "Density") +  # Substituindo expression() por strings comun
  theme_minimal(base_size = 15)
g_1p <- ggplotly(g_1)

Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.
ℹ The deprecated feature was likely used in the ggplot2 package.
  Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.

g_1p

# Exibir as médias
cat("Média das amostras:", mean(vx), "\nMédia da população:", E, "\n")

Média das amostras: 5.385775 
Média da população: 5.4 

A probabilidade de ser maior do que 5:

# Inicializando a variável para contar quantas médias são maiores que 5
nmaior <- 0

# Loop para verificar quais médias são maiores que 5
for (i in 1:length(vx)) {
  if (vx[i] > 5) {
    nmaior <- nmaior + 1
  }
}

# Calculando a probabilidade
nmaior <- nmaior / length(vx)

# Exibindo o resultado
cat("Probabilidade de ser maior do que 5:", nmaior, "\nValor teórico: 0.88\n")

Probabilidade de ser maior do que 5: 0.862 
Valor teórico: 0.88

Exemplo: Seja a variável aleatória com distribuição de probabilidade: P(X=1)=0,4; P(X=5)=0,3; P(X=10)=0,3. Uma amostra com 100 observações é sorteada. Qual é a probabilidade de que a média amostral ser menor do que 5?

Vamos sortear várias amostras de tamanho n=100 e verificar qual a probabilidade da média dessa amostra ser maior do que 5.

# Definindo os valores de X e P
X <- c(1, 5, 10)
P <- c(0.4, 0.3, 0.3)

# Inicializando as variáveis de Esperança e Esperança do quadrado
E <- 0
E2 <- 0

# Calculando a Esperança (E) e Esperança do quadrado (E2)
for (i in 1:length(X)) {
  E <- E + X[i] * P[i]
  E2 <- E2 + (X[i]^2) * P[i]
}

# Calculando a Variância (V)
V <- E2 - E^2

# Exibindo os resultados
cat("Esperança:", E, "\nVariância:", V, "\n")

Esperança: 4.9 
Variância: 13.89 

# Simulações
n <- 100
ns <- 1000
vx <- numeric(ns)

# Realizando as simulações
for (s in 1:ns) {
  A <- sample(X, n, replace = TRUE, prob = P)  # Amostragem com reposição usando as probabilidades
  vx[s] <- mean(A)
}

# Calculando P(X < 5)
nmenor <- sum(vx < 5) / length(vx)

# Exibindo o resultado
cat("P(X < 5) =", nmenor, "\nValor teórico: 0.60\n")

P(X < 5) = 0.594 
Valor teórico: 0.60

library(ggplot2)

# Definindo os valores de X e P
X <- c(1, 5, 10)
P <- c(0.4, 0.3, 0.3)
E <- sum(X * P)  # Esperança já calculada

# Simulações
n <- 100
ns <- 1000
vx <- numeric(ns)

# Realizando as simulações
for (s in 1:ns) {
  A <- sample(X, n, replace = TRUE, prob = P)  # Amostragem com reposição
  vx[s] <- mean(A)
}

# Plotando o histograma com ggplot2
p <- ggplot(data.frame(vx), aes(x = vx)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "#0504aa", alpha = 0.7) +
  labs(x = "Mean of X", y = "Density") +  # Substituindo expression() por strings comun
  theme_minimal(base_size = 15)
ggplotly(p)

# Exibindo a média amostral e a média populacional
cat("Média das amostras:", mean(vx), "\nMédia da população:", E, "\n")

Média das amostras: 4.90764 
Média da população: 4.9