From d85f14420c3cc5a9278fa60fa131687c32a99b10 Mon Sep 17 00:00:00 2001 From: tereom Date: Wed, 20 Nov 2024 15:06:16 -0600 Subject: [PATCH] Cambiar a nuevo pipe --- 10-bootstrap-parametrico.Rmd | 30 ++++++------- 12-mas-hipotesis.Rmd | 81 ++++++++++++++++++++---------------- 14-intro-bayesiana.Rmd | 8 ++-- 3 files changed, 63 insertions(+), 56 deletions(-) diff --git a/10-bootstrap-parametrico.Rmd b/10-bootstrap-parametrico.Rmd index f65daba..203afff 100644 --- a/10-bootstrap-parametrico.Rmd +++ b/10-bootstrap-parametrico.Rmd @@ -60,7 +60,7 @@ Ahora optimizamos (revisa que el método converge): ```{r} res <- optim(c(0, 0.5), log_p, control = list(fnscale = -1, maxit = 1000), method = "Nelder-Mead") res$convergence -est_mle <- tibble(parametro = c("media", "sigma"), estimador = res$par) %>% +est_mle <- tibble(parametro = c("media", "sigma"), estimador = res$par) |> column_to_rownames(var = "parametro") ``` @@ -94,7 +94,7 @@ log_p_boot <- crear_log_p(muestra_bootstrap) res_boot <- optim(c(0, 0.5), log_p_boot, control = list(fnscale = -1, maxit = 1000), method = "Nelder-Mead") res_boot$convergence -est_mle_boot <- tibble(parametro = c("media", "sigma"), estimador = res_boot$par) %>% +est_mle_boot <- tibble(parametro = c("media", "sigma"), estimador = res_boot$par) |> column_to_rownames(var = "parametro") est_mle_boot ``` @@ -129,15 +129,15 @@ reps_boot Ya ahora podemos estimar error estándar: ```{r} -error_est <- reps_boot %>% group_by(parametro) %>% +error_est <- reps_boot |> group_by(parametro) |> summarise(ee_boot = sd(estimador_boot)) error_est ``` Así que nuestra estimación final sería: ```{r} -bind_cols(est_mle, error_est) %>% - mutate(across(where(is.numeric), round, 3)) %>% +bind_cols(est_mle, error_est) |> + mutate(across(where(is.numeric), round, 3)) |> select(parametro, estimador, ee_boot) ``` Si usamos la rutina estándar de `R` (dejaremos para después explicar cómo @@ -170,7 +170,7 @@ log_p <- crear_log_p(muestra) # máxima verosimilitud res <- optim(c(0, 0.5), log_p, control = list(fnscale = -1, maxit = 1000), method = "Nelder-Mead") res$convergence -est_mle <- tibble(parametro = c("media", "sigma"), estimador = res$par) %>% +est_mle <- tibble(parametro = c("media", "sigma"), estimador = res$par) |> column_to_rownames(var = "parametro") est_mle ``` @@ -199,7 +199,7 @@ las estimaciones individuales - en algunos casos pueden estar correlacionadas. P examinar este comportamiendo visualizando las replicaciones bootstrap ```{r} -ggplot(reps_boot %>% pivot_wider(names_from = parametro, values_from = estimador_boot), +ggplot(reps_boot |> pivot_wider(names_from = parametro, values_from = estimador_boot), aes(x = media, y = sigma)) + geom_point(alpha = 0.5) + coord_equal() ``` @@ -315,11 +315,11 @@ Separamos por Cena y Comida, dividimos entre número de personas y probamos ajustando un modelo para cada horario: ```{r} -propinas <- read_csv("data/propinas.csv") %>% +propinas <- read_csv("data/propinas.csv") |> mutate(cuenta_persona = cuenta_total / num_personas) -propinas_mle <- propinas %>% - group_by(momento) %>% - summarise(est_mle = list(tidy(MASS::fitdistr(cuenta_persona, "lognormal")))) %>% +propinas_mle <- propinas |> + group_by(momento) |> + summarise(est_mle = list(tidy(MASS::fitdistr(cuenta_persona, "lognormal")))) |> unnest(est_mle) propinas_mle ``` @@ -329,11 +329,11 @@ para ver cómo calcular media y desviación estándar de las distribuciones orig Ahora verificamos el ajuste: ```{r} -g_1 <- ggplot(propinas %>% filter(momento == "Cena"), aes(sample = cuenta_persona)) + +g_1 <- ggplot(propinas |> filter(momento == "Cena"), aes(sample = cuenta_persona)) + geom_qq(dparams = list(mean = propinas_mle$estimate[1], sd = propinas_mle$estimate[2]), distribution = stats::qlnorm) + ylim(c(0, 20)) + geom_abline() + labs(subtitle = "Cena") -g_2 <- ggplot(propinas %>% filter(momento == "Comida"), aes(sample = cuenta_persona)) + +g_2 <- ggplot(propinas |> filter(momento == "Comida"), aes(sample = cuenta_persona)) + geom_qq(dparams = list(mean = propinas_mle$estimate[3], sd = propinas_mle$estimate[4]), distribution = stats::qlnorm) + ylim(c(0, 20)) + geom_abline() + labs(subtitle = "Comida") @@ -447,9 +447,9 @@ rep_boot <- function(rep, simular, crear_log_p, pars, n){ } set.seed(8934) reps_boot <- map(1:500, ~ rep_boot(.x, simular_modelo, crear_log_p, est_mle, - n = length(muestra))) %>% + n = length(muestra))) |> bind_rows -reps_boot %>% mutate(across(everything(), round, 2)) %>% head() +reps_boot |> mutate(across(everything(), round, 2)) |> head() ``` El optimizador encontró resultados que no tienen sentido: diff --git a/12-mas-hipotesis.Rmd b/12-mas-hipotesis.Rmd index a544bfd..aec2be9 100644 --- a/12-mas-hipotesis.Rmd +++ b/12-mas-hipotesis.Rmd @@ -82,11 +82,11 @@ Obtenemos: ```{r} set.seed(29) -muestra_edomex <- read_csv("data/enlace.csv") %>% - filter(estado == "ESTADO DE MEXICO") %>% +muestra_edomex <- read_csv("data/enlace.csv") |> + filter(estado == "ESTADO DE MEXICO") |> sample_n(180) -resumen <- muestra_edomex %>% - summarise(media = mean(mate_6), s = sd(mate_6), n = n()) %>% +resumen <- muestra_edomex |> + summarise(media = mean(mate_6), s = sd(mate_6), n = n()) |> mutate(ee = s / sqrt(n)) resumen ``` @@ -95,10 +95,11 @@ La hipótesis nula es que la media poblacional del Estado de México es igual a 454. Calculamos el valor-$p$ usando la prueba de Wald: ```{r} -dif <- (resumen %>% pull(media)) - 454 -ee <- resumen %>% pull(ee) +dif <- (resumen |> pull(media)) - 454 +ee <- resumen |> pull(ee) w <- dif / ee -p <- 2 * (1 - pnorm(abs(w))) + +p <- 2 * (1 - pt(abs(w), 179)) p ``` y vemos que esta muestra es consistente con la media nacional. No tenemos @@ -318,41 +319,47 @@ Supongamos por ejemplo que los datos que obtenemos son: ```{r, echo = FALSE, include = FALSE} set.seed(31521) -datos_clasif <- tibble(x = rep(1, 135), y = rep(1, 135)) %>% - bind_rows(tibble(x = rep(0, 65), y = rep(1, 65))) %>% - bind_rows(tibble(x = rep(1, 40), y = rep(0, 40))) %>% - bind_rows(tibble(x = rep(1 ,100), y = rep(1, 100))) %>% - sample_n(300) %>% - mutate(caso = as.character(1: 300)) %>% +datos_clasif <- tibble(x = rep(1, 135), y = rep(1, 135)) |> + bind_rows(tibble(x = rep(0, 65), y = rep(1, 65))) |> + bind_rows(tibble(x = rep(1, 40), y = rep(0, 40))) |> + bind_rows(tibble(x = rep(1 ,100), y = rep(1, 100))) |> + sample_n(300) |> + mutate(caso = as.character(1: 300)) |> select(caso, x, y) ``` ```{r} -datos_clasif %>% head +datos_clasif |> head() ``` Como explicamos arriba, nos interesa la diferencia. Calculamos $d$: ```{r} -datos_clasif <- datos_clasif %>% +datos_clasif <- datos_clasif |> mutate(d = x - y) -datos_clasif %>% head +datos_clasif |> head() + + +datos_clasif |> + summarise(sd_x = sd(x), + sd_y = sd(y), + sd_d = sd(d)) ``` Y ahora calculamos la media de $d$ (y tasa de correctos de cada clasificador:) ```{r} medias_tbl <- - datos_clasif %>% summarise(across(where(is.numeric), mean, .names = "{col}_hat")) + datos_clasif |> summarise(across(where(is.numeric), mean, .names = "{col}_hat")) d_hat <- pull(medias_tbl, d_hat) medias_tbl ``` Ahora necesitamos calcular el error estándar. Como explicamos arriba, hacemos ```{r} -ee <- datos_clasif %>% - mutate(d_hat = mean(d)) %>% - mutate(dif_2 = (d - d_hat)) %>% - summarise(ee = sd(dif_2) / sqrt(n())) %>% +ee <- datos_clasif |> + mutate(d_hat = mean(d)) |> + mutate(dif_2 = (d - d_hat)) |> + summarise(ee = sd(dif_2) / sqrt(n())) |> pull(ee) ee ``` @@ -362,7 +369,7 @@ Y ahora podemos calcular la estadística $W$ y el valor p correspondiente: ```{r} w <- d_hat / ee valor_p <- 2 * (1 - pnorm(abs(w))) -c(w = w, valor_p = valor_p) %>% round(3) +c(w = w, valor_p = valor_p) |> round(3) ``` @@ -421,7 +428,7 @@ hipótesis nula. Esto se explica en la siguiente gráfica: log_verosim <- function(p){ 75 * log(p) + (120 - 75) * log(1 - p) } -verosim_tbl <- tibble(p = seq(0.4, 0.7, 0.01)) %>% +verosim_tbl <- tibble(p = seq(0.4, 0.7, 0.01)) |> mutate(log_verosim = log_verosim(p)) ggplot(verosim_tbl, aes(x = p, y = log_verosim)) + geom_line() + @@ -479,7 +486,7 @@ lambda <- function(n, x, p_0 = 0.5){ lambda } lambda_obs <- lambda(n_volados, 75, 0.5) -sims_tbl <- tibble(sim_x = simulados_nula) %>% +sims_tbl <- tibble(sim_x = simulados_nula) |> mutate(lambda = map_dbl(sim_x, ~ lambda(n_volados, .x, p_0 = 0.5))) ggplot(sims_tbl, aes(x = lambda)) + geom_histogram(binwidth = 0.7) + @@ -514,7 +521,7 @@ crear_log_p <- function(x){ # crear log verosim para dos muestras normales independientes. log_p <- function(params){ mu <- params[1] - log_vero <- dnorm(x, mean = mu, sd = 1, log = TRUE) %>% sum + log_vero <- dnorm(x, mean = mu, sd = 1, log = TRUE) |> sum() log_vero } } @@ -536,7 +543,7 @@ bajo la nula ```{r, cache = TRUE} sims_nula <- map(1:10000, ~ rnorm(n_muestra, 8, 1)) lambda_nula_sim <- map_dbl(sims_nula, ~ lambda_calc(.x, crear_log_p)) -tibble(lambda = lambda_nula_sim) %>% +tibble(lambda = lambda_nula_sim) |> ggplot(aes(x = lambda)) + geom_histogram() + geom_vline(xintercept = lambda, colour = "red") ``` @@ -553,9 +560,9 @@ resultan ser aproximadamente **$\chi$-cuadrada con 1 grado de libertad** (ji-cua checar para el último ejemplo: ```{r} -teorica <- tibble(x = seq(0.1, 10, 0.01)) %>% +teorica <- tibble(x = seq(0.1, 10, 0.01)) |> mutate(f_chi_1 = dchisq(x, df = 1)) -tibble(lambda = lambda_nula_sim) %>% +tibble(lambda = lambda_nula_sim) |> ggplot() + geom_histogram(aes(x = lambda, y = ..density..), binwidth = 0.1) + geom_line(data = teorica, aes(x = x, y = f_chi_1), colour = "red") ``` @@ -564,7 +571,7 @@ O mejor, con una gráfica de cuantiles de las simulaciones vs la téorica: ```{r} -tibble(lambda = lambda_nula_sim) %>% +tibble(lambda = lambda_nula_sim) |> ggplot(aes(sample = lambda)) + geom_qq(distribution = stats::qchisq, dparams = list(df = 1)) + geom_qq_line(distribution = stats::qchisq, dparams = list(df = 1)) @@ -606,8 +613,8 @@ crear_log_p <- function(x, y){ mu_2 <- params[2] sigma_1 <- params[3] sigma_2 <- params[4] - log_vero_1 <- dnorm(x, mean = mu_1, sd = sigma_1, log = TRUE) %>% sum - log_vero_2 <- dnorm(y, mean = mu_2, sd = sigma_2, log = TRUE) %>% sum + log_vero_1 <- dnorm(x, mean = mu_1, sd = sigma_1, log = TRUE) |> sum() + log_vero_2 <- dnorm(y, mean = mu_2, sd = sigma_2, log = TRUE) |> sum() log_vero <- log_vero_1 + log_vero_2 #se suman por independiencia log_vero } @@ -622,8 +629,8 @@ crear_log_p_nula <- function(x, y){ mu <- params[1] sigma_1 <- params[2] sigma_2 <- params[3] - log_vero_1 <- dnorm(x, mean = mu, sd = sigma_1, log = TRUE) %>% sum - log_vero_2 <- dnorm(y, mean = mu, sd = sigma_2, log = TRUE) %>% sum + log_vero_1 <- dnorm(x, mean = mu, sd = sigma_1, log = TRUE) |> sum() + log_vero_2 <- dnorm(y, mean = mu, sd = sigma_2, log = TRUE) |> sum() log_vero <- log_vero_1 + log_vero_2 #se suman por independiencia log_vero } @@ -702,7 +709,7 @@ Y graficamos la distribución de referencia junto con el valor de $\lambda$ que obtuvimos: ```{r} -tibble(lambda = lambda_sim) %>% +tibble(lambda = lambda_sim) |> ggplot(aes(x = lambda)) + geom_histogram() + geom_vline(xintercept = lambda, colour = "red") @@ -718,7 +725,7 @@ una vez más que la distribución de referencia es cercana a una $\chi$-cuadrada con un grado de libertad. ```{r} -tibble(lambda = lambda_sim) %>% +tibble(lambda = lambda_sim) |> ggplot(aes(sample = lambda)) + geom_qq(distribution = stats::qchisq, dparams = list(df = 1)) + geom_qq_line(distribution = stats::qchisq, dparams = list(df = 1)) @@ -926,8 +933,8 @@ donde $\Phi$ es la función acumulada de la normal estándar. La gráfica de la función potencia es entonces ```{r} -potencia_tbl <- tibble(mu = seq(450, 550, 1)) %>% - mutate(beta = pnorm((505 - mu)/13)) %>% # probabilidad de rechazar +potencia_tbl <- tibble(mu = seq(450, 550, 0.5)) |> + mutate(beta = pnorm((505 - mu)/13)) |> # probabilidad de rechazar mutate(nula_verdadera = factor(mu >= 515)) # nula verdadera ggplot(potencia_tbl, aes(x = mu, y = beta, colour = nula_verdadera)) + geom_line() diff --git a/14-intro-bayesiana.Rmd b/14-intro-bayesiana.Rmd index 98efc0a..1511e39 100644 --- a/14-intro-bayesiana.Rmd +++ b/14-intro-bayesiana.Rmd @@ -101,10 +101,10 @@ Este es un ejemplo completo, aunque muy simple, de inferencia bayesiana. La estrategia de **inferencia bayesiana implica tomar decisiones basadas en las probabilidades posteriores.** -```{block, type='ejercicio'} -¿Cuál sería la estimación de máxima verosimilitud para este problema? ¿Cómo -cuantificaríamos la incertidumbre en la estimación de máxima verosimilitud? -``` + + + + Finalmente, podríamos hacer predicciones usando la **posterior predictiva**. Si ${X}_{nv}$ es una nueva tirada adicional de la moneda que estamos usando, nos