Preliminares.

La simulación es una de las herramientas computacionales más útiles y constituye una pieza fundamental del arsenal de métodos que se usan en Computación Científica. Simular los resultados de un experimento antes de llevarlo a cabo permite a menudo detectar problemas de diseño del experimento en fases preliminares de su desarrollo, evitando así pérdidas de recursos. Aquí vamos a empezar a ver los rudimentos de las técnicas de simulación que podemos implementar con R, aprovechando los ejercicios de probabilidad y variables aleatorias que estamos estudiando.

Matrices para las simulaciones.

Las matrices son un tipo especial de tablas de R (de tipo matrix), que tienen la propiedad de que todos los elementos de la tabla son del mismo tipo: todos números, todos factores, etc.

Por ejemplo, este código fabrica una matriz de con los números del 1 al 36:

(M = matrix(1:36, nrow=4) )
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    1    5    9   13   17   21   25   29   33
## [2,]    2    6   10   14   18   22   26   30   34
## [3,]    3    7   11   15   19   23   27   31   35
## [4,]    4    8   12   16   20   24   28   32   36

Fíjate en que R rellena la matriz por columnas. Para que lo haga al revés ponemos:

(M = matrix(1:36, nrow=4, byrow = TRUE) )
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    1    2    3    4    5    6    7    8    9
## [2,]   10   11   12   13   14   15   16   17   18
## [3,]   19   20   21   22   23   24   25   26   27
## [4,]   28   29   30   31   32   33   34   35   36

En una simulación lo habitual es repetir un mismo experimento un cierto número de veces (en general. muchas veces). Y por eso surge la necesidad de almacenar los resultados de cada una de esas repeticiones (iteraciones). Las matrices nos sirven muchas veces para esto, ya que podemos almacenar los resultados de cada repetición en una fila de la matriz, cuando todos esos resultados son del mismo tipo. Si los resultados no son del mismo tipo, debemos usar un data.frame.

La forma más habitual de fabricar una matriz es partiendo de un vector. Arriba hemos usado el vector 1:36, pero puedes usar cualquier vector, siempre que las dimensiones cuadren. Por ejemplo, si fabricamos un vector con 100 números al azar:

set.seed(2017)
(valores = sample(-30:30, size = 100, replace = TRUE))
##   [1]  26   2  -2 -13  16  17 -28  -4  -2 -14  11 -30 -29  -4   0  -7  -6
##  [18]  13  27  20   8  -7  27  16 -17  25   7 -11 -26 -15  23 -29 -11 -18
##  [35]  25   9 -23 -24 -20   2   3   6 -18 -26  29   1 -27  11  17  15   2
##  [52]  23  -9   8  15 -13 -19  29   2  26 -15  -8   6  15  14  -2 -20 -24
##  [69]   7  10  30  12 -21  23  29  -8  -7   1  -5  26 -22  22   1  26 -16
##  [86] -19 -21 -20  -8 -25  22   9   1   3  21 -23  20 -17  24 -20

los podemos usar para crear una matriz 10x10:

(A = matrix(valores, nrow = 10, byrow = TRUE))
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]   26    2   -2  -13   16   17  -28   -4   -2   -14
##  [2,]   11  -30  -29   -4    0   -7   -6   13   27    20
##  [3,]    8   -7   27   16  -17   25    7  -11  -26   -15
##  [4,]   23  -29  -11  -18   25    9  -23  -24  -20     2
##  [5,]    3    6  -18  -26   29    1  -27   11   17    15
##  [6,]    2   23   -9    8   15  -13  -19   29    2    26
##  [7,]  -15   -8    6   15   14   -2  -20  -24    7    10
##  [8,]   30   12  -21   23   29   -8   -7    1   -5    26
##  [9,]  -22   22    1   26  -16  -19  -21  -20   -8   -25
## [10,]   22    9    1    3   21  -23   20  -17   24   -20

o una matriz 4x25:

(A = matrix(valores, nrow = 25, byrow = TRUE))
##       [,1] [,2] [,3] [,4]
##  [1,]   26    2   -2  -13
##  [2,]   16   17  -28   -4
##  [3,]   -2  -14   11  -30
##  [4,]  -29   -4    0   -7
##  [5,]   -6   13   27   20
##  [6,]    8   -7   27   16
##  [7,]  -17   25    7  -11
##  [8,]  -26  -15   23  -29
##  [9,]  -11  -18   25    9
## [10,]  -23  -24  -20    2
## [11,]    3    6  -18  -26
## [12,]   29    1  -27   11
## [13,]   17   15    2   23
## [14,]   -9    8   15  -13
## [15,]  -19   29    2   26
## [16,]  -15   -8    6   15
## [17,]   14   -2  -20  -24
## [18,]    7   10   30   12
## [19,]  -21   23   29   -8
## [20,]   -7    1   -5   26
## [21,]  -22   22    1   26
## [22,]  -16  -19  -21  -20
## [23,]   -8  -25   22    9
## [24,]    1    3   21  -23
## [25,]   20  -17   24  -20

Comprueba que el comando

(A = matrix(valores, ncol = 4, byrow = TRUE))

produce el mismo resultado.

Si aplicas c a una matriz recuperas el vector, pero cuidado: ¡leído por columnas!

c(A)
##   [1]  26  16  -2 -29  -6   8 -17 -26 -11 -23   3  29  17  -9 -19 -15  14
##  [18]   7 -21  -7 -22 -16  -8   1  20   2  17 -14  -4  13  -7  25 -15 -18
##  [35] -24   6   1  15   8  29  -8  -2  10  23   1  22 -19 -25   3 -17  -2
##  [52] -28  11   0  27  27   7  23  25 -20 -18 -27   2  15   2   6 -20  30
##  [69]  29  -5   1 -21  22  21  24 -13  -4 -30  -7  20  16 -11 -29   9   2
##  [86] -26  11  23 -13  26  15 -24  12  -8  26  26 -20   9 -23 -20

La función t (de transpose) intercambia las filas y columnas de una matriz.

La función dim permite ver, pero también cambiar las dimensiones de la matriz.

M
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    1    2    3    4    5    6    7    8    9
## [2,]   10   11   12   13   14   15   16   17   18
## [3,]   19   20   21   22   23   24   25   26   27
## [4,]   28   29   30   31   32   33   34   35   36
t(M)
##       [,1] [,2] [,3] [,4]
##  [1,]    1   10   19   28
##  [2,]    2   11   20   29
##  [3,]    3   12   21   30
##  [4,]    4   13   22   31
##  [5,]    5   14   23   32
##  [6,]    6   15   24   33
##  [7,]    7   16   25   34
##  [8,]    8   17   26   35
##  [9,]    9   18   27   36
dim(M)
## [1] 4 9
dim(M) = c(3, 12)
M
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,]    1   28   20   12    4   31   23   15    7    34    26    18
## [2,]   10    2   29   21   13    5   32   24   16     8    35    27
## [3,]   19   11    3   30   22   14    6   33   25    17     9    36

Vamos a ver en ejemplos otras funciones matriciales muy útiles para las simulaciones:

(W = matrix(1:12, nrow = 3))
##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12
colSums(W)
## [1]  6 15 24 33
rowSums(W)
## [1] 22 26 30
colMeans(W)
## [1]  2  5  8 11
rowMeans(W)
## [1] 5.5 6.5 7.5

Las matrices, en el fondo, son vectores. Si aplicas una función vectorial a una matriz se aplica elemento a elemento.

W + 3
##      [,1] [,2] [,3] [,4]
## [1,]    4    7   10   13
## [2,]    5    8   11   14
## [3,]    6    9   12   15
W^2
##      [,1] [,2] [,3] [,4]
## [1,]    1   16   49  100
## [2,]    4   25   64  121
## [3,]    9   36   81  144
sqrt(W)
##          [,1]     [,2]     [,3]     [,4]
## [1,] 1.000000 2.000000 2.645751 3.162278
## [2,] 1.414214 2.236068 2.828427 3.316625
## [3,] 1.732051 2.449490 3.000000 3.464102

Recuerda también que podemos aplicar condiciones lógicas a las matrices y que se aplicarán elemento a elemento, devolviendo una matriz de valores TRUE/FALSE

W > 5
##       [,1]  [,2] [,3] [,4]
## [1,] FALSE FALSE TRUE TRUE
## [2,] FALSE FALSE TRUE TRUE
## [3,] FALSE  TRUE TRUE TRUE

Finalmente recuerda que los valores TRUE equivalen a 1, mientras que los FALSE equivalen a 0. Para saber cuántos TRUE hay en cada fila basta con usar rowSums.

rowSums(W > 5)
## [1] 2 2 3

Herramientas para las simulaciones.

La función replicate

Esta función facilita mucho el trabajo de construcción de simulaciones con R. La estructura básica es:

repeticiones = replicate(n = MuchasVeces,{
  # Aquí entre las dos llaves escribimos 
  # las instrucciones de R necesarias
  # para fabricar el resultado de cada
  # repetición del experimento.
})

Por ejemplo, vamos a usar replicate para simular el primer juego del Caballero de Mére. Recuerda que en cada partida tiramos un dado cuatro veces. Así que una partida se puede conseguir con:

sample(1:6, size = 4, replace = TRUE)
## [1] 6 6 1 2

Ahora usamos replicate para simular 7 partidas de ese juego:

numRepeticiones = 7
repeticiones = replicate(n = numRepeticiones,{
  partida = sample(1:6, size = 4, replace = TRUE)
})

¿Qué tipo de objeto hemos obtenido como resultado?

class(repeticiones)
## [1] "matrix"
dim(repeticiones)
## [1] 4 7
repeticiones
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    4    1    6    3    5    1    5
## [2,]    1    1    2    1    3    2    4
## [3,]    3    1    6    3    4    4    4
## [4,]    4    1    4    2    3    3    2

Como ves, replicate ha fabricado una matriz, colocando cada partida en una columna. Si prefieres ver cada partida en una fila basta con usar t para trasponer la matriz:

(repeticiones = t(repeticiones))
##      [,1] [,2] [,3] [,4]
## [1,]    4    1    3    4
## [2,]    1    1    1    1
## [3,]    6    2    6    4
## [4,]    3    1    3    2
## [5,]    5    3    4    3
## [6,]    1    2    4    3
## [7,]    5    4    4    2

Las partidas ganadoras son aquellas que incluyen un 6. Puedes localizarlas usando:

repeticiones == 6
##       [,1]  [,2]  [,3]  [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,]  TRUE FALSE  TRUE FALSE
## [4,] FALSE FALSE FALSE FALSE
## [5,] FALSE FALSE FALSE FALSE
## [6,] FALSE FALSE FALSE FALSE
## [7,] FALSE FALSE FALSE FALSE

Y sumando por filas para ver si hay algún 6 en cada una de las filas:

(filasCon6 = (rowSums(repeticiones == 6) > 0))
## [1] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE

En la Sección 3 del Tutorial 3 puedes ver otra forma de hacer esta simulación, sin usar replicate.

Ejercicio 1.

  • Usa estas ideas para simular \(100000\) partidas del primer juego del Caballero De Méré. Haz una tabla de frecuencias relativas de partidas ganadoras / perdedoras.
  • Calcula la probabilidad teórica de ganar en una partida de ese juego y compárala con los resultados de la simulación.

Ejercicio 2.

  • Vamos a hacer el Ejercicio 1 del Tutorial 4 sobre variables aleatorias. Lo que queremos es simular \(n = 1000000\) (un millón) tiradas de dos dados, y calcular la tabla de frecuencias relativas de la variable aleatoria \[X = \{\text{suma de los dos dados}\}\]

  • ¿Cuál era la tabla teórica de valores y probabilidades de esta variable?

  • ¿Cuál es la media de esas sumas (las simuladas)? ¿Y cuál era la media teórica?

La función sample usando probabilidades.

Hasta ahora hemos visto ejemplos de uso de sample como este:

(dado = sample(1:6, 100, replace = TRUE))
##   [1] 2 2 5 2 5 3 2 3 6 6 2 2 2 2 1 6 6 5 6 3 2 4 3 5 5 2 4 1 1 3 6 1 3 4 4
##  [36] 5 4 3 4 4 5 4 5 2 3 2 3 6 3 2 4 2 1 5 1 5 2 3 5 4 3 3 1 1 2 5 2 1 4 5
##  [71] 6 6 3 5 5 3 2 2 1 1 5 1 4 5 6 2 2 5 2 1 3 6 4 6 2 5 2 4 5 5

que nos permite simular 100 tiradas de un dado “honesto”. Es decir, que los seis posibles resultados son equiprobables. Pero a veces necesitaremos simular experimentos en los que el dado está cargado y no hay equiprobabilidad. Supongamos, por ejemplo, que queremos simular un dado en el que el 6 es 5 veces más probable que cualquiera de los otros resultados. Podemos usar sample, per añadimos un argumento adicional probs para indicar esas probabilidades así:

(dadoCargado = sample(1:6, 100, replace = TRUE, prob = c(1, 1, 1, 1, 1, 5)))
##   [1] 1 6 3 6 5 6 6 6 4 5 6 3 6 5 6 6 6 1 1 1 5 6 6 6 2 6 6 5 6 2 4 2 2 6 6
##  [36] 5 6 6 6 5 4 4 1 6 6 6 6 6 6 5 3 6 6 6 5 1 6 3 3 4 6 2 6 6 3 5 4 2 6 6
##  [71] 6 6 5 6 6 6 2 1 6 1 3 2 6 6 5 6 2 6 6 6 6 6 4 6 6 6 1 6 5 1

Como ves, él resultado contiene muchos más seises que en el caso de un dado honesto.

La instrucción if / else.

El ejemplo 3.5.1 del libro (pág. 67) dice lo siguiente:

Supongamos que tenemos dos urnas, la primera con 3 bolas blancas y dos negras, y la segunda con 4 bolas blancas y 1 negra. Para extraer una bola lanzamos un dado. Si el resultado es 1 o 2 usamos la primera urna; si es cualquier otro número usamos la segunda urna. ¿cuál es la probabilidad de obtener una bola blanca?

Vamos a aprender a simular este tipo de situaciones. Pero como ves, el proceso que conduce a saber si la bola es blanca o negra tiene dos etapas. Y la segunda etapa depende del resultado de la primera. Así que necesitamos una forma de decirle a R algo como “si pasa esto haz tal cosa y si no pasa haz tal otra”.

La instrucción if...else sirve precisamente para eso. Su estructura general es esta:


if(se cumple cierta condicion){
  vamos por el camino A
} else {
  vamos por el camino B
}

Por ejemplo, para decidir el color de la bola que extraemos en el ejemplo que nos ocupa podríamos usar if...else combinado con la versión de sample con probabilidades que acabamos de ver así:

(dado = sample(1:6, 1))
## [1] 3
bola = if(dado < 3){
  print("Usamos la urna 1")
  sample(c("b", "n"), 1, prob = c(3, 2)) # urna 1
} else {
  print("Usamos la urna 2")
  sample(c("b", "n"), 1, prob = c(4, 1)) # urna 1
}
## [1] "Usamos la urna 2"
bola
## [1] "b"

Ejecuta varias veces el código para convencerte de que funcione. Hemos usado la función print para que sea más fácil esa comprobación, aunque en la simulación final no la incluiremos.

Ejercicio 2.

  • Usando replicate y el fragmento de código que acabamos de ver, simula 10000 extracciones de la bola en el ejemplo anterior.
  • Haz una tabla de frecuencias relativas de bola blanca y negra.
  • Calcula la probabilidad teórica y compárala con la simulación.

Combinatoria con R.

Al realizar ejercicios de probabilidad nos hemos visto varias veces en la necesidad de pensar en la lista completa de permutaciones o combinaciones de una colección de objetos. Muchas veces nos basta con saber cuántos elementos tiene esa lista. Pero en algunas simulaciones vamos a querer construir la lista explícitamente. Vamos a ver brevemente algunos ejemplos sencillos.

Usaremos la librería gtools de R. Si no la tienes instalada el siguiente comando la instalará (siempre que tengas conexión, claro).

if(!require(gtools)) install.packages(gtools)
## Loading required package: gtools

Para construir una lista de permutaciones (es decir, que el orden importa) de 4 elementos tomados de 3 en 3 usaremos:

permutations(4,3)
##       [,1] [,2] [,3]
##  [1,]    1    2    3
##  [2,]    1    2    4
##  [3,]    1    3    2
##  [4,]    1    3    4
##  [5,]    1    4    2
##  [6,]    1    4    3
##  [7,]    2    1    3
##  [8,]    2    1    4
##  [9,]    2    3    1
## [10,]    2    3    4
## [11,]    2    4    1
## [12,]    2    4    3
## [13,]    3    1    2
## [14,]    3    1    4
## [15,]    3    2    1
## [16,]    3    2    4
## [17,]    3    4    1
## [18,]    3    4    2
## [19,]    4    1    2
## [20,]    4    1    3
## [21,]    4    2    1
## [22,]    4    2    3
## [23,]    4    3    1
## [24,]    4    3    2

Son 24 porque \(4\cdot 3\cdot 2 = 24\). Este tipo de factoriales incompletos se puede calcular en R con prod:

prod(4:2)
## [1] 24

Si lo que queremos es la lista de combinaciones (el orden no importa) de 4 elementos tomados de 3 en 3 usamos:

combinations(4, 3)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    2    4
## [3,]    1    3    4
## [4,]    2    3    4

Son muchas menos (la sexta parte), porque:

choose(4, 3)
## [1] 4

Un detalle importante. Si quieres usar otros vectores como base de las permutaciones y combinaciones puedes hacerlo con un argumento adicional. Por ejemplo, ¿cuántos codones hay que tengan todos sus nucleótidos distintos?

permutations(n = 4, r = 3, v = c("A", "T", "G", "C"))
##       [,1] [,2] [,3]
##  [1,] "A"  "C"  "G" 
##  [2,] "A"  "C"  "T" 
##  [3,] "A"  "G"  "C" 
##  [4,] "A"  "G"  "T" 
##  [5,] "A"  "T"  "C" 
##  [6,] "A"  "T"  "G" 
##  [7,] "C"  "A"  "G" 
##  [8,] "C"  "A"  "T" 
##  [9,] "C"  "G"  "A" 
## [10,] "C"  "G"  "T" 
## [11,] "C"  "T"  "A" 
## [12,] "C"  "T"  "G" 
## [13,] "G"  "A"  "C" 
## [14,] "G"  "A"  "T" 
## [15,] "G"  "C"  "A" 
## [16,] "G"  "C"  "T" 
## [17,] "G"  "T"  "A" 
## [18,] "G"  "T"  "C" 
## [19,] "T"  "A"  "C" 
## [20,] "T"  "A"  "G" 
## [21,] "T"  "C"  "A" 
## [22,] "T"  "C"  "G" 
## [23,] "T"  "G"  "A" 
## [24,] "T"  "G"  "C"

Este último ejemplo nos lleva a la pregunta natural ¿y si quiero permitir repeticiones para obtener la lista completa de codones posibles? De nuevo, basta con un argumento adicional:

permutations(n = 4, r = 3, v = c("A", "T", "G", "C"), repeats.allowed = TRUE)
##       [,1] [,2] [,3]
##  [1,] "A"  "A"  "A" 
##  [2,] "A"  "A"  "C" 
##  [3,] "A"  "A"  "G" 
##  [4,] "A"  "A"  "T" 
##  [5,] "A"  "C"  "A" 
##  [6,] "A"  "C"  "C" 
##  [7,] "A"  "C"  "G" 
##  [8,] "A"  "C"  "T" 
##  [9,] "A"  "G"  "A" 
## [10,] "A"  "G"  "C" 
## [11,] "A"  "G"  "G" 
## [12,] "A"  "G"  "T" 
## [13,] "A"  "T"  "A" 
## [14,] "A"  "T"  "C" 
## [15,] "A"  "T"  "G" 
## [16,] "A"  "T"  "T" 
## [17,] "C"  "A"  "A" 
## [18,] "C"  "A"  "C" 
## [19,] "C"  "A"  "G" 
## [20,] "C"  "A"  "T" 
## [21,] "C"  "C"  "A" 
## [22,] "C"  "C"  "C" 
## [23,] "C"  "C"  "G" 
## [24,] "C"  "C"  "T" 
## [25,] "C"  "G"  "A" 
## [26,] "C"  "G"  "C" 
## [27,] "C"  "G"  "G" 
## [28,] "C"  "G"  "T" 
## [29,] "C"  "T"  "A" 
## [30,] "C"  "T"  "C" 
## [31,] "C"  "T"  "G" 
## [32,] "C"  "T"  "T" 
## [33,] "G"  "A"  "A" 
## [34,] "G"  "A"  "C" 
## [35,] "G"  "A"  "G" 
## [36,] "G"  "A"  "T" 
## [37,] "G"  "C"  "A" 
## [38,] "G"  "C"  "C" 
## [39,] "G"  "C"  "G" 
## [40,] "G"  "C"  "T" 
## [41,] "G"  "G"  "A" 
## [42,] "G"  "G"  "C" 
## [43,] "G"  "G"  "G" 
## [44,] "G"  "G"  "T" 
## [45,] "G"  "T"  "A" 
## [46,] "G"  "T"  "C" 
## [47,] "G"  "T"  "G" 
## [48,] "G"  "T"  "T" 
## [49,] "T"  "A"  "A" 
## [50,] "T"  "A"  "C" 
## [51,] "T"  "A"  "G" 
## [52,] "T"  "A"  "T" 
## [53,] "T"  "C"  "A" 
## [54,] "T"  "C"  "C" 
## [55,] "T"  "C"  "G" 
## [56,] "T"  "C"  "T" 
## [57,] "T"  "G"  "A" 
## [58,] "T"  "G"  "C" 
## [59,] "T"  "G"  "G" 
## [60,] "T"  "G"  "T" 
## [61,] "T"  "T"  "A" 
## [62,] "T"  "T"  "C" 
## [63,] "T"  "T"  "G" 
## [64,] "T"  "T"  "T"

Y son 64 porque \(4^3 = 64\), claro.

Las funciones rbinom, runif y rnorm.

Estas tres funciones han aparecido en las clases de teoría y queremos destacar aquí su utilidad para hacer simulaciones, mediante tres ejemplos.

Ejemplo 1:

La probabilidad de que el próximo nacimiento en un hospital sea de un niño es 0.52. En los próximos mil nacimientos que se produzcan en dicho hospital, vamos a calcular la probabilidad de que haya más de 540 niños.

Por un lado podríamos usar la distribución binomial, con \(n = 1000, p=0.52\) para calcular ese valor. Pero lo que vamos a hacer aquí es simular esos 1000 nacimientos en un hospital y ver cuantos niños obtenemos.

rbinom(n = 1, size = 1000, prob = 0.52)
## [1] 525

Ese número podría representar lo que ha ocurrido en un hospital concreto. ¿Y si queremos pensar en 300 hospitales? Usando replicate

numHospitales = 300
(hospitales = replicate(n = numHospitales, {
  rbinom(n = 1, size = 1000, prob = 0.52)
}))
##   [1] 509 521 539 515 523 537 515 509 520 514 539 522 543 492 503 494 522
##  [18] 525 497 502 499 511 504 509 519 515 554 497 522 540 520 524 520 515
##  [35] 501 515 505 517 519 520 510 510 508 511 513 497 531 542 541 519 536
##  [52] 521 515 513 494 493 545 521 510 488 517 528 528 485 478 522 528 534
##  [69] 518 528 525 549 508 526 503 525 545 539 512 524 511 518 520 522 524
##  [86] 527 524 514 513 504 536 560 523 546 527 500 528 481 508 494 495 503
## [103] 505 513 548 516 493 527 504 521 497 523 484 513 494 513 516 519 520
## [120] 508 519 527 507 507 507 520 522 544 523 498 523 535 525 515 506 496
## [137] 526 512 520 515 501 507 517 501 524 508 514 534 520 534 505 548 537
## [154] 531 526 545 522 534 517 496 516 528 508 519 505 515 514 519 517 514
## [171] 519 505 527 511 512 513 514 536 536 498 514 546 528 507 560 523 506
## [188] 506 527 521 515 528 530 538 494 528 547 527 516 495 501 496 508 499
## [205] 486 542 548 524 498 497 527 539 546 520 534 519 510 535 526 514 545
## [222] 541 522 520 508 518 520 538 535 529 512 542 536 509 515 507 520 534
## [239] 521 518 512 501 511 524 503 512 564 565 538 542 505 491 527 561 506
## [256] 515 498 500 526 518 543 517 518 522 492 523 520 523 529 533 506 520
## [273] 551 534 501 514 512 528 510 524 503 523 529 511 537 508 534 512 501
## [290] 511 554 525 526 491 503 552 540 523 509 525

¿En qué fracción de esos hospitales se supera la cifra de 540 niños?

sum(hospitales > 540) / numHospitales
## [1] 0.1

La probabilidad calculada con la binomial es:

1 - pbinom(540, size = 1000, prob = 0.52)
## [1] 0.09715473

Nota: se podría hacer una simulación mucho más “básica”, en la que generamos explícitamente, para cada hospital, la lista de 1000 nacimientos, con probabilidad 0.52 de que sea niño y luego sumamos el número de niños. Pero una vez que sabemos que esa variable es binomial, es mucho más sencillo usar rbinom.

Ejemplo 2:

En el Capítulo 3 vimos este ejemplo:
Elegimos un número \(x\) al azar en el intervalo \([0, 1]\). ¿Cuánto vale la siguiente probabilidad? \[P\left(\frac{1}{3}\leq x\leq \frac{2}{3}\right)\] Para simular esa idea de elegir “un número \(x\) al azar en un intervalo” podemos usar runif. Y como estamos haciendo una simulación, en lugar de elegir un valor de \(x\) vamos a elegir muchos:

n = 1000000
x = runif(n, min = 0, max = 1)

Veamos unos cuantos:

head(x)
## [1] 0.7819173 0.5646282 0.4363120 0.7096982 0.2659250 0.4300874

Y ahora nos preguntamos cuántos de estos puntos cumplen: \[\frac{1}{3}\leq x\leq \frac{2}{3}\] El código es sencillo:

sum((1/3 <= x) & (x <= 2/3)) / n
## [1] 0.333141

Y se aproxima mucho al valor \(1/3\).

Ejemplo 3:

Hemos visto en teoría que Si \(X_1\sim N(\mu_1,\sigma_1)\) y \(X_2\sim N(\mu_2,\sigma_2)\) son variables normales independientes, su suma es de nuevo una variable normal de tipo \[ N\left(\mu_1+\mu_2,\sqrt{\sigma_1^2+\sigma_2^2}\right). \] Vamos a comprobar esto usando rnorm para simular dos muestras de las variables \(X1\) y \(X2\) y representado su función de densidad junto con la densidad teórica de la suma.

n = 100000

mu1 = 5
sigma1 = 0.4
muestraX1 = rnorm(n, mean = mu1, sd = sigma1)

mu2 = -3
sigma2 = 0.6
muestraX2 = rnorm(n, mean = mu2, sd = sigma2)

suma = muestraX1 + muestraX2

plot(density(suma), col="tan", main="Suma de normales", xlab="", ylab="", lwd= 15)
legend("topright", c("Muestra de la suma","Suma teórica"), col=c("tan", "red"), lwd=c(15, 5))
curve(dnorm(x, mean = mu1 + mu2, sd = sqrt(sigma1^2 + sigma2^2)), 
      from = -1, to = 5, add = TRUE, n, col="red", lwd= 5)

Como puede verse, la coincidencia es muy buena.

Vamos a aprovechar para simular la mezcla de estas dos normales.

mezcla = c(muestraX1, muestraX2)
plot(density(mezcla), col="tan", main="Mezcla (mixture) de normales", xlab="", ylab="", lwd= 10)

El resultado de mexclar estas dos normales es una población bimodal (porque sus medias están bastante separadas).

El bucle for de R.

Para aprender más sobre el bucle for, puedes leer la Sección 3 del tutorial 4.

Nuestro primer ejemplo de bucle for recorre los números del 1 al 5. Para cada uno de ellos se usa print y paste0 para escribir una frase en la que aparece el número y su cuadrado.

for(i in 1:5){
  print(paste0("El número i es ", i, " y su cuadrado es ", i^2))
}

El segundo ejemplo es más elaborado. Vamos a fabricar una secuencia de ADN aleatorio y vamos a contar cuántos codones de parada contiene.

bases = c("A", "T", "G", "C")

codonesSTOP = c("TGA", "TAG", "TAA")

secuencia = paste0(sample(bases, 6000, replace = TRUE), collapse = "")

# Usamos substring para ver las primeras 60 bases

substring(secuencia, first = 1, last = 60)

# Inicialmente el número de codones de parada es 0.
cuantosStop = 0

# Usamos el bucle for para recorrer la secuencia.
for(i in 1:(nchar(secuencia) - 2)){
  # Leemos el siguiente codón de la secuencia
  (codon = substring(secuencia, i, i + 2))
  # Analizamos si es un codón de parada
  (esSTOP = codon %in%  codonesSTOP)
  # Y si lo es aumentamos la cuenta en 1
  if(esSTOP){
    cuantosStop = (cuantosStop + 1)
    }
}
cuantosStop

Una pregunta natural es ¿son muchos? En clase hemos visto el ejemplo del genoma de un bacteriófago que era de longitud 5386 bp. ¿Cuántos codones de parada contiene? Los genomas naturales de virus ¿contienen un número mayor o menor de codones de parada que un “genoma aleatorio” como el que hemos fabricado? Este tipo de simulaciones en las que generamos secuencias aleatorias y las comparamos con secuencias naturales son frecuentes en Bioinformática y han permitido avances científicos importantes. Más adelante en el curso veremos como se pueden usar, por ejemplo, para estimar el número de genes en un genoma bacteriano que acabamos de secuenciar.