3 Ciclos y condicionales

Vamos a analizar los ciclos (se encuentran en el manual de R por si gustas.

##For Al for lo alimentas con una lista de elementos y él (o ella) realizan la operación indicada con todos los elementos de la lista. Es decir el for recorre de uno por uno los elementos de una lista y les aplica una instrucción.

La estructura es como sigue: \[ \begin{equation} \textrm{ for } \Big( \overbrace{i}^\text{Nombre de tu variable} \textrm{ in } \underbrace{1:10}_\text{Lista de elementos}\Big) \overbrace{ \left \{\textrm{Házle algo a $i$} \right \} }^\text{Instrucción que aplicarle a cada elemento} \end{equation} \] Por ejemplo una función que imprima los cuadrados de los primeros 10 números:

for (i in 1:10) {
  print(i^2)
}

[1] 1
[1] 4
[1] 9
[1] 16
[1] 25
[1] 36
[1] 49
[1] 64
[1] 81
[1] 100

##While El while es una instrucción peligrosa. Un while realiza una instrucción indefinidamente mientras se cumpla una condición. La estructura es como sigue: \[ \begin{equation} \textrm{ while } \overbrace{\Big( \textrm{Condición fabulosa} \Big)}^\text{Cosas que tienen que cumplirse para seguir operando} \underbrace{ \left \{ \textrm{Cosas por hacer} \right \}}_\text{Instrucciones} \end{equation} \] Por ejemplo, mientras nuestro número $i$ sea $\leq 5$ le sumamos $1$. Por ejemplo, mientras nuestro número $i$ sea $\leq 5$ le sumamos $1$.

i <- 1
while (i <= 5) {
  print(i)
  i <- i + 1
}

##IMPORTANTE Es muy importante que pusiéramos i <- i + 1 ya que esto obliga a que cada vez que da una vuelta la computadora sume el valor de $i$. Sin ésta instrucción la $i$ valdría siempre 1 y jamás saldríamos del loop (jamás llegaría a valer 5). Si alguna vez quedas atrapado en un loop puedes usar el botón de STOP que tiene R junto a la consola o la tecla escape del teclado.

##If-else (condicionales) En el manual podemos encontrar la sección de condicionales. Los condicionales evalúan el camino que debe de seguir el código según una condición. Los condicionales tienen la siguiente estructura:

\[ \begin{align*} if (Condición) \left \{ \right. \\ \\ \text{Cosas por hacer si ocurre la condición} \\ \\ \left. \right \} else \left \{ \right. \\ \\ \text{Cosas por hacer si no se cumple} \\ \\ \left. \right \} \\ \\ \end{align*} \]

¡Vamos al ejemplo!

Asignemos, primero, el valor de 2 a $i$:

i <- 2

Luego pongamos el condicional

if (i == 5) { # Nota el doble signo de igual '=='.
  i <- i^2 # Si sólo pones uno, R hace que i=5.
}

Veamos cuánto vale $i$

print(i)

[1] 10

Hagamos ahora $i = 5$ y veamos qué pasa cuando atraviesa el condicional:

i <- 5
if (i == 5) {
  i <- i^2
}
print(i)

[1] 25

Hagamos ahora un condicional más complicado: pongamos que si $i = 5$ entonces haga $i^2$ pero si $i \neq 5$ entonces al valor de $i$ le sume $1$:

i <- 1
if (i == 5) {
  i <- i^2
} else {
  i <- i + 1
}
print(i)

[1] 2

##And-or (Operadores lógicos)

¿Qué pasa si tienes varias condiciones que necesitas se cumplan a la vez? ¡Para eso están los operadores lógicos! .

##And

Por ejemplo, supongamos queremos usar dos condiciones dentro del if. ¡Es muy fácil! Basta con escribir Condición 1 & Condición 2.

Ejemplo:

i <- 1
j <- 7
if ((i == 1) & (i < j)) {
  i <- i + j
}
print(i)

[1] 8

Por otro lado si ahora hacemos $i > j$:

i <- 11
j <- 7
if ((i == 1) & (i < j)) {
  i <- i + j
}
print(i)

[1] 11

Y si hacemos que $i = 1$ pero $i > j$:

i <- 1
j <- 0
if ((i == 1) & (i < j)) {
  i <- i + j
}
print(i)

[1] 1

##Or

El or se usa en el caso de que querramos que se cumpla al menos una de dos condiciones del if. Es decir, si tenemos dos condiciones el or se cumple cuando se cumple una de ellas o cuando se cumplen ambas. Por ejemplo:

Cuando $i = 1$ con $i > j$:

i <- 1
j <- 1 / 2
if ((i == 1) | (i < j)) {
  i <- i + j
}
print(i)

[1] 1.5

Cuando $i \neq 1$ pero $ i < j $:

i <- 21
j <- 22
if ((i == 1) | (i < j)) {
  i <- i + j
}
print(i)

[1] 43

O bien cuando $i \neq 1$ e $i > j$

i <- 121
j <- 22
if ((i == 1) | (i < j)) {
  i <- i + j
}
print(i)

[1] 121

La siguiente tabla resume cuándo se cumplen las condiciones:

\[ \begin{array}{ccccc} Condición 1 & Condición 2 & And & Or \\ \hline Si & Si & Si & Si \\ Si & No & No & Si \\ No & Si & No & Si \\ No & No & No & No \\ \end{array} \]

##Ejercicio 1 Crea los comandos necesarios para calcular la media y desviacion estandar del siguiente vector:

numeros <- c(
  7.65688984, 0.45416281, -0.53197482, -11.68901517,
  -0.22092715, -6.65860576, -0.96411401, -0.04875882,
  -0.88076032, -9.47716275, 12.48699956, 58.37690942,
  0.75332369, -0.07644519, -0.47168251, 0.04574367,
  0.21158367, 6.57919350, 1.61654489, -144.28602691
)

Tus resultados deberían ser:

[1] "Media: -4.356206118"

[1] "Desviación: 35.8220144835353"

Nota Para el ejercicio puedes usar cualquier función de R excepto: mean, sd, var.

##RECORDATORIO Por si no lo recuerdas, aquí están las definiciones de media y desviacion estándar. Si bien no es la única forma pues ¡hay varias definiciones equivalentes!.

Aquí consideraremos $X$ como una variable con observaciones para $N$ individuos. Es decir: $X = (x_1,x_2,\cdots,x_n)$.

##Media \[ \begin{equation} \textrm{Media de }X = \frac{x_1 + x_2 + \cdots + x_n}{n} \end{equation} \] ##Desviación estándar \[ \begin{equation} \textrm{Desviación Estándar de }X = \sqrt{\textrm{Media de }X^2 - \Big(\textrm{Media de }X\Big)^2} \end{equation} \]

##Ejercicio 2 Sin correr el siguiente pedazo de código en R, estima cuánto valdrá $k$ al final:

k <- 3
for (i in 1:6) {
  if (i > k || k == 3) {
    k <- k^2
  } else if (i == 3 & k == 7) {
    k <- k - 2
  } else if (k > i & i < 5) {
    k <- k * i / 2
  } else if (k > i & i >= 5) {
    k <- k + 1
  } else {
    k <- k / 2
  }
}

##Ejercicio 3 Un grupo de investigadores tienen tres vectores de datos sobre individuos: sexo, edad y exposición (horas) a humo de tabaco expo.

sexo <- c(
  "Hombre", "Mujer", "Mujer", "Mujer", "Hombre",
  "Mujer", "Hombre", "Hombre"
)
edad <- c(28, 12, 77, 32, 46, 53, 17, 20, 88)
expo <- c(1, 0, 1.5, 2.2, 2, 5, 1.01, 3.2)

Ellos saben que por cada hora de exposición el riesgo relativo de enfermedad cardiovascular es de $1.025$ para hombres menores a 45 y $1.032$ para mujeres de la misma edad. Para mayores de 45, el riesgo es $1.052$ en caso de hombres y $1.066$ en caso de mujeres.

Los investigadores hicieron el siguiente código para estimar los riesgos de cada uno de los individuos. Ayúdalos a que su código funcione:

#Hay n personas: para cada una hay que calcular su riesgo 

n      <- length(sexo)
riesgo <- c()
while (persona < n){
  
  #Checar la edad de la persona
  if (edad[persona] < 45){
    
    #Checar el sexo
    if (sexo[persona] = Hombre){
      
      riesgo[persona] <- expo[persona]*1.025
      
    } else {
      
      riesgo[persona] <- expo[persona]*1.032
      
    }
    
  } else {
    
    #Checar el sexo
    if (sexo[persona] = Hombre){
      
      riesgo[persona] <- expo[persona]*1.052
      
    } else {
      
      riesgo[persona] <- expo[persona]*1.066
      
    }
    
  }
  
  
  
}

Para que cheques que funcione, te dejo la respuesta. El riesgo es:

[1] 1.02500 0.00000 1.59900 2.27040 2.10400 5.33000 1.03525

##Advertencias y otras cosas poco intuitivas

Es importante entender cómo funcionan las computadoras para poder simular (y entender los problemas de la simulación). Aunque los números son infinitos, las computadoras no tienen una cantidad infinita de dígitos. Por ejemplo, nosotros (humanos) podemos representar: \[ \begin{equation} 1 - 0.000000001 = 0.99999999 \end{equation} \] La computadora no puede hacerlo:

1 - 0.000000001

[1] 1

Tampoco puede representar números muy grandes:

exp(1000)

[1] Inf

Mientras que para los humanos no hay ``un número positivo más chico’’ (si dices, por ejemplo, que $0.00000000001$ es el más chico de todos los positivos (no cero), siempre puedes dividirlo entre $2$: $0.00000000001/2$ y obtener un número más pequeño) para las computadoras sí hay. Eso quiere decir que cuando hacemos una operación la computadora NO da la respuesta correcta sólo su mejor aproximación. A veces su mejor aproximación es la respuesta correcta:

sqrt(100)

[1] 10

Otras veces está cerca:

sqrt(12345678.12345678^2)

[1] 12345678

Pero…

##Donde fallan estas cosas

Intuitivamente, los decimales que nos acabamos de comer en el inciso anterior no importan ¡son simples decimales! El siguiente ejemplo muestra que sí importan.

Este ejemplo calcula una función recursivamente. ¿Puedes explicar qué estamos haciendo?

ejemplo <- c()
a <- 2.701
for (i in 1:100) {
  if (i == 1) {
    ejemplo[i] <- 10
  } else {
    ejemplo[i] <- ejemplo[i - 1] +
      a * ejemplo[i - 1] * (1 - ejemplo[i - 1] / 100)
  }
}

print(ejemplo[100])

[1] 125.8633

El ejemplo anterior resulta en un maravilloso resultado de ejemplo[100]. Redondeemos a dos decimales el valor de $a$ para que sea $2.70$. Intuitivamente, el valor debería estar cerca y ser ciento y tantos. Pues no…

ejemplo2 <- c()
a <- 2.70
for (i in 1:100) {
  if (i == 1) {
    ejemplo2[i] <- 10
  } else {
    ejemplo2[i] <- ejemplo2[i - 1] +
      a * ejemplo2[i - 1] * (1 - ejemplo2[i - 1] / 100)
  }
}

print(ejemplo2[100])

[1] 87.09381

¡Resulta que con cambiar un decimal, el resultado cambió hasta ejemplo2[100]! La gráfica siguiente muestra como varían los valores coincidiendo al inicio y alejándose después:

##Números pseudoaleatorios Para simular necesitamos generar números aleatorios. La única forma que tenemos de hacerlo (actualmente) es por medio de isótopos radiactivos que decaen aleatoriamente. ¡Si tienes uno guardado por ahí es el momento de usarlo!

Como no es muy bueno que tengamos por ahí material radiactivo, los matemáticos han generado números que se conocen como pseudoaleatorios. Estas son funciones (como la del apartado anterior) que si conoces el valor inicial (en el caso pasado, $a$) las funciones son tan alocadas que los números que resultan de ella ‘’parecen aleatorios’’.

##Ejercicio 4 Considera los siguientes dos fragmentos de código. Analiza los resultados. ¿Cuál de ellos es un mejor generador pseudoaleatorio y por qué?

a <- 2
aleatorio <- c()
for (i in 1:100) {
  aleatorio[i] <- a + i
}

[1] 3 4 5 6 7 8

aleatorio2 <- c()
for (i in 1:100) {
  if (i == 1) {
    aleatorio2[i] <- 0.9
  } else {
    aleatorio2[i] <- aleatorio2[i - 1] +
      2.81 * aleatorio2[i - 1] * (1 - aleatorio2[i - 1] / 17)
  }
}

[1]  0.900000  3.295112 10.759652 21.858156  4.305520 13.339891

##Aleatoreidad en R

Nuestra semilla

Si un día despiertas con ganas de tener 10 números aleatorios, en R ¡puedes hacerlo!:

runif(10)

 [1] 0.47 0.34 0.33 0.48 0.25 0.31 0.03 0.40 0.72 0.22

Además puedes especificar la distribución. Por ejemplo, ahora sacaremos 7 números aleatorios de una normal estándar:

rnorm(7)

[1]  1.08 -0.72 -0.33 -0.57  0.58  0.55  0.40

O bien 8 valores de una exponencial con parámetro 5:

rexp(8, rate = 5)

[1] 0.18 0.14 0.02 0.44 0.03 0.23 0.76 0.25

Vamos entonces a simular 1,000 números alearorios de una normal estándar:

y <- rnorm(1000)

Un comando muy útil es la función summary que resume los cuantiles principales de la distribución así como el mínimo, el máximo y el promedio.

summary(y)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-3.181468 -0.654028 -0.002774 -0.016848  0.644352  3.633669

Al momento de pedir ayuda para el comando rnorm nos podemos dar cuenta de otras funciones interesantes relacionadas con la normal:

?rnorm

Podemos, estimar, por ejemplo, la densidad acumulada de una normal estándar en el 0. Es decir, ¿a qué percentil de una normal corresponde el 0?

pnorm(0)

[1] 0.5

Podemos hacer lo mismo para los números entre -10 y 10:

pnorm(-10:10)

 [1] 7.619853e-24 1.128588e-19 6.220961e-16 1.279813e-12 9.865876e-10
 [6] 2.866516e-07 3.167124e-05 1.349898e-03 2.275013e-02 1.586553e-01
[11] 5.000000e-01 8.413447e-01 9.772499e-01 9.986501e-01 9.999683e-01
[16] 9.999997e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
[21] 1.000000e+00

Igualmente podemos graficar cómo se ven:

plot(-10:10, pnorm(-10:10))

O bien graficar la función de distribución de una normal entre -10 y 10:

plot(-10:10, dnorm(-10:10))

¿Notas cómo nos faltan números en medio? Es porque el comando -10:10 recorre los números entre -10 y 10 de 1 en 1.

-10:10

 [1] -10  -9  -8  -7  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6   7   8
[20]   9  10

Para hacer más refinada la cantidad de puntos podemos hacer ahora una nueva secuencia pero yendo de 0.1 en 0.1:

x <- seq(-10, 10, 0.1)

Puedes ver cómo se guardaron estos valores (este documento nada más muestra los primeros 9 porque es desperdiciar mucho espacio poner los 201 valores que hizo R)

[1] -10.0  -9.9  -9.8  -9.7  -9.6  -9.5  -9.4  -9.3  -9.2

¡La gráfica ahora se ve genial!

plot(x, dnorm(x))

##Las semillas En el apartado anterior dijimos que los números de R no eran aleatorios sino pseudoaleatorios y que estos se generaban por medio de una función. Cuando estamos haciendo investigación con simulaciones, para que nuestro estudio sea reproducible, aunque usemos números aleatorios, debemos usar siempre los mismos. La semilla se asegura de ello. Para poner una semilla usa el comando set.seed y pon dentro un entero.

set.seed(1234)

Obtengamos un número aleatorio normal:

rnorm(1)

[1] -1.207066

rnorm(1)

[1] 0.2774292

Volvamos a poner la semilla y saquemos un tercero:

set.seed(1234)
rnorm(1)

[1] -1.207066

¿Notas que es el mismo número que al inicio?

##Ejercicio 5

Considera una población cuyo peso se distribuye normal con media 69 y desviación estándar 4.7. Esa misma población tiene una altura normal con media 1.8 $m^2$ y desviación estándar de $0.05$. Simula su índice de masa corporal. Calcula la media y desviación estándar del mismo. ¡No te olvides de usar una semilla!

##Ejercicio 6

Un grupo de investigadores ha decidido que el índice de masa corporal tiene una distribución Cauchy y han simulado el índice de masa corporal como sigue:

set.seed(6207)
# Simular IMC
imc <- rcauchy(100, 25, 0.9)

# Calcular la media
mean(imc)

El código es correcto. Pero la hipótesis de la distribución Cauchy no. Calcula la media varias veces usando diferentes semillas ¿Cuál es el problema? ¿Ocurre lo mismo si calculas la mediana?