Frank Harrell a făcut acest comentariu lămuritor la BlueSky zilele trecute:
„În centrul eșecului CLT (teorema limită centrală) este necesitatea ca media și abaterea standard să fie independente. Cu distribuțiile asimetrice, acestea sunt departe de a fi independente. N=50.000 poate fi mult prea mic pentru ca CLT să funcționeze suficient de bine.”
La fel ca și alți oameni, mi-a trecut imediat prin minte că wow, 50.000 este ceva mai mult de 30. Eram conștient, desigur, că teorema centrală a limitei poate avea nevoie de dimensiuni destul de mari ale eșantionului pentru ca mediile din distribuțiile oblice să devină „aproximativ normale”, pentru o definiție pragmatică a „aproximativ”. Dar nu mă gândisem la atât de mare. Așa că am făcut câteva simulări pentru a ajuta la imprimarea acestei impresii în mintea mea.
Am început cu construirea asta sim_clt() funcția în R. Ei bine, de fapt, a început mult mai simplu decât asta, dar iată versiunea finală care îmi permite să fac câteva lucruri pe care vreau să le fac cu aceste simulări.
Sarcina de bază a acestei funcții este de a lua o populație (generată separat), de a extrage un eșantion de o dimensiune specificată din acea populație și de a calcula media eșantionului. Apoi, faceți acest lucru de mai multe ori (implicit 10.000) pentru a examina distribuția acelor mijloace. Ieșirile includ o diagramă qq-normal; statistica Anderson-Darling generată ca și cum ar fi folosită pentru un test de normalitate (cu cât valoarea sa este mai mare în acest caz, cu atât distribuția mediei eșantionului este mai puțin normală); și proporția intervalelor de încredere, construite presupunând normalitatea mediei eșantionului (adică afirmația de bază a teoremei limitei centrale), care conțin de fapt media adevărată.
library(nortest)
library(tidyverse)
library(actuar)
library(glue)
library(scales)
set.seed(123)
today_reps <- 10000
N <- 1e6
sim_clt <- function(
x,
n = 30,
reps = today_reps,
replace = TRUE,
plot = TRUE,
conf = 0.95,
...
) {
true_mean <- mean(x)
samples <- replicate(
reps,
sample(x = x, size = n, replace = replace),
simplify = FALSE
)
means <- sapply(samples, mean)
covered <- sapply(samples, function(s) {
# rely on asymptotic normality to estimate a confidence interval and check
# for coverage for each sample
se <- stats::sd(s) / sqrt(length(s))
ci <- mean(s) + c(-1, 1) * qnorm((1 - conf) / 2 + conf) * se
ci(1) <= true_mean & true_mean <= ci(2)
})
if (plot) {
qqnorm(means, bty = "l", ...)
qqline(means, col = "steelblue")
grid()
}
return(list(
ad_stat = nortest::ad.test(means)$statistic,
coverage = mean(covered)
))
}
Utilizarea funcției este cam așa. În această linie de cod, creez 10.000 de mostre simulate dintr-o populație de un milion de puncte dintr-o distribuție normală standard. Fiecare eșantion are doar 5 dimensiuni:
sim_clt(rnorm(N), n = 5, main = "Mean from normal dist, n = 5")
Asta îmi oferă informații că statistica Anderson Darling a celor 10.000 de eșantionuri rezultate înseamnă 0,45, că un interval de încredere de 95% creat cu o distribuție normală și folosind deviația standard a eșantionului conține media adevărată 87% din timp, iar acest grafic QQ-normal:
Folosesc destul de deliberat o distribuție normală pentru statistica medie a eșantionului pentru a produce intervalul de încredere, mai degrabă decât la distribuție. Distribuția t ar fi adecvată în acest caz particular (unde se știe că populația este distribuită normal) și ar oferi un interval de încredere mai larg, cu o acoperire mult mai apropiată de 95%. Dar mă prefac că trebuie să mă bazez în schimb pe normalitatea asimptotică. Scopul meu este să testez cât de bună (în termeni pragmatici) este ipoteza că statistica în sine (adică media eșantionului) este distribuită în mod normal datorită teoremei limitei centrale. Inferența folosind varianța eșantionului ca estimare a varianței populației va funcționa în construirea unei erori standard pe măsură ce dimensiunea eșantionului devine mai mare; folosirea la distribuție atunci când se întâmplă să știm că se aplică ar fi înșelăciune.
De aceea, acoperirea este de doar 87%, nu de 95% pe care am avea-o dacă am fi construit un interval de încredere cu la distribuție.
Un lucru pe care îl vedem din graficul QQ-normal de mai sus este că linia dreaptă ne spune că distribuția mediei eșantionului, chiar și din acest eșantion minuscul de 5, este distribuită în mod normal. Desigur că știm asta teoretic. Se dovedește a fi destul de diferit cu o dimensiune a eșantionului de 5 față de o distribuție anormală, log-normală:
Într-adevăr, chiar și atunci când creștem dimensiunea eșantionului la 100, tot nu obținem o distribuție normală satisfăcătoare a eșantionului înseamnă:
Cele două comploturi au fost făcute cu asta:
lnpop <- exp(rnorm(N))
sim_clt(lnpop, n = 5, main = "Mean from log-normal dist, n = 5")
sim_clt(lnpop, n = 100, main = "Mean from log-normal dist, n = 100")
Este în regulă dacă mă uit la un experiment la un moment dat, dar ceea ce îmi doream cu adevărat era o gamă de distribuții diferite, încercate sistematic pentru diferite dimensiuni de eșantion, fiecare combinație rulată de 10.000 de ori. Așadar, asta face următoarea bucată de cod: definește unele dimensiuni de eșantion, face o listă de funcții care generează descrieri și le etichetează și rulează simularea pentru fiecare combinație.
Distribuțiile pe care le folosesc sunt:
- distribuție normală standard (adică medie 0, abatere standard 1).
- standard (adică media 0, abaterea standard 1 pentru distribuția normală subiacentă) log normal
- o varianță log-normală mai mare (distribuția normală de bază are medie 0, abaterea standard de 2)
- o distribuție Cauchy cu varianță infinită – acesta este un caz cunoscut în care teorema limită centrală nu funcționează și indiferent de dimensiunea eșantionului, distribuția mediei nu tinde niciodată spre normal. Acesta este extremul opus al cazului în care distribuția populației este un simplu normal.
- o distribuție Pareto cu forma de 2,1 și scară de 1 – acesta a fost ales ca un caz provocator, dar nu imposibil pentru teorema limită centrală
- o distribuție exponențială standard
- un amestec 60/40 de distribuție exponențială standard și distribuție normală standard
- un amestec de 90/10 dintr-o distribuție normală standard și o distribuție normală anormală – pentru a testa cazul unui amestec de două normale, care, în general, l-am găsit un caz interesant și uneori supărător.
Iată codul pentru toate acestea:
ns <- c(5:30, 50, 100, 300, 1000, 2000, 5000, 10000, 20000, 40000)
sims <- list()
pop_dists <- list(
"rnorm(N)" = function() rnorm(N),
"exp(rnorm(N))" = function() exp(rnorm(N)),
"exp(rnorm(N, sd = 2))" = function() exp(rnorm(N, sd = 2)),
"rcauchy(N)" = function() rcauchy(N),
"rpareto(N, shape = 2.1, scale = 1)" = function() {
rpareto(N, shape = 2.1, scale = 1)
},
"rexp(N)" = function() rexp(N),
"c(rexp(N*0.6), exp(rnorm(N*0.4)))" = function() {
c(rexp(N * 0.6), exp(rnorm(N * 0.4)))
},
"c(rnorm(N*0.9,0,1), rnorm(N*0.1,3,3))" = function() {
c(rnorm(N * 0.9, 0, 1), rnorm(N * 0.1, 3, 3))
}
)
all_sims <- lapply(names(pop_dists), function(lbl) {
pop <- pop_dists((lbl))()
results <- lapply(ns, function(n) sim_clt(pop, n = n, plot = FALSE))
tibble(
n = ns,
ad_stat = sapply(results, `((`, "ad_stat"),
coverage = sapply(results, `((`, "coverage"),
pop_dist = lbl
)
}) |>
bind_rows()
Acum, evident, nu ne vom uita la toate diagramele QQ-normale ale tuturor acestor distribuții și dimensiuni ale eșantionului, ci mai degrabă la unele statistici rezumative. Prima mea idee a fost să mă uit la statisticile Anderson-Darling, care măsoară cât de diferită este o distribuție de o distribuție de referință dată. Având același număr de repetări ale fiecărei simulări, am gândit că acestea ar fi comparabile; și, în special, le-am putea compara cu cazul în care populația de bază este ea însăși normală, așa că știm că, indiferent de dimensiunea eșantionului, media eșantionului va fi, de asemenea, normală.
Acest grafic următor arată această comparație, ordonând distribuțiile de la cel mai bine comportat (după standardele teoremei limită centrale) la cel mai rău (pe care știam dinainte că va fi Cauchy):
Deci, putem vedea deja că dimensiunile eșantionului pentru multe dintre aceste distribuții trebuie într-adevăr să fie de mii sau zeci de mii pentru ca eșantionul să fie la fel de bine comportat, în ceea ce privește o distribuție normală, ca în cel mai bun caz. Cele mai rele două cazuri sunt distribuția Pareto și varianța mai mare log-normală.
Cât de rău este acest lucru pragmatic? Folosim teorema limită centrală nu pentru că ne pasă direct de forma distribuției statisticii noastre, ci din cauza inferențelor pe care vrem să le tragem. Deci, funcția cu care am început această postare pe blog calculează și un interval de încredere bazat pe normalitatea asimptotică.
Deci aici vedem rezultatele. Practic, nu atât de grozav. Pentru ca intervalele noastre de încredere de 95% să aibă o acoperire de aproximativ 95%, în unele dintre aceste distribuții avem într-adevăr nevoie de mii sau zeci de mii de observații în eșantionul nostru:
Iată aceeași diagramă mărită la doar 80% și mai mult. Rețineți că media eșantionată din distribuția Cauchy nici măcar nu apare pe diagramă:
Iată codul care a produs acele diagrame:
normal_band <- all_sims |>
filter(pop_dist == "rnorm(N)") |>
pull(ad_stat) |>
quantile(probs = c(0.025, 0.975))
p1 <- all_sims |>
mutate(pop_dist = fct_reorder(pop_dist, ad_stat)) |>
ggplot(aes(x = n, y = ad_stat)) +
facet_wrap(~pop_dist, nrow = 2) +
annotate(
"rect",
xmin = min(ns),
xmax = Inf,
ymin = normal_band(1),
ymax = normal_band(2),
alpha = 0.5,
fill = "orange"
) +
geom_point() +
scale_x_log10(label = comma, breaks = ns(ns >= 30 | ns %in% c(5, 10, 20))) +
scale_y_log10(label = comma) +
theme(
panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text = element_text(face = "plain", size = 9)
) +
labs(
x = "Sample size",
y = "Anderson-Darling statistic",
subtitle = glue(
"Points represent average from {comma(today_reps)} simulations of given sample size. Shaded area covers 95% of values from a normal distribution."
),
title = "Increasing sample size and growing effectiveness of the central limit theorem - Anderson-Darling statistic"
)
print(p1)
p2 <- all_sims |>
mutate(pop_dist = fct_reorder(pop_dist, ad_stat)) |>
ggplot(aes(x = n, y = coverage)) +
facet_wrap(~pop_dist, nrow = 2) +
annotate(
"rect",
xmin = min(ns),
xmax = Inf,
ymin = 0.95,
ymax = 1,
alpha = 0.5,
fill = "orange"
) +
geom_point() +
scale_x_log10(label = comma, breaks = ns(ns >= 30 | ns %in% c(5, 10, 20))) +
scale_y_continuous(label = percent) +
theme(
panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text = element_text(face = "plain", size = 9)
) +
labs(
x = "Sample size",
y = "Coverage of a 95% confidence interval based on t distribution",
subtitle = glue(
"Points represent average from {comma(today_reps)} simulations of given sample size. Shaded area shows 95% and higher, as desired."
),
title = "Increasing sample size and growing effectiveness of the central limit theorem - confidence interval coverage"
)
print(p2)
p3 <- p2 +
labs(
y = "Coverage of a 95% confidence interval based on t distributionnY axis truncated to only start at 80%; some points excluded because of that."
) +
coord_cartesian(ylim = c(0.8, 1))
print(p3)
Deci asta e cu adevărat. Desigur, profesorul Harrell a avut dreptate așa cum știam eu că va avea. Și nu a fost greu să vină cu câteva distribuții distorsionate care au fost problematice. Nici măcar amestecul a două distribuții normale nu a funcționat genial. Deci, da, aveți grijă cu dimensiunile eșantionului atunci când extrageți din distribuțiile dificile ale populației de bază.
