Manualele spun că LOOCV are cea mai scăzută părtinire, dar cea mai mare variație, comparativ cu 10 și 5 ori. Am codificat un CV K-Fold de la zero pentru a învăța să îl testați pe date simulate 🔍📊 — și la n=1000 rezistă. La n=100? Nu atât de mult. 🤔
Imaginea de mai sus a fost generată prin chatGPT. Am încărcat tot textul acestei postări pe blog și i-a cerut să genereze un desen animat. Foarte impresionant! În trecut, era o greșeală de ortografie și un farfurie de text, dar acum cuvinte coezive pe imagine. Doar wow.
Motivații
Validarea încrucișată este un pas atât de crucial în învățarea automată (și în metodele tradiționale), încât în prezent este încorporat în modele de sklearn sau tidymodele ușor de utilizat, fără a fi nevoie să construim unul de la zero. Ca și în cealaltă experiență de învățare a mea, cel mai bun mod de a învăța conceptul (altul decât învățarea conceptului 🤣) este să îl codific de la zero și să vezi cum funcționează! În K-Fold CV, datele de antrenament sunt împărțite în K bucăți; modelul este antrenat de K ori, de fiecare dată ținând o bucată diferită. Performanța este mediată pentru toate pliurile K, oferind o estimare mai stabilă. Un caz special este Leave-One-Out CV (LOOCV), în care fiecare observație individuală servește ca set propriu de validare. Este minuțios, dar costisitor din punct de vedere informatic. Mi s-a spus că, părtinire LOOCV < de 10 ori < de 5 ori; în timp ce varianța LOOCV > 10 ori > 5 ori. Este adevărat? De asemenea, ce se întâmplă cu repetările, asta chiar reduce varianța? Să le verificăm.
Obiective:
Simulați date
library(tidyverse) set.seed(1) n <- 1000 x <- rnorm(n) w <- rnorm(n) y <- 0.5*x^2 + -0.5*w + 0.3*w*x + rnorm(n) df <- tibble(x,y,w) idx <- sample(1:n, size=0.8*n) train <- df(idx, ) test <- df(-idx, )
Codul de mai sus simulează un set de date cu 1000 de observații, unde variabila răspuns y este generat pe baza unui proces cunoscut de generare a datelor care implică predictori x şi w. Setul de date este apoi împărțit într-un set de antrenament (80%) și un set de testare (20%). Să vizualizăm.
df |> mutate(w_cut = cut_interval(w, n=5)) |> ggplot(aes(x=x, y=y, color=w_cut, group=w_cut)) + geom_point(alpha=0.5) + theme_bw() + geom_smooth(method = "gam", se=F)
Uau, vizualizare foarte interesantă în care relațiile nu sunt cu siguranță liniare aici. Este o formă de interacțiune între x şi w. Să vedem dacă putem recupera procesul de bază de generare a datelor folosind validarea încrucișată K-Fold.
Validare încrucișată K-Fold de la zero
folds <- 5
segment_portion <- nrow(train)/folds
formula_list <- list(as.formula("y~x"),as.formula("y~I(x^2)"),as.formula("y~I(x^2)+w+w:x"),as.formula("y~I(x^3)+w+w:x"),
as.formula("y~w:x"),as.formula("y~w"),as.formula("y~x+w+x:w"),as.formula("y~I(x^2)+w:x"),
as.formula("y~I(x^2)+w"))
cv_log <- tibble()
for (formula in formula_list) {
print(formula)
predict_log <- y_log <- vector(mode="numeric",length=segment_portion*folds)
start <- 1
end <- segment_portion
for (fold in 1:folds) {
val_i <- train(start:end,)
train_i <- train(-c(start:end),)
model_i <- lm(formula,train_i)
predict_i <- predict(model_i, val_i)
predict_log(start:end) <- predict_i
y_log(start:end) <- val_i$y
start <- end + 1
end <- start + segment_portion - 1
}
val_df <- tibble(predict=predict_log,y=y_log) |>
mutate(formula=deparse(formula))
cv_log <- cv_log |>
bind_rows(val_df)
}
## y ~ x
## y ~ I(x^2)
## y ~ I(x^2) + w + w:x
## y ~ I(x^3) + w + w:x
## y ~ w:x
## y ~ w
## y ~ x + w + x:w
## y ~ I(x^2) + w:x
## y ~ I(x^2) + w
Bine, ceea ce am făcut mai sus este o implementare manuală a validării încrucișate K-Fold. Parcurgem fiecare formulă din lista noastră și, pentru fiecare formulă, împărțim datele de antrenament în 5 ori. Pentru fiecare pliu, antrenăm modelul pe celelalte 4 falduri și îl validăm pe pliul curent. Stocăm predicțiile și valorile reale pentru evaluare ulterioară.
Practic, vrem să vedem care formulă are cel mai mic RMSE în pliuri. Să calculăm asta în continuare. Din formula DGP, știm că cel mai bun model ar trebui să fie y~I(x^2)+w+w:x. Să vedem dacă putem recupera asta folosind K-Fold CV.
Evaluarea RMSE
cv_log |> group_by(formula) |> summarize(rmse = sqrt(mean((y-predict)^2))) |> arrange(rmse) |> mutate(rmse = format(rmse, digits = 8)) ## # A tibble: 9 × 2 ## formula rmse #### 1 y ~ I(x^2) + w + w:x 1.0397338 ## 2 y ~ I(x^2) + w 1.1008693 ## 3 y ~ I(x^2) + w:x 1.1607225 ## 4 y ~ I(x^2) 1.2144857 ## 5 y ~ x + w + x:w 1.2826723 ## 6 y ~ I(x^3) + w + w:x 1.2912046 ## 7 y ~ w 1.3578304 ## 8 y ~ w:x 1.3732478 ## 9 y ~ x 1.4451807
Aici funcția noastră de pierdere este RMSE de atunci y este o dată continuă și încercăm să anticipăm asta. Formula cu cel mai mic RMSE este într-adevăr y~I(x^2)+w+w:xcare se potrivește cu procesul de bază de generare a datelor. OK, cel puțin, acum suntem capabili să recuperăm DGP-ul de bază folosind validarea încrucișată de 5 ori. Dar există o diferență între 5 ori, 10 ori sau chiar LOOCV? Dacă există o diferență, cum o evaluăm? În trecut, am putut evalua părtinirea și varianța pe baza unui ATE adevărat, dar ce naiba este un RMSE adevărat !?!
Pentru a verifica dacă afirmația manuală (bias LOOCV < 10 ori < 5 ori; varianța LOOCV > 10 ori > 5 ori) se menține, am rulat o mică simulare cu ajutorul lui Claude Sonnet 5. Deoarece controlăm procesul de generare a datelor, putem compara „formula corectă” (presupunând că formula corectă are cel mai mic RMSE de date simulat față de cel mai mic RMSE de mai sus față de un RMSE simulat de mai sus) estimat dintr-un set mare de teste (n=10000) — suficient de mare, prin legea numerelor mari, pentru a fi tratat ca adevăr de bază. Valoarea mediei între simulări oferă o părtinire (cât de departe este CV-ul de eroarea adevărată) și varianță (cât de mult variază estimarea CV-ului de la un eșantion la altul) pentru fiecare metodă. Este legal acest lucru? 🤔 Dacă afirmația manuală este corectă, ar trebui să putem observa părtinire LOOCV < de 10 ori < de 5 ori; varianţă LOOCV > de 10 ori > de 5 ori. Să vedem dacă putem observa asta în simularea de mai jos.
Comparați modelele candidate
# set.seed(1)
# k-fold CV RMSE for a given formula and dataset (k = n gives LOOCV)
cv_rmse <- function(data, formula, k) {
n <- nrow(data)
folds <- sample(rep(1:k, length.out = n))
preds <- numeric(n)
for (i in 1:k) {
train_i <- data(folds != i, )
val_i <- data(folds == i, )
model_i <- lm(formula, train_i)
preds(folds == i) <- predict(model_i, val_i)
}
sqrt(mean((data$y - preds)^2))
}
# "true" RMSE: fit on train, evaluate on a large fresh draw from the DGP
true_rmse <- function(train, formula, n_test = 10000) {
x <- rnorm(n_test); w <- rnorm(n_test)
y <- 0.5*x^2 - 0.5*w + 0.3*w*x + rnorm(n_test)
test <- tibble(x, y, w)
model <- lm(formula, train)
sqrt(mean((test$y - predict(model, test))^2))
}
formula_true <- as.formula("y ~ I(x^2) + w + w:x")
n_sim <- 500
n_train <- 100
results <- vector("list", n_sim)
for (s in 1:n_sim) {
x <- rnorm(n_train); w <- rnorm(n_train)
y <- 0.5*x^2 - 0.5*w + 0.3*w*x + rnorm(n_train)
train_s <- tibble(x, y, w)
results((s)) <- tibble(
sim = s,
true_err = true_rmse(train_s, formula_true),
loocv = cv_rmse(train_s, formula_true, k = n_train),
cv5 = cv_rmse(train_s, formula_true, k = 5),
cv10 = cv_rmse(train_s, formula_true, k = 10)
)
}
sim_df <- bind_rows(results)
sim_long <- sim_df |>
pivot_longer(cols = c(loocv, cv5, cv10), names_to = "method", values_to = "cv_estimate")
sim_long |>
group_by(method) |>
summarize(
mean_cv_estimate = mean(cv_estimate),
mean_true_error = mean(true_err),
bias = mean(cv_estimate - true_err),
variance = var(cv_estimate),
.groups = "drop"
) |>
arrange(bias) |>
mutate(variance = format(variance, digit = 8))
## # A tibble: 3 × 5
## method mean_cv_estimate mean_true_error bias variance
##
## 1 loocv 1.00 1.00 0.000703 0.00053227867
## 2 cv10 1.00 1.00 0.000945 0.00053192050
## 3 cv5 1.00 1.00 0.00124 0.00053176441
Uau, uitându-ne la rezultate, putem vedea că afirmația manuală se menține. LOOCV are cea mai mică părtinire, dar cea mai mare varianță. CV-ul de 10 ori este între ele, iar CV-ul de 5 ori are cea mai mare prejudecată, dar cea mai mică variație. Dar, am observat că a trebuit să ne creștem cifra la 8 pentru a vedea diferența de varianță, este într-adevăr minuscul. În principal pentru că n=1000 nostru este deja destul de mare, deci varianța este deja destul de mică. Dacă reducem n la 100, cum ar arăta?
## # A tibble: 3 × 5 ## method mean_cv_estimate mean_true_error bias variance #### 1 loocv 1.02 1.02 -0.00400 0.0063905148 ## 2 cv10 1.02 1.02 -0.00138 0.0065753049 ## 3 cv5 1.03 1.02 0.00244 0.0068098389
!?!?!?!?! 🤷♂️ cu n=100, euristica ordinului de prejudecată și varianță nu mai rezistă? De ce este asta? Nu știu. Dacă o faci, te rog să-mi spui. Am crescut chiar și n_sim-ul la 1000, dar tot același model. Interesant cum a trebuit să împing n până la 850 pentru a observa din nou ordinea de variație a manualului, deși din nou este destul de minuscul.
Verificați pe set de testare
test |>
mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
mutate(residual = y-predict) |>
mutate(res_square = residual^2) |>
pull(res_square) |>
mean() |>
sqrt()
## (1) 1.018281
Bine! RMSE pe setul de testare este destul de similar cu seturile noastre medii de validare! 🙌 Să vizualizăm y-ul prezis și real al modelului de pe setul de testare.
test |>
mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
ggplot(aes(x=predict, y=y)) +
geom_point(alpha=0.5) +
theme_bw() +
geom_smooth(method = "lm") +
labs(title="Predicted vs Actual y on Test Set", x="Predicted y", y="Actual y")
summary(lm(y ~ predict, data=test |> mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test))))
##
## Call:
## lm(formula = y ~ predict, data = mutate(test, predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"),
## train), test)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.98185 -0.61410 -0.05514 0.62883 2.28837
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.08581 0.08390 1.023 0.308
## predict 0.80885 0.07562 10.696 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.007 on 198 degrees of freedom
## Multiple R-squared: 0.3662, Adjusted R-squared: 0.363
## F-statistic: 114.4 on 1 and 198 DF, p-value: < 2.2e-16
Panta aici este 0,81, nu 1 – chiar și cu formula corectă, coeficienții sunt încă estimați dintr-un eșantion finit, zgomotos, deci predicția este într-adevăr „semnal adevărat + eroare de estimare”. Regresarea y pe o versiune zgomotoasă a lui însuși trage panta sub 1 – un efect cunoscut numit atenuare (aceeași idee ca și eroarea în variabile). Mai multe date de antrenament micșorează acea eroare de estimare, astfel încât panta ar trebui să se apropie de 1. Interceptarea rămâne aproape de 0, deoarece eroarea de estimare nu este părtinitoare sistematic într-o singură direcție – doar zgomotoasă – deci predicțiile încă se centrează corect în medie. Dacă acest lucru este adevărat, atunci dacă creștem n-ul nostru la 10000, ar trebui să vedem că panta este mai aproape de 1 și interceptarea mai aproape de 0. Să testăm această teorie.
set.seed(1)
n <- 10000
x <- rnorm(n)
w <- rnorm(n)
y <- 0.5*x^2 + -0.5*w + 0.3*w*x + rnorm(n)
df <- tibble(x,y,w)
idx <- sample(1:n, size=0.8*n)
train <- df(idx, )
test <- df(-idx, )
formula <- as.formula("y ~ I(x^2) + w + w:x")
folds <- 10
segment_portion <- nrow(train)/folds
predict_log <- y_log <- vector(mode="numeric",length=segment_portion*folds)
start <- 1
end <- segment_portion
for (fold in 1:folds) {
val_i <- train(start:end,)
train_i <- train(-c(start:end),)
model_i <- lm(formula,train_i)
predict_i <- predict(model_i, val_i)
predict_log(start:end) <- predict_i
y_log(start:end) <- val_i$y
start <- end + 1
end <- start + segment_portion - 1
}
val_df <- tibble(predict=predict_log,y=y_log)
val_df |>
summarize(rmse = sqrt(mean((y-predict)^2))) |>
arrange(rmse) |>
mutate(rmse = format(rmse, digits = 8))
## # A tibble: 1 × 1
## rmse
##
## 1 1.0126176
test |>
mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
mutate(residual = y-predict) |>
mutate(res_square = residual^2) |>
pull(res_square) |>
mean() |>
sqrt()
## (1) 0.9884315
test |>
mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test)) |>
ggplot(aes(x=predict, y=y)) +
geom_point(alpha=0.5) +
theme_bw() +
geom_smooth(method = "lm") +
labs(title="Predicted vs Actual y on Test Set", x="Predicted y", y="Actual y")
summary(lm(y ~ predict, data=test |> mutate(predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"), train), test))))
##
## Call:
## lm(formula = y ~ predict, data = mutate(test, predict = predict(lm(as.formula("y ~ I(x^2) + w + w:x"),
## train), test)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.08711 -0.67088 0.01239 0.69015 2.86620
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.02778 0.02527 1.099 0.272
## predict 0.96577 0.02369 40.762 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9884 on 1998 degrees of freedom
## Multiple R-squared: 0.454, Adjusted R-squared: 0.4538
## F-statistic: 1662 on 1 and 1998 DF, p-value: < 2.2e-16
Iată-l! RMSE pe setul de testare este destul de similar cu seturile noastre medii de validare! 🙌 Panta este acum 0,97, mult mai aproape de 1, iar interceptarea este 0,02, mult mai aproape de 0. Acest lucru confirmă teoria noastră că, cu mai multe date de antrenament, eroarea de estimare scade, ceea ce duce la predicții mai bune.
Oportunități de îmbunătățire
- aplicați validarea încrucișată de la zero la reglarea hiperparametrului, cum ar fi xgboost hyperparam, lasso/ridge/elasticnet/glmnet
- aflați despre reeșantionarea imbricată care utilizează un strat suplimentar de reeșantionare care separă activitățile de reglare de procesul utilizat pentru a estima eficacitatea modelului. vezi aici
- învăţaţi să codificaţi de la zero cv grupate şi stratificate
- aflați puțin mai multe despre CV-urile din serii de timp
- trebuie să încercăm n_sim de 10000 și să vedem dacă ordinea de polarizare și varianță se menține pentru n = 100, deoarece simularea este destul de lentă, nu am rulat-o până la finalizare. Actualizare: am completat un n_sim de 10000 de n=100 și LOOCV a avut încă cea mai mică varianță.
Lecții învățate
- Am învățat că trebuie să folosim
Ipentru a indica faptul că dorim să includem termeni polinomi în formulă, altfel R îi va interpreta ca termeni de interacțiune. - Învăţat
deparseeste util pentru a converti un obiect formulă într-un șir de caractere în scopuri de înregistrare.
Daca iti place acest articol:
