Cuaderno de trabajo de probabilidades

TI<-proc.time()

Un dado de \(k\) caras

En R podemos simular un dado de \(k\) caras fácilmente:

rDado<-function(n,k){
  return(ceiling(k*runif(n)))
}

Por ejemplo, si queremos observar el lanzamiento de 10 dados de 20 caras, podemos invocar:

replicate(10,rDado(1,20))

 [1]  4  1 13 20 16  4 18  7 20 14

o podríamos aprovechar la manera como escribimos la función para obtener un resultado similar…

rDado(10,20)

 [1]  1 12 14 16 15 20  9 17 15 18

En los diversos experimentos que realizaremos hoy, vamos a trabajar con dados de seis caras.

Podemos escribir una función para simular un dado de seis caras, que sirva para encapsular la función anterior y nos permita trabajar con los diversos experimentos planteados.

d6<-function(n){
  return(rDado(n,6))
}

Ahora podemos utilizar la función para hacer varios lanzamientos de dados de seis, por ejemplo siete dados de seis

d6(7)

[1] 4 5 2 2 1 1 1

Ahora vamos a atacar desde el punto de vista empírico los experimentos mencionados en clase.

Experimentos

La idea es repetir el experimento suficientes veces hasta tener una idea razonable de cuanto valen las diversas probabilidades y cual es el espacio de resultados.

Es conveniente determinar una buena forma de visualizar cada corrida de simulación para tener una idea de los diversos valores.

Evidentemente la visualización correcta para cada caso depende de una multiplicidad de factores, como son:

• El tipo de experimento y su espacio de resultados.
• La cantidad de veces que se repite el experimento.

E1: Un dado de seis

Este experimento es relativamente sencillo, el espacio de resultados consiste en el conjunto de los valores \(\{1,2,3,4,5,6\}\) y las probabilidades de cada valor son todas iguales a \(\frac{1}{6}\) ya que el dado es justo.

Para mostrar esto podemos hacer el experimento varias veces.

Pocas veces

E11<-d6(60)
table(E11)/length(E11)

E11
        1         2         3         4         5         6 
0.1000000 0.2333333 0.1666667 0.2166667 0.1833333 0.1000000 

summary(E11)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    2.00    3.50    3.45    5.00    6.00 

stripchart(E11,method="stack",pch=19)
abline(h=1-0.05)

points(mean(E11),1-0.1,pch=17,col="red")

hist(E11,col="lightgreen",breaks=1:7-0.5,prob=TRUE)

boxplot(E11,col="lightgreen")

Muchas veces

E12<-d6(60000)
table(E12)/length(E12)

E12
        1         2         3         4         5         6 
0.1665667 0.1660167 0.1642667 0.1670667 0.1681667 0.1679167 

summary(E12)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   2.000   4.000   3.508   5.000   6.000 

hist(E12,col="lightgreen",breaks=1:7-0.5,prob=TRUE)

boxplot(E12,col="lightgreen")

E2: Suma de tres dados de seis

rSumaDe3<-function(n){
  replicate(n,sum(d6(3)))
}

Este experimento es un poco más complicado, el espacio de resultados consiste en el conjunto de los valores \(\{3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18\}\) y las probabilidades de cada valor se pueden calcular trabajando sobre el espacio equiprobable de las tripletas \((x_1,x_2,x_3)\) con \(x_i\in\{1,2,3,4,5,6\}\) para \(i\in\{1,2,3\}\).

Para estimar las probabilidades repetir el experimento varias veces.

Pocas veces

E21<-rSumaDe3(216)
table(E21)/length(E21)

E21
          3           4           5           6           7 
0.004629630 0.009259259 0.027777778 0.055555556 0.060185185 
          8           9          10          11          12 
0.097222222 0.148148148 0.101851852 0.087962963 0.138888889 
         13          14          15          16          17 
0.101851852 0.069444444 0.041666667 0.037037037 0.004629630 
         18 
0.013888889 

summary(E21)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.00    8.00   10.00   10.55   13.00   18.00 

stripchart(E21,method="stack",pch=19,ylim=c(0.8,2),xlim=c(3,18))
abline(h=1-0.05)

points(mean(E21),1-0.1,pch=17,col="red")

hist(E21,col="lightgreen",breaks=3:19-0.5,prob=TRUE)

boxplot(E21,col="lightgreen")

Como se puede apreciar, ya el diagrama de puntos no es particularmente útil, pues en realidad nos da casi la misma información que el histograma. ### Muchas veces

E22<-rSumaDe3(60000)
table(E22)/length(E22)

E22
          3           4           5           6           7 
0.004250000 0.013966667 0.028566667 0.045333333 0.068816667 
          8           9          10          11          12 
0.097333333 0.113816667 0.125983333 0.126266667 0.115866667 
         13          14          15          16          17 
0.097016667 0.069133333 0.047916667 0.027783333 0.013716667 
         18 
0.004233333 

summary(E22)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3.00    8.00   11.00   10.51   13.00   18.00 

hist(E22,col="lightgreen",breaks=3:19-0.5,prob=TRUE)

boxplot(E22,col="lightgreen")

E3: Número de veces que sale el seis en veinte lanzamientos

De nuevo, el experimento es más complicado que el anterior. Debería estar claro que el espacio de resultados es el conjunto de los números \(\{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20\}\) pero ahora deberíamos notar que incluso una cantidad enorme de réplicas del experimento, no revela el espacio de resultados completo, por lo que para estimar de manera razonable las probabilidades de manera empírica, tenemos que irnos a números verdaderamente grandes.

rTotalDe6En20<-function(n){
  replicate(n,sum(d6(20)==6))
}

m<-600000
ti<-proc.time()
E32<-rTotalDe6En20(m)
tf<-proc.time()
t<-(tf-ti)[3]
tiempolargo<-((6^20/m)*t)/(3600*24*365.25)
table(E32)/length(E32)

E32
           0            1            2            3            4 
2.587667e-02 1.039450e-01 1.987567e-01 2.384750e-01 2.018067e-01 
           5            6            7            8            9 
1.296633e-01 6.480167e-02 2.557833e-02 8.296667e-03 2.258333e-03 
          10           11           12           13 
4.583333e-04 7.000000e-05 1.166667e-05 1.666667e-06 

summary(E32)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.000   3.000   3.332   4.000  13.000 

hist(E32,col="lightgreen",breaks=0:21-0.5,prob=TRUE)

boxplot(E32,col="lightgreen")

De hecho, para esperar observar una vez el resultado \(20\), debemos repetir el experimento \(6^{20}=\) 3,656,158,440,062,976 veces. Esto ya resulta impráctico ya que el resultado observado se tardó 34.05 segundos, por lo que la simulación larga se tardaría aproximadamente 6574.9 años.

E4: Número de lanzamientos hasta que sale seis por primera vez

El experimento ahora es genuinamente complicado (bueno, no realmente) y requiere entender el uso de la instrución while. Como R es un lenguaje de programación podemos utilizarlo para simular el experimento 4 de manera directa.

El espacio de resultados ahora es \(\{1,2,3,\ldots\}=\mathbb{N}\). Como este espacio de resultados es infinito, la idea de la tabla de valores ya no es factible, por lo que obtener aproximaciones a las probabilidades de manera empírica solo es práctico para valores relativamente pequeños de el número de lanzamientos hasta observar el primer seis.

rNHP6<-function(){
  listo<-FALSE
  i<-0
  while(!listo){
    i<-i+1
    listo<-(d6(1)==6)
  }
  return(i)
}
rNumeroHastaPrimer6<-function(n){
  return(replicate(n,rNHP6()))
}

E42<-rNumeroHastaPrimer6(1000)
table(E42)/length(E42)

E42
    1     2     3     4     5     6     7     8     9    10    11 
0.173 0.147 0.119 0.101 0.069 0.075 0.047 0.039 0.043 0.033 0.024 
   12    13    14    15    16    17    18    19    20    21    22 
0.021 0.015 0.012 0.018 0.006 0.005 0.008 0.005 0.014 0.003 0.004 
   23    24    25    26    27    30    31    32    35 
0.003 0.003 0.005 0.002 0.001 0.002 0.001 0.001 0.001 

summary(E42)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   2.000   4.000   5.884   8.000  35.000 

hist(E42,col="lightgreen",breaks=0:(max(E42)+1)-0.5,prob=TRUE)

boxplot(E42,col="lightgreen")

E5: Número de lanzamientos hasta que salen dos seis consecutivos

Para simular el experimento 5, deben modificar la función que simula el experimento 4 de manera cuidadosa.

El espacio de resultados es \(\{2,3,\ldots\}=\mathbb{N}\setminus\{1\}\). Estimar las probabilidades por simulación es difícil si no se tiene una idea clara de cuán grande debe ser el número de réplicas.

Una dificultad adicional en este problema es que, a diferencia de los experimentos anteriores, saber si el próximo lanzamiento concluye el experimento es algo que depende del lanzamiento actual, si no fue seis sabemos que el experimento no puede terminar en el próximo lanzamiento, mientras que si el lanzamiento actual fue un seis entonces el experimento termina en el próximo lanzamiento con probabilidad \(\frac{1}{6}\).

E6: Número de lanzamientos hasta que sale el octavo seis

Para simular el experimento 6, deben modificar la función que simula el experimento 4 de manera cuidadosa.

E7: Número de valores distintos que se observan en 6 lanzamientos

E8: Número de valores distintos que se observan en 30 lanzamientos

E9: Número de lanzamientos hasta que se repite un valor de manera consecutiva

E10: Número de veces hasta que se observan todos los valores del dado

Proyecto de laboratorio

1. Reportar los diversos descriptores numéricos mencionados en clase para cada experimento. Puede ser interesante leer las ayudas de summary, quantile,median, mean, sd, var, min y max entre otras. Para el coeficiente de variación y el rango intercuartil deben escribir funciones cdv y ric que dada una muestra cualquiera, determinen sus valores de acuerdo a la definición. NO deben usar las funciones prefabricadas de R.

2. Realizar una visualización informativa donde se comparen los resultados teóricos y los resultados empíricos para cada experimento. Para poder hacer esto es fundamental entender el espacio de resultados de cada experimento y poder encontrar la probabilidad teórica de cada resultado.

3. Especificar, cuando sea posible, que variable aleatoria sirve para modelar cada uno de los experimentos.

4. Reportar los valores esperados teóricos de todos los experimentos.

La entrega debe ser un documento .Rmd que siga el formato de este documento (de hecho, podrían simplemente rellenar los espacios vacíos de este documento).

El tiempo de ejecución debe ser razonable, es decir, el número de réplicas debe ser suficientemente grande para que las visualizaciones sean informativas pero suficientemente pequeño para no tardarse una eternidad a la hora de revisarlo. Nota: La diversidad de visualizaciones y las subsecciones de Pocas veces y Muchas veces son solo para dar una idea de lo que se puede hacer y no deben considerarse como una propuesta de visualización informativa, en particular, ni los ejes, ni los títulos fueron etiquetados apropiadamente por tratarse este documento de una guía de laboratorio con la que deben practicar. Entre los otros documentos que ya han visto pueden ver ejemplos de como etiquetar apropiadamente los gráficos y como decidir que presentación es pertinente. En cualquier caso, recuerden que tienen un libro de texto. Advertencia: Los experimentos 9 y 10 son difíciles de resolver desde el punto de vista teórico.

TF<-proc.time()
tT<-(TF-TI)[3]

Este documento tardó 57.32 segundos en ser generado.