TP3 - Variables discretas

En este práctico vamos a trabajar con variables aleatorias discretas. Vamos a realizar experimentos simulados usando las capacidades de R (don’t panic ya que les daremos las líneas para eso), vamos a construir tibbles que almacenen la información generada, y luego vamos a usar los verbos de tidyverse para analizar los datos. A partir de esos datos, vamos a graficar histogramas y otras cosas útiles usando ggplot. Con todo eso, vamos a discutir algunos conceptos fundamentales de probabilidad y estadística.

Así que… manos a la obra!!

Carga de bibliotecas

library(magrittr)
library(purrr)
library(tidyverse)

Ejercicio 1

Las siguientes líneas de código generan N realizaciones de una variable aleatoria binaria (por ejemplo, una moneda) cuyos valores posibles son 0 y 1:

N = 50
set.seed(101)
realizaciones = base::sample(x=0:1, size=N, replace=T, prob=c(0.5, 0.5))

realizaciones
#>  [1] 1 1 0 0 1 1 0 1 0 0 0 0 0 0 1 0 0 1 1 1 0 0 1 0 0 0 1 1 1 0 1 1 1 1 0 0 1 0
#> [39] 1 0 1 1 0 1 1 1 1 0 0 1

Vamos a almacenar los datos de este “experimento” en un tibble:

datos = tibble(id=1:N, valor=realizaciones)

Armamos el tibble con dos columnas (variables del dataset). La primera (id) sirve para identificar cada realización. La segunda (valor) contiene las realizaciones.

datos
#> # A tibble: 50 x 2
#>       id valor
#>    <int> <int>
#>  1     1     1
#>  2     2     1
#>  3     3     0
#>  4     4     0
#>  5     5     1
#>  6     6     1
#>  7     7     0
#>  8     8     1
#>  9     9     0
#> 10    10     0
#> # ... with 40 more rows

  1. Usar mutate para agregar una nueva columna (piense un nombre apropiado, evite usar valor2 por ejemplo) en la que los valores de la variable se transformen de modo que un 0 corresponda a la palabra “ceca” y un 1 a “cara”.

    Ayuda: La función case_when() puede ser muy útil dentro del mutate.

  2. Obtener el número de apariciones de cada valor posible en el experimento (es decir, la frecuencia absoluta), tomando los valores de la columna que más le guste, y guardar el resultado en un nuevo tibble.

    Ayuda: Puede ser muy útil la función count().

    Luego, agregar una nueva columna que contenga la frecuencia relativa observada. (Notar que la cantidad de realizaciones está almacenada en la variable de entorno N)

    Nota: No confundir variables de un dataset (es decir, columnas de un tibble) con variables del entorno R (por ejemplo, la cantidad de datos N, u otros parámetros que nos permitan generar el dataset).

  3. A partir del tibble obtenido en el ítem anterior, realizar un gráfico de barras que muestre la frecuencia asociada a cada valor (es decir, un histograma). ¿Qué significa el resultado obtenido? ¿Son equiprobables los valores posibles?

  4. Repetir los pasos anteriores para un nuevo experimento que tiene muchas más realizaciones que el anterior.

    N = 5001
set.seed(101)
realizaciones = base::sample(x=0:1, size=N, replace=T, prob=c(0.5, 0.5))
datos = tibble(id = 1:N, valor=realizaciones)

    ¿Cambian sus conclusiones respecto a la pregunta del ítem anterior?

  5. A partir del último tibble generado (el grande), obtener una nueva columna (usando mutate) que contenga la suma parcial sobre las filas de resultados igual a alguno de los dos valores posibles (0 o 1). Por ejemplo, si las primeras filas corresponden a valores 0, 1, 0, 1, 1, 0, las primeras filas de la nueva columna deben valer 0, 1, 1, 2, 3, 3.

    Ayuda: Investigar las funciones sum() y cumsum().

  6. Realizar un gráfico que muestre la columna obtenida en el ítem anterior (en el eje y) versus el número de fila (en el eje x). Recordar que el número de fila está contenido en la columna id. ¿Es posible sacar alguna conclusión respecto a la probabilidad de aparición de cada valor a partir de este gráfico?

  7. Para pensar: Obtener una columna que contenga la frecuencia relativa parcial, es decir, la frecuencia relativa calculada a medida que se consideran nuevas realizaciones de la variable. Luego grafícar y discutir la pregunta de siempre. ¿A partir de qué momento cree que puede concluir con certeza cuál es la probabilidad asociada a cara y ceca? (Pensar y explicitar el criterio de certeza!!).

    Nota: Puede ser útil modificar el gráfico filtrando una cierta cantidad de datos al principio. ¿Cambian sus conclusiones al hacer esto?

Ejercicio 2

Considerar un experimento similar al anterior pero usando una moneda que no está balanceada (sus dos valores posibles tienen diferentes chances):

N = 5000
set.seed(101)
realizaciones = base::sample(x=0:1, size=N, replace=T, prob=c(0.2, 0.8))
datos2 = tibble(id = 1:N, valor=realizaciones)

Usando las herramientas discutidas en el ejercicio anterior, estimar la probabilidad de cada valor y cuántas realizaciones son necesarias para establecerla con “certeza”.

Ejercicio 3

Ahora vamos a hacer un experimento en el cuál se tiran 5 monedas justas (equiprobables) en 10 oportunidades.

Lo hago en dos pasos. Primero defino un tibble que identifica cada tirada y cada moneda mediante etiquetas numéricas:

# armo la estructura del experimento
N_monedas = 5
N_tiradas = 10

datos = list(id_tirada = 1:N_tiradas, id_moneda = 1:N_monedas) %>%
  cross_df() %>%
  arrange(id_tirada)

datos
#> # A tibble: 50 x 2
#>    id_tirada id_moneda
#>        <int>     <int>
#>  1         1         1
#>  2         1         2
#>  3         1         3
#>  4         1         4
#>  5         1         5
#>  6         2         1
#>  7         2         2
#>  8         2         3
#>  9         2         4
#> 10         2         5
#> # ... with 40 more rows

Luego “lleno” el tibble con los datos propiamente dichos (es decir las realizaciones), usando mutate y sample:

# realizo el experimento
set.seed(101)
datos %<>% mutate(valor = base::sample(x=0:1, size=N_monedas*N_tiradas, replace=T, prob=c(0.5, 0.5)))

datos
#> # A tibble: 50 x 3
#>    id_tirada id_moneda valor
#>        <int>     <int> <int>
#>  1         1         1     1
#>  2         1         2     1
#>  3         1         3     0
#>  4         1         4     0
#>  5         1         5     1
#>  6         2         1     1
#>  7         2         2     0
#>  8         2         3     1
#>  9         2         4     0
#> 10         2         5     0
#> # ... with 40 more rows

  1. Obtener, usando summarise y group_by, un nuevo tibble que contenga la cantidad de caras obtenidas en cada tirada de 5 monedas (llame X a dicha variable).

  2. A partir de ese tibble, armar uno nuevo que contenga la frecuencia de aparición de cada resultado posible. Recuerde que, si se tiran 5 monedas, la cantidad de caras iguales que puede obtener en cada tirada es un número entre 0 y 5 (incluyendo ambos valores, 0 y 5).

    Luego, verificar que la suma de frecuencias absolutas corresponde a la cantidad de tiradas, y obtener la frecuencia relativa para cada valor posible.

  3. Grafique el histograma para ambas frecuencias.

    Para pensar: ¿El gráfico contiene todos los valores posibles de la variable X? ¿Cómo puede incluirlos? ¿Cuál sería la frecuencia observada para dichos valores faltantes?

  4. El resultado de este tipo de experimentos se describe mediante la distribución binomial. Vamos a comparar el resultado del experimento con la distribución. Para esto, vamos a usar una función en R que nos da la frecuencia relativa esperada, en función de la cantidad de monedas tiradas en cada oportunidad, y de la probabilidad de cara y ceca:

    frec_relat = dbinom(0:N_monedas, size=N_monedas, prob=0.5)

frec_relat
#> [1] 0.03125 0.15625 0.31250 0.31250 0.15625 0.03125

    Vamos a construir un tibble que contenga los resultados esperados para nuestro experimento:

    binom = tibble(X = 0:N_monedas, f_teo = frec_relat)

binom
#> # A tibble: 6 x 2
#>       X  f_teo
#>   <int>  <dbl>
#> 1     0 0.0312
#> 2     1 0.156 
#> 3     2 0.312 
#> 4     3 0.312 
#> 5     4 0.156 
#> 6     5 0.0312

    A partir de este tibble, grafique la distribución.

  5. Compare la distribución con el histograma obtenido a partir del experimento. Para esto, obtenga un tibble que contenga las frecuencias relativas obtenidas en nuestro experimento y también las frecuencias teóricas según la distribución binomial (piense diferentes maneras de hacer esto).

    ¿Qué podríamos modificar para aumentar la similitud entre ambos? Verifícar computacionalmente.

Ejercicio 4 (optativo)

Este ejercicio es optativo. Pueden resolverlo en cualquier momento. Nuestra recomendación es intentarlo ya que repasa varias técnicas que serán útiles más adelante.

El objetivo principal es introducir varias funciones de tidyverse que permiten manipular varios datasets a la vez. En los primeros dos ejercicios, generamos un dataset a partir del cuál obteníamos cierta información. ¿Pero qué pasa si necesitamos obtener varias realizaciones de un mismo dataset, para estudiar, por ejemplo, la distribución de parámetros obtenidos? En el ejercicio 3 trabajamos de esa manera, generando un tibble que contenia varias repeticiones de un mismo experimento. En este ejercicio vamos a trabajar con varios verbos que nos permiten facilitar ese tipo de análisis, y que son imprescindibles para casos más complejos.

Comencemos por definir nuestro problema. Queremos estudiar una moneda con una cierta probabilidad de salir cara igual a p y queremos realizar un experimento donde la moneda se tira N veces. A su vez, este experimento será repetido N_rep veces. Definimos nuestros parámetros:

p = 0.5
N = 10
N_rep = 10

Vamos a guardar esos parámetros en una estructura cómoda para trabajar: una lista.

par <- list(p = p, N = N, N_rep = N_rep)

par
#> $p
#> [1] 0.5
#> 
#> $N
#> [1] 10
#> 
#> $N_rep
#> [1] 10

Usar una lista nos permite reducir la cantidad de variables de entorno que produce nuestro código. A su vez, cada componente de la lista tiene un nombre, lo cuál facilita la legibilidad del código (esto no significa que no sea necesario incluir una descripción de los mismos). De hecho, podríamos omitir la definición de las variables p, N y N_rep, y definir directamente la lista de parámetros. El resultado final es el mismo:

par <- list(p = 0.5,  # probabilidad de cara
            N = 10,   # cantidad de monedas tiradas en cada experimento
        N_rep = 10)   # cantidad de repeticiones del experimento

par
#> $p
#> [1] 0.5
#> 
#> $N
#> [1] 10
#> 
#> $N_rep
#> [1] 10

  1. Antes de simular nuestros datos, vamos a transformar nuestros parámetros (contenidos en la lista par) en un tibble mediante la función cross_df():

    sims = par %>% 
  cross_df() %>%
  mutate(id = 1:n(), .before = 1)

sims
#> # A tibble: 1 x 4
#>      id     p     N N_rep
#>   <int> <dbl> <dbl> <dbl>
#> 1     1   0.5    10    10

    Este tibble contiene una sola fila. ¿Qué ocurre si nuestra lista contiene más de un valor de cada parámetro? Por ejemplo, si incluímos dos valores posibles para la probabilidad:

    par2 <- list(p = c(0.1, 0.5), # dos monedas distintas
             N = 10,
         N_rep = 10)

par2
#> $p
#> [1] 0.1 0.5
#> 
#> $N
#> [1] 10
#> 
#> $N_rep
#> [1] 10

    Verificar el resultado que genera cross_df() para este conjunto de parámetros:

    sims2 <- par2 %>%
  cross_df() %>%
  mutate(id = 1:n(), .before = 1)

    Probar con otros valores posibles de los mismos e interpretar el resultado.

  2. Para realizar la simulación, construya un pipe que permita evaluar la función base::sample() según los parámetros contenidos en sims. ¿Cuál de los cuatro pipes de magrittr es el adecuado en este caso?

    sims %??% # escribir el pipe correcto
  base::sample(x    = 0:1,
               size = N*N_rep,
               prob = c(1-p, p),
            replace = T)

    Luego, emplee dicho pipe para simular los datos correspondientes a alguna de las combinaciones de parámetros correspondientes a lista par2. Para esto, investigue el verbo slice().

    par2 %>% 
  cross_df() %>%
  slice(??) %??% # escribir correctamente
  base::sample(x    = 0:1,
               size = N*N_rep,
               prob = c(1-p, p),
            replace = T)

  3. El siguiente bloque de código permite almacenar los datos que generamos en forma de un tibble anidado (nested):

    sims %<>%
  group_by(id) %>% # notar este pipe, y el uso de base::sample dentro de mutate
  mutate(data = base::sample(x    = 0:1,
                             size = N*N_rep,
                             prob = c(1-p, p),
                          replace = T) %>%  # <-- notar el uso de pipes
                                            #     dentro de mutate
           # guardo los datos en una columna llamada values.
           # el punto refiere a los datos pasados a través
           # del pipe previo
           tibble(values = .) %>%
           # el resultado final es una lista de tibbles
           list()
           )

sims
#> # A tibble: 1 x 5
#> # Groups:   id [1]
#>      id     p     N N_rep data              
#>   <int> <dbl> <dbl> <dbl> <list>            
#> 1     1   0.5    10    10 <tibble [100 x 1]>

    Verifique el resultado accediendo al contenido de la columna data:

    sims$data

    Verifique que el contenido de dicha columna corresponde a una lista. Luego, analice el contenido del primer elemento de dicha lista:

    sims$data[[1]]

    Compruebe que dicho contenido corresponde a un tibble de una única columna. Por último, repita el procedimiento usando los parámetros contenidos en par2 (o en sims2). ¿Cuántos tibbles resultan almacenados en la columna data?

    Nota: Antes de continuar, chequee que la cantidad de elementos dentro de cada elemento de sims$data corresponda a la esperada. Tenga esto en cuenta para el resto del ejercicio, para no exceder la capacidad de memoria de su computadora.

  4. La utilidad de los tibbles anidados es que nos permiten modificarlos y extraer información de ellos simultáneamente. De este modo, podemos tener diferentes versiones de un dataset (en nuestro caso, las correspondientes a diferentes parámetros que definen el experimento, pero también podría ser una versión del dataset con outliers y otra sin, o bien diferentes versiones correspondientes a diferentes métodos para extraer outliers) y obtener un mismo análisis a partir de cada una ellas, para después comparar sus resultados.

    Vamos a comenzar modificando los tibbles anidados. Para esto, vamos a incorporar la estructura de las tiradas, la cuál estaba definida por los parámetros N (cantidad de monedas) y N_rep (cantidad de repeticiones):

    sims %<>%
mutate(data = pmap(list(N, N_rep, data),
                   ~ data[[1]] %>%
                     mutate(id      = 1:(N*N_rep),               # id del dato, es el nro de fila
                            # rep permite repetir secuencias de dos maneras diferentes
                            # usando each o times según el caso
                            id_rep  = rep(1:N_rep, each = N),    # id de la repeticion, va de 1 a N_rep
                            id_coin = rep(1:N, times = N_rep),   # id de la moneda dentro de una repetición
                            .before = 1)
                  ))

sims$data[[1]]
#> # A tibble: 100 x 4
#>       id id_rep id_coin values
#>    <int>  <int>   <int>  <int>
#>  1     1      1       1      1
#>  2     2      1       2      1
#>  3     3      1       3      0
#>  4     4      1       4      0
#>  5     5      1       5      1
#>  6     6      1       6      0
#>  7     7      1       7      0
#>  8     8      1       8      0
#>  9     9      1       9      0
#> 10    10      1      10      0
#> # ... with 90 more rows

    Investigue la familia de funciones map del paquete purrr y analice (a partir del resultado del bloque de código anterior) cómo se mapean los parámetros N y N_rep a cada elemento dentro de sims$data.

  5. Piense de qué manera podría haber incorporado la estructura de las tiradas a los tibbles anidados en el momento de generar los datos. A continuación va una pista:

    sims %>%
  mutate(data = tibble(id      = ??,
                       id_rep  = ??,
                       id_coin = ??,
                       values  = base::sample(x    = 0:1,
                                              size = N*N_rep,
                                              prob = c(1-p, p),
                                           replace = T)
                       ) %>%
           list())

    Compare con el bloque de códigos provisto en el ítem (c).

  6. El siguiente bloque de código permite extraer la cantidad total de caras obtenidas en cada tibble anidado dentro de la columna data:

    sims %<>%
  mutate(suma = pmap(list(data),
                     ~ data[[1]] %$%
                         # sumo sobre la columna values de cada tibble
                         sum(values)) %>%
                  unlist())

sims
#> # A tibble: 1 x 6
#> # Groups:   id [1]
#>      id     p     N N_rep data                suma
#>   <int> <dbl> <dbl> <dbl> <list>             <int>
#> 1     1   0.5    10    10 <tibble [100 x 4]>    47

    Modifique el bloque de código para obtener el promedio de caras obtenido dentro de cada tibble. Verifique los resultados obtenidos mediante un pipe apropiado a partir de cada tibble individual (Ayuda: recuerde la sintaxis sims$data[[j]], donde j es un índice para acceder al elemento deseado dentro de data).

    Nota: Por simplicidad, en este ítem no estamos considerando al estructura de tiradas y repeticiones dentro de cada tibble.

  7. Interprete el resultado generado por el siguiente bloque de código:

    sims %<>%
  mutate(distrib_p = pmap(list(N, data),
                          ~ data[[1]] %>%
                              group_by(id_rep) %>%
                              summarise(p = sum(values) / N)))

sims
#> # A tibble: 1 x 7
#> # Groups:   id [1]
#>      id     p     N N_rep data                suma distrib_p        
#>   <int> <dbl> <dbl> <dbl> <list>             <int> <list>           
#> 1     1   0.5    10    10 <tibble [100 x 4]>    47 <tibble [10 x 2]>

    Luego, analice la lista de parámetros pasada a la función pmap(): list(N, data) y discuta si son todos necesarios.

    Ayuda: Evalúe cada elemento contenido en la columna distrib_p mediante la sintaxis sims$distrib_p[[j]], donde j es un índice para acceder al elemento j-ésimo.

  8. Interprete el resultado generado por el siguiente bloque de código:

    sims %<>%
  mutate(desvio = pmap(list(distrib_p),
                       ~ distrib_p[[1]] %$%
                           sd(p)) %>%
                    unlist())

sims
#> # A tibble: 1 x 8
#> # Groups:   id [1]
#>      id     p     N N_rep data                suma distrib_p         desvio
#>   <int> <dbl> <dbl> <dbl> <list>             <int> <list>             <dbl>
#> 1     1   0.5    10    10 <tibble [100 x 4]>    47 <tibble [10 x 2]>  0.157

  9. Para finalizar, construya un código que le permita estimar el desvío estándar asociado a la estimación de p en función de N, para N_rep fijo. Ayuda: Para empezar, puede considerar N entre 10 y 100 en pasos de 10, p igual a 0.5 y N_rep igual a 500.

  10. Bonus track: Analice el rol de N_rep en la curva obtenida en el ítem anterior.

    Sugerencia: Construya una simulación (a partir de un conjunto de parámetros apropiado) que le permita realizar dicho análisis reutilizando todo lo visto en el ejercicio.

    Nota: Para comenzar, analice dos valores de N_rep, por ejemplo, 50 y 500 (o 50 y 5000). Recuerde verificar que no se exceda la capacidad de memoria de su computadora.