Simulación

Uso básico de R

Esta sección sirve como repaso o introducción, por lo cual participantes que ya tienen familiaridad en R y lo han utilizado recientemente no necesitan revisar esta parte sino pueden proceder directamente a la sección de paralelismo.

Instalación

R es una herramienta muy útil y versátil para el cómputo científico, en especial todo lo relacionado a probabilidad y estadística. Funciona en Windows, Linux & Mac OS. Está disponible en línea de forma gratuita en www.r-project.org.

El material de estudio en esta página incluye las instrucciones requeridas para diferentes operaciones necesarias. Conviene hacer esta ventana de navegador delgada para que ocupe la mitad de la pantalla por máximo y colocar a su lado otra ventana que tenga R ejecutándose.

Cada código de R (que se muestra en esta página en letra azul) conviene copiarlo a R y probarlo uno mismo, modificando y explorando hasta que quede claro qué y por qué está sucediendo en la herramienta.

Al salir de R (con la instrucción quit(), la herramienta pregunta si se desea guardar la sesión (dice Save workspace image? [y/n/c]:. Si se guarda (presionando la tecla "y" por "yes"), las mismas variables ya definidas seguirán disponibles la próxima vez que R se inicia desde esa misma carpeta (mientras "c" significa "cancelar y no salir" y "n" significa "salir sin guardar"). Lo de guardar datos de sesiones es práctico durante el curso ya que algunos datos se usan en múltiples apartados y así no se tendrán que redefinir entre sesiones de trabajo.

La asignación de valores a variables en R y la consulta de ellas se realiza de la siguiente forma:

> X <- 3
> X
[1] 3
> X = 4
> X
[1] 4
> 3 -> X
> X
[1] 3

El texto de arriba incluye la salida de R, por lo cual incluye los símbolos ">" que marcan los inicios de instrucciones en R. No se teclean esos símbolos. El usuario escribe lo de X <- 3 mientras R ya produjo el ">" para indicar que está listo para recibir instrucciones. Solamente se ocupa escribir la parte en letra azul.

La creación de un vector en R y la consulta de ello se realiza de la siguiente forma:

> datos <- c(4, 2, 4, 5)
> n <- length(datos)
> datos
[1] 4 2 4 5
> n
[1] 4
> min(datos)
[1] 2
> max(datos)
[1] 5
> sum(datos)
[1] 15
> ord <- sort(datos)
> ord[1]
[1] 2
> ord[n]
[1] 5

Otras rutinas útiles para analizar vectores son las siguientes:

> data <- c(1, 4, 2, 4, 2, 5, 6, 7, 4, 76, 3, 2, 5, 6, 7)
> mean(data)
[1] 8.933333
> median(data)
[1] 4
> summary(data)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
1.000   2.500   4.000   8.933   6.000  76.000
> quantile(data, 0.9)
90%
7
> table(data)
data
1  2  3  4  5  6  7 76
1  3  1  3  2  2  2  1
> var(data)
[1] 347.781
> sd(data)
[1] 18.64889
> otra <- c(1, 5, 3, 5, 3, 65, 4, 6, 3, 6, 3, 56, 4)
> length(otra)
[1] 13
> length(data)
[1] 15
> otra <- c(otra, 2, 4)
> cor(data, otra)
[1] -0.08119101

Para crear una matríz, se ocupa la rutina matrix. Toma como argumentos los elementos por acomodar por columnas (la opción byrow=TRUE cambia el orden de llenado) y las dimensiones. Por ejemplo, M <- matrix(rep(0, 9), nrow=3, ncol=3) produce una matriz tres por tres con puros ceros:

> M <- matrix(rep(0, 9), nrow=3, ncol=3)
> M
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

Para comparaciones entre vectores, las subrutinas all y any sirven para determinar si todos o algunos de los elementos son iguales:

> a <- c(1, 2, 3)
> b <- c(1, 2, 4)
> a == b
[1] TRUE TRUE FALSE
> all(a == b)
[1] FALSE
> any(a == b)
[1] TRUE

La aritmética funciona en grandes rasgos como uno esperaría, con el detalle que las operaciones vectoriales y matriciales de álgebra lineal contienen símbolos de porcentaje mientras los que operan por elementos tienen la notación estándar:

> a <- 3
> b <- 4
> c <- c(5, 6, 7)
> d <- c(8, 10, 12)
> a + b
[1] 7
> a - b
[1] -1
> a * b
[1] 12
> a / b
[1] 0.75
> a^b
[1] 81
> b^a
[1] 64
> a**b
[1] 81
> a + c
[1]  8  9 10
> c + d
[1] 13 16 19
> c * d
[1] 40 60 84
> c %*% d
[,1]
[1,]  184
> e <- matrix(c(1, 2, 3, 4), nrow=2)
> e
[,1] [,2]
[1,]    1    3
[2,]    2    4
> f <- matrix(c(2, 4, 6, 8), nrow=2)
> f
[,1] [,2]
[1,]    2    6
[2,]    4    8
> e + f
[,1] [,2]
[1,]    3    9
[2,]    6   12
> e - f
[,1] [,2]
[1,]   -1   -3
[2,]   -2   -4
> e * f
[,1] [,2]
[1,]    2   18
[2,]    8   32
> e %*% f
[,1] [,2]
[1,]   14   30
[2,]   20   44

El redondeo de enteros a decimales se hace con las tres reglas: hacia abajo con la función piso floor, hacia arriba con la función techo ceiling y al entero más cercano con round. La división entera se logra con %/% y el residuo (también llamado modulo, que en R funciona hasta con valores decimales) con %%:

> d <- 7.3
> floor(d)
[1] 7
> ceiling(d)
[1] 8
> round(d)
[1] 7
> d / 2
[1] 3.65
> d %/% 2
[1] 3
> d %% 2
[1] 1.3

Las funciones matemáticas tienen sus nombres típicos y están disponibles en R básico sin cargar paquetes adicionales:

> x <- 123.4
> sqrt(x)
[1] 11.10856
> exp(x)
[1] 3.907861e+53
> sin(x)
[1] -0.7693905
> cos(x)
[1] -0.6387787
> tan(x)
[1] 1.204471
> tanh(2.5)
[1] 0.9866143
> log(x)
[1] 4.815431
> log(x, 10)
[1] 2.091315
> log(x, 2)
[1] 6.947199

Para permutaciones y combinaciones, existen subrutinas para el factorial y la coeficiente binomial:

> factorial(5)
[1] 120
> choose(9, 3)
[1] 84
> factorial(9) / (factorial(6) * factorial(3))
[1] 84

Para manipular secuencias, existen las siguientes subrutinas prácticas:

> s <- seq(4, 100, 2)
> s
[1]   4   6   8  10  12  14  16  18  20  22  24  26  28  30  32  34  36  38  40
[20]  42  44  46  48  50  52  54  56  58  60  62  64  66  68  70  72  74  76  78
[39]  80  82  84  86  88  90  92  94  96  98 100
> r <- rep(13, 20)
> r
[1] 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
> i <- 2:12
> i
[1]  2  3  4  5  6  7  8  9 10 11 12
> c <- cut(s, breaks=4)
> table(c)
c
(3.9,28]  (28,52]  (52,76] (76,100]
      13       12       12       12
> tail(s)
[1]  90  92  94  96  98 100
> tail(s, 10)
[1]  82  84  86  88  90  92  94  96  98 100
> head(s)
[1]  4  6  8 10 12 14
> head(s, 8)
[1]  4  6  8 10 12 14 16 18

Para expresiones condicionales, repeticiones, y lógica, sirve lo siguiente:

> x <- 3
> y <- NULL
> if (x > 5) { y <- x * 3 } else { y <- 6 - x }
> y
[1] 3
> z <- NULL
> if (x %% 2 == 1 & y %% 2 == 1) { z <- 2 } else { z <- 3 }
> z
[1] 2
> if (z != 1) { print("no es uno") }
[1] "no es uno"
> if (x == 2 | y == 2 | z == 2) { print("hay un dos") } else { print("no hay un dos") }
[1] "hay un dos"
> for (i in 1:5) { print(2**i) }
[1] 2
[1] 4
[1] 8
[1] 16
[1] 32
> while (x > 0) { print("quito uno"); x <- x - 1 }
[1] "quito uno"
[1] "quito uno"
[1] "quito uno"
> x
[1] 0

Subrutinas

Cuando algún cálculo ocupa ser reutilizado en múltiples ocasiones, posiblemente con algunas variaciones, conviene encapsularlo en una subrutina como una función:

nombre <- function(p1, p2, ... ){...; return(resultado)}

Conviene probar la creación y llamadas a subrutinas propias — hay un buen tutorial en inglés [1].

> potdos <- function(x) 2**x
> potdos(3)
[1] 8
> eucl <- function(x1, y1, x2, y2) { dx <- x1 - x2; dy <- y1 - y2; return(sqrt(dx^2 + dy^2)) }
> eucl(1, 2, 3, 4)
[1] 2.828427

Además del ; se puede utilizar un salto de línea para dividir el contenido de un bloque de instrucciones en una subrutina o dentro de una expresión condicional tipo if o un ciclo for o while.

Números pseudoaleatorios

Generación pseudoaleatoria refiere a la creación de datos que parecen ser al azar. Verdaderamente no lo son si se generan de forma mecánica/aritmética. El sitio web random.org proporciona datos que están basados en bits generados por ruido atmosférico que son "más aleatorios" que cualquier cosa que se genere de forma mecánica.

Para números pseudoaleatorios uniformemente distribuidos, se usa runif, con parámetros opciones en el caso que no se quiera que sean entre cero y uno sino de otro rango de valores. Para crear valores de verdad o falso, 50–50, por ejemplo, se puede condicionar al valor uniforme. La subrutina sample sirve para muestreo y permutaciones.

> runif(1)
[1] 0.8386876
> runif(5)
[1] 0.7410066 0.5553238 0.8469588 0.4764683 0.7612842
> runif(2, 4, 10)
[1] 8.608724 7.326382
> runif(1) < 0.5
[1] TRUE
> runif(1) < 0.5
[1] FALSE
> runif(1) < 0.5
[1] TRUE
> runif(5) < 0.5
[1] TRUE TRUE FALSE FALSE TRUE
> sample(1:10, 3)
[1] 7 4 3
> sample(1:10, 3)
[1] 3 6 1
> sample(1:10)
[1] 5 4 3 1 9 6 10 2 8 7
> sample(c(2, 4, 6), 10, prob=c(0.1, 0.2, 0.7), replace=TRUE)
[1] 6 6 6 6 2 2 6 2 4 6

R es capaz de producir números pseudoaleatorios en una gran gama de distributiones (c.f. catálogo), además de contar con las funciones de densidad y cumulativas de las distribuciones. Ejemplos útiles incluyen las distribución normal y exponencial.

> rnorm(10) # media cero, desv. est. 1
[1] -0.82017580 -0.17005184 0.19908914 0.53885462 -0.13823311 0.08619352
[7] -0.55169478 0.43233720 1.64800753 -1.30523199
> rnorm(10, mean=5, sd=1.3)
[1] 3.463890 3.637588 3.634250 6.278641 6.297803 5.133415 6.750187 4.553651
[9] 4.419553 3.758398
> rexp(2)
[1] 2.1036209 0.3027257
> rexp(2, rate=3)
[1] 0.1174889 0.3573015

Conjuntos de datos

Para traer datos a R, lo más conveniente es utilizar archivos de texto, o importar desde hojas de cálculo en el formato CSV (comma separated values), aunque en R se puede utilizar cualquier separador.

Para probar, colocamos los siguientes datos en un archivo llamado ejemplo.csv con un editor de texto básico tipo notepad o emacs:

Nombre Estatura Gatos Perros
Elisa 184 2 0
Gonzalo 177 0 1
Jorge 180 1 1
Mariana 158 0 0
Carlos 164 1 0

> datos <- read.csv("ejemplo.csv", header=TRUE, sep=" ", stringsAsFactors=FALSE)
> datos
   Nombre Estatura Gatos Perros
1   Elisa      184     2      0
2 Gonzalo      177     0      1
3   Jorge      180     1      1
4 Mariana      158     0      0
5  Carlos      164     1      0
> datos$Estatura <- as.numeric(datos$Estatura)
> datos$Gatos <- as.numeric(datos$Gatos)
> datos$Perros <- as.numeric(datos$Perros)
> sum(datos$Gatos)
[1] 4
> sum(datos$Gatos) > sum(datos$Perros)
[1] TRUE
> datos$Gatos == datos$Perros
[1] FALSE FALSE  TRUE  TRUE FALSE
> datos[datos$Nombre == "Elisa",]
  Nombre Estatura Gatos Perros
1  Elisa      184     2      0
> prom <- mean(datos$Estatura)
> datos[datos$Estatura > prom,]
   Nombre Estatura Gatos Perros
1   Elisa      184     2      0
2 Gonzalo      177     0      1
3   Jorge      180     1      1
> datos[datos$Estatura > prom,]$Gatos
[1] 2 0 1
> datos <- rbind(datos, c("Tania", 175, 0, 2))
> datos
   Nombre Estatura Gatos Perros
1   Elisa      184     2      0
2 Gonzalo      177     0      1
3   Jorge      180     1      1
4 Mariana      158     0      0
5  Carlos      164     1      0
6   Tania      175     0      2
> datos$Estudiante <- c(FALSE, TRUE, TRUE, TRUE, FALSE, FALSE)
> datos
   Nombre Estatura Gatos Perros Estudiante
1   Elisa      184     2      0      FALSE
2 Gonzalo      177     0      1       TRUE
3   Jorge      180     1      1       TRUE
4 Mariana      158     0      0       TRUE
5  Carlos      164     1      0      FALSE
6   Tania      175     0      2      FALSE
> datos[3,]
  Nombre Estatura Gatos Perros Estudiante
3  Jorge      180     1      1       TRUE
> datos[2:4,]
   Nombre Estatura Gatos Perros Estudiante
2 Gonzalo      177     0      1       TRUE
3   Jorge      180     1      1       TRUE
4 Mariana      158     0      0       TRUE
> datos[2:4,3]
[1] 0 1 0
> datos$Estudiante <- as.factor(datos$Estudiante)
> datos$Estudiante
[1] FALSE TRUE  TRUE  TRUE  FALSE FALSE
Levels: FALSE TRUE

Gráficas sencillas

R contiene bastantes funcionalidades ya en su forma básica para graficar información. La subrutina plot crea una gráfica (será de puntos si no se solicita que sea de líneas con type="l"). Se le pueden agregar cosas en ese mismo dibujo llamando a points (para agregar puntos), lines (para líneas entre puntos definidos) o abline para rectas. El título se coloca con el argumento main:

> x <- 1:20
> y <- sin(x)
> plot(x, y)
> plot(x, y, type="l")
> lines(x, cos(x), col="red")
> abline(h=0, lwd=5, col="green")
> plot(x, y)
> points(x + 1, cos(x + 1), pch=15, col="green")
> plot(x, y, main="Texto arriba", xlab="Etiqueta", ylab="Otra etiqueta")

Se puede realizar la redirección de gráficas a archivos, por ejemplo para compartirlos o cuando uno no está trabajando en un ambiente gráfico que pueda abrir ventanas:

> png(filename="salida.png")
> plot(x, y, pch=16)
> graphics.off()

Si no se define un archivo de salida para la gráfica y R no está ejecutándose en un ambiente gráfico, se crea por automático un archivo llamado Rplots.pdf con las gráficas generadas.

Diagramas caja-bigote

> d1 <- c(10, 32, 34, 24, 49, 42, 89)
> d2 <- c(30, 23, 45, 24, 75, 34, 12, 56, 33)
> boxplot(d1)
> boxplot(d2)
> boxplot(d1, d2)

Histogramas

Existe una rutina en R que se llama hist. Llamando a hist(datos), R produce lo que se le hace un histograma y lo dibuja. Se le puede indicar a la rutina a cuántas cubetas realizar la división: hist(datos, breaks=5), lo que produce cinco divisiones dentro del intervalo y por ende seis barras.

> datos <- c(1.6, 4.6, 2.6, 3.6, 5.6, 6.6, 3.5, 2.2, 4.4, 5.2, 5.4, 7.6, 5.8, 4.4, 6.4)
> hist(datos)
> hist(datos, breaks=2)
> histograma <- hist(datos)
> histograma
$breaks
[1] 1 2 3 4 5 6 7 8

$counts
[1] 1 2 2 3 4 2 1

$density
[1] 0.06666667 0.13333333 0.13333333 0.20000000 0.26666667 0.13333333 0.06666667

$mids
[1] 1.5 2.5 3.5 4.5 5.5 6.5 7.5

$xname
[1] "datos"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
> histograma$density
[1] 0.06666667 0.13333333 0.13333333 0.20000000 0.26666667 0.13333333 0.06666667
> sum(histograma$density)
[1] 1
> plot(histograma, freq=FALSE)

Instalación de paquetes adicionales

Al utilizar la herramienta install.packages para incorporar una nueva funcionalidad a R, es necesario determinar de dónde se van a descargar:

--- Please select a CRAN mirror for use in this session ---
HTTPS CRAN mirror

1: 0-Cloud [https]                 2: Algeria [https]              
3: Australia (Canberra) [https]    4: Australia (Melbourne) [https]
5: Australia (Perth) [https]       6: Austria [https]              
7: Belgium (Ghent) [https]         8: Brazil (PR) [https]          
9: Brazil (RJ) [https]            10: Brazil (SP 1) [https]        
11: Bulgaria [https]               12: Canada (MB) [https]          
13: Chile 1 [https]                14: Chile 2 [https]              
15: China (Beijing) [https]        16: China (Hefei) [https]        
17: China (Lanzhou) [https]        18: Colombia (Cali) [https]      
19: Czech Republic [https]         20: Denmark [https]              
21: Estonia [https]                22: France (Lyon 1) [https]      
23: France (Lyon 2) [https]        24: France (Marseille) [https]   
25: France (Montpellier) [https]   26: France (Paris 2) [https]     
27: Germany (Göttingen) [https]    28: Germany (Münster) [https]    
29: Greece [https]                 30: Iceland [https]              
31: India [https]                  32: Indonesia (Jakarta) [https]  
33: Ireland [https]                34: Italy (Padua) [https]        
35: Japan (Tokyo) [https]          36: Malaysia [https]             
37: Mexico (Mexico City) [https]   38: New Zealand [https]          
39: Norway [https]                 40: Philippines [https]          
41: Russia (Moscow) [https]        42: Serbia [https]               
43: Spain (A Coruña) [https]       44: Spain (Madrid) [https]       
45: Sweden [https]                 46: Switzerland [https]          
47: Taiwan (Chungli) [https]       48: Turkey (Denizli) [https]     
49: Turkey (Mersin) [https]        50: UK (Bristol) [https]         
51: UK (Cambridge) [https]         52: UK (London 1) [https]        
53: USA (CA 1) [https]             54: USA (IA) [https]             
55: USA (IN) [https]               56: USA (KS) [https]             
57: USA (MI 1) [https]             58: USA (OR) [https]             
59: USA (TN) [https]               60: USA (TX 1) [https]           
61: Vietnam [https]                62: (HTTP mirrors)


Selection: 37

Si el sistema operativo suele abrir una ventana para la selección aunque el usuario prefiere utilizar el terminal, se puede forzar esto con chooseCRANmirror(graphics=FALSE).

Al instalar más paquetes, una vez que ya se haya pasado por este paso, R ya recordará la selección del servidor (hasta que se cierre la sesión).

Actualizado el 16 de noviembre del 2022.
https://satuelisa.github.io/simulation/tutorial.html