Skip to content

Commit

Permalink
Cambiar a nuevo pipe
Browse files Browse the repository at this point in the history
  • Loading branch information
@tereom
tereom committed Nov 20, 2024
1 parent 6ce5d8e commit d85f144
Show file tree
Hide file tree
Showing 3 changed files with 63 additions and 56 deletions.
30 changes: 15 additions & 15 deletions 10-bootstrap-parametrico.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")
```

Expand Down Expand Up @@ -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
```
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
```
Expand Down Expand Up @@ -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()
```
Expand Down Expand Up @@ -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
```
Expand All @@ -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")
Expand Down Expand Up @@ -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:
Expand Down
81 changes: 44 additions & 37 deletions 12-mas-hipotesis.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
```
Expand All @@ -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
Expand Down Expand Up @@ -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
```
Expand All @@ -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)
```


Expand Down Expand Up @@ -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() +
Expand Down Expand Up @@ -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) +
Expand Down Expand Up @@ -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
}
}
Expand All @@ -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")
```
Expand All @@ -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")
```
Expand All @@ -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))
Expand Down Expand Up @@ -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
}
Expand All @@ -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
}
Expand Down Expand Up @@ -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")
Expand All @@ -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))
Expand Down Expand Up @@ -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()
Expand Down
8 changes: 4 additions & 4 deletions 14-intro-bayesiana.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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?
```
<!-- ```{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
Expand Down

0 comments on commit d85f144

Please sign in to comment.