Skip to content

Commit

Permalink
Deploying to gh-pages from @ 5c31f8b 🚀
Browse files Browse the repository at this point in the history
  • Loading branch information
@tereom
tereom committed Nov 28, 2024
1 parent 17b1939 commit 64f3bc6
Show file tree
Hide file tree
Showing 91 changed files with 8,811 additions and 1,062 deletions.
294 changes: 147 additions & 147 deletions 01-exploratorio.md
View file

Large diffs are not rendered by default.

32 changes: 16 additions & 16 deletions 10-bootstrap-parametrico.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ res$convergence
```

``` r
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 @@ -116,7 +116,7 @@ res_boot$convergence
```

``` r
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 @@ -176,7 +176,7 @@ 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
```
Expand All @@ -192,8 +192,8 @@ 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)
```

Expand Down Expand Up @@ -252,7 +252,7 @@ res$convergence
```

``` r
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 @@ -314,7 +314,7 @@ 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 @@ -459,11 +459,11 @@ 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 @@ -484,11 +484,11 @@ 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 @@ -630,9 +630,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))) %>%
bind_rows
reps_boot %>% mutate(across(everything(), round, 2)) %>% head()
n = length(muestra))) |>
bind_rows()
reps_boot |> mutate(across(everything(), round, 2)) |> head()
```

```
Expand Down
79 changes: 47 additions & 32 deletions 12-mas-hipotesis.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,11 +77,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 @@ -98,15 +98,16 @@ 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
```

```
## [1] 0.8410806
## [1] 0.8413082
```
y vemos que esta muestra es consistente con la media nacional. No tenemos
evidencia en contra de que la media del estado de México es muy similar
Expand Down Expand Up @@ -347,7 +348,7 @@ Supongamos por ejemplo que los datos que obtenemos son:


``` r
datos_clasif %>% head
datos_clasif |> head()
```

```
Expand All @@ -366,9 +367,9 @@ 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()
```

```
Expand All @@ -382,12 +383,26 @@ datos_clasif %>% head
## 5 5 0 1 -1
## 6 6 1 0 1
```

``` r
datos_clasif |>
summarise(sd_x = sd(x),
sd_y = sd(y),
sd_d = sd(d))
```

```
## # A tibble: 1 × 3
## sd_x sd_y sd_d
## <dbl> <dbl> <dbl>
## 1 0.393 0.309 0.539
```
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
```
Expand All @@ -402,10 +417,10 @@ 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 @@ -420,7 +435,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 @@ -483,7 +498,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 @@ -540,7 +555,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 @@ -584,7 +599,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 @@ -611,7 +626,7 @@ bajo la nula
``` r
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 @@ -636,9 +651,9 @@ 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 @@ -650,7 +665,7 @@ 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 @@ -696,8 +711,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 @@ -713,8 +728,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 @@ -835,7 +850,7 @@ 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 @@ -859,7 +874,7 @@ 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 @@ -1068,8 +1083,8 @@ 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
Binary file modified 12-mas-hipotesis_files/figure-html/unnamed-chunk-44-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading

0 comments on commit 64f3bc6

Please sign in to comment.