domingo, 24 de julio de 2011

Expresiones y asignación

### Expresiones y asignación
 
# una expresión es un comando ingresado en la línea de comandos, evaluado por R,
# que devuelve (imprime) un resultado, y luego descartado.
2 + 3 # una expresión
[1] 5 # el resultado de la expresión evaluada (salida)
 
# Una asignación consiste en asignarle un nombre a un objeto que puede ser el
# resultado de una expresión, u otro tipo de objeto.
# en R los caracteres utilizados para asignar objetos son "<-" y "->" y se
# utilizan de la siguiente manera
x <- 2 + 3 # no imprime el resultado, sino que lo asigna a un objeto de nombre "x"
 
# Escribiendo en la línea de comandos el nombre del objeto que acabamos de generar
x # devuelve el resultado de la expresión asignada, que debería aparecer impreso
  # de la siguiente manera:
[1] 5
 
# "->" hace lo mismo que "<-", pero asigna hacia el otro lado:
2 + 3 -> x
 
# comandos incompletos
x <-
2 + 3 # R pone un signo de + en el lugar del prompt normal
      # (el signo de ">" que vemos en la consola previo a ingresar un
      # comando), para indicar que hay que
      # completar el comando.
 
x # si invocamos al objeto x, el resultado va a ser el mismo.
  # R no toma en cuenta los saltos de línea ni los espacios en blanco
  # para interpretar las expresiones (pero los saltos de línea sí son
  # importantes para diferenciar expresiones consecutivas).
 
## El nombre del objeto se puede usar adentro de nuevos comandos. por ejemplo:
x + 1
x * 2
y <- x * 4 + 5
y

Created by Pretty R at inside-R.org

Distribuciones, muestreos y números aleatorios

Dentro de los paquetes base se encuentra "stats", el cual es central para el trabajo con estadísica dentro de R. Una de las utilidades más provechosas tal vez sea la serie de funciones incluidas para trabajar con distribuciones de probabilidad, en particular para generar números (pseudo) aleatorios. En esta sección se van a mostrar algunos generadores de números aleatorios y la forma general de usar estos tipos de funciones.

También vamos a darle especial atención a la función "sample", ya que es de enorme utilidad para hacer muestreos dentro de grupos de objetos.

Por último se hace una mención a los métodos internos de R para generar números aleatorios, así como la forma de "plantar semillas" para reproducir resultados idénticos.

####################################################
## Números aleatorios en R
 
## Podemos ver la lista de distribuciones incluidas en el paquete stats:
?Distributions
 
#### Distibución Normal
?Normal # o
?dnorm  # o
?rnorm  # etc...
## En la ayuda de R se pueden ver las cuatro funciones incluidas para la
## distribución normal:
## dnorm : Función de Densidad de Probabilidad (PDF en inglés).
## pnorm : Función de Distribución de Probabilidad (CDF en inglés; probabilidad acumulada).
## qnorm : Función de Cuantiles.
## rnorm : Generador aleatorio de números con distribución normal.
 
## Las cuatro funciones tienen argumentos que especifican los parámetros
## de la distribución: mean y sd (mu y sigma). Tienen los valores por
## defecto 0 y 1 (lo que se denota por la presencia de los "=").
## Otros parámetros varían ya que son específicos de cada función.
 
## Para crear 1000 números aleatorios con promedio 0 y desvío estándar 1:
x <- rnorm(1000, mean=0, sd=1)
## Claro está que el promedio es de la población y no de cada número en sí...
## Para corroborar que el generador funciona bien:
mean(x) ## Debería ser cercano a 0
sd(x)   ## y este cercano a 1
 
## De acuerdo a la ley fuerte de los grandes números (?), cuánto mayor
## es la cantidad de elementos de x, más cercanos a 0 y 1 son estos
## valores... invitamos al usuario a que explore lo que ocurre con rnorm.
 
## Podemos hacer un resumen de los números generados:
summary(x)
## O un histograma:
hist(x)
## O una curva de densidad:
dens <- density(x)
plot(dens)
## O ambos...
hist(x, freq=FALSE, add=TRUE)
## Nótese el uso de los argumentos "freq" y "add" para cambiar las
## opciones por defecto (así grafica densidades y no borra el gráfico
## anterior). En esta parte no voy a profundizar en lo que refiere
## a herramientas gráficas, pero se retomará el tema más adelante.
 
#### Sobre las otras tres funciones:
## dnorm es la función con forma de campana que la mayoría conocemos:
curve(dnorm(x), from=-3.5, to=3.5)
 
## pnorm es la función que describe P(Y <= X) (probabilidad de que Y sea
## menor o igual a X, siendo X un cuantil e Y una variable aleatoria
## tomada de una distribución normal:
curve(pnorm(x), from=-3.5, to=3.5)
# Naturalmente Y está entre 0 y 1
 
## qnorm matemáticamente es la función inversa a pnorm (f^-1(x), no 1/f(x)).
curve(qnorm(x), from=-1, to=1)
# Se paree a:
plot(sort(rnorm(1e3)))
 
## Puedo combinar dnorm y rnorm para hacer un lindo gráfico:
x <- rnorm(1000)
hist(x, freq=FALSE)
curve(dnorm(x), add=TRUE, col='blue', lwd=2)
lines(density(x), col='red', lwd=2)
 
## La función qqnorm sirve para comparar los cuantiles de un vector
## cualquiera con la distribución de cuantiles esperada para una
## distribución normal:
?qqnorm
 
## Ejemplo:
x <- rnorm(100)
qqnorm(x)
## Usando qqline agregamos una recta que pasa por el primer y tercer
## cuartiles (correspondientes a los percentiles 25% y 75%). 
qqline(x, col='red', lwd=2)
## ¿Qué pasa si hago una prueba similar con un x generado con runif?
## Comparar con la curva de pnorm, ¿Se parecen?¿Tiene lógica?
 
 
 
####################################################
#### Distribución Uniforme
?Uniform # o
?dunif
 
## Cómo dice el nombre, la distribución uniforme es equiprobable para
## todos los valores incluídos en un intervalo dado. Por lo tanto sólo
## necesita un mínimo y un máximo como parámetros.
## Al igual que todos los generadores de números de "stats", el primer
## argumento es el número de elementos, mientras que los otros dos son
## el mín. y máx. mencionados. Para tomar mil valores aleatorios con
## distribución U(0,1):
x <- runif(1000, 0, 1)
hist(x) # Bastante horizontal, ¿no?
 
 
 
####################################################
#### Muestreos aleatorios
?sample
 
## La función sample hace muestreos de los elementos contenidos en el
## primer argumento. Por defecto muestrea sin reposición todos los
## elementos posibles:
sample(1:5)
 
## En la ayuda se pueden ver varios argumentos para modificar este
## comportamiento:
sample(1:5, replace=TRUE) # con reposición
sample(1:5, 2)            # muestrea sólo dos valores, sin reposición
sample(1:5, size=2)       # ídem
sample(1:5, 2, replace=TRUE) # los anteriores combinados
 
## Los mismos resultados se pueden obtener si usamos vectores "character":
sample(c('Juan', 'Mauro', 'Arnold', 'Bill Dance'), 3, replace=TRUE)
sample(letters, 4)
 
## Algo que puede confundir: si le damos solamente un número como entrada,
## puede tener comportamientos inesperados:
sample(10)    # idéntico a pedir "sample(1:10)"
sample(10.3)  # redondea 10.3 --> 10, así que es igual que el anterior
              # es idéntico a pedir: sample(floor(10.3))
sample(-10.3) # el comportamiento es distinto ahora, por ser negativo
 
## Otra opción útil puede ser usar el argumento "prob" de la función.
## Este argumento permite asignar distintos valores de peso a cada
## uno de los elementos del vector a muestrear.
## Por ejemplo, si queremos muestrear los valores 5 10 y 15 pero de tal
## forma que el valor 5 sea 3 veces más frecuente que el 15, y que el
## el 10 sea 2 veces más frecuente que el 15, entonces podemos ingresar
muestreo <- sample(c(5, 10, 15), size=6000, prob=c(3, 2, 1), replace=FALSE)
x <- table(muestreo)
x # x tiene por nombres los valores 5 10 y 15
x['5'] / x['15'] # ¿Es el valor que esperábamos?
x['5'] / x['10'] # ¿Es el valor que esperábamos?
# las comillas indican el nombre del elemento
# en el resultado el 5 indica el nombre del elemento, de forma similar
# a cuando pedimos que muestre x.
 
## Nótese que:
## prob es un sólo argumento, pero eso no equivale a un sólo número, si no
## a un sólo objeto (en este caso un vector, con tres elementos).
## prob puede aceptar valores decimales, incluso pueden asignarse
## probabilidades, aunque da error si alguno es negativo.
## Internamente los valores son siempre convertidos a un valor
## de probabilidad (entre 0 y 1); si el usuario ingresa lo siguigente:
## prob=c(A, B, C)
## Entonces transforma los valores así:
## a = A / (A + B + C)
## b = B / (A + B + C)
## c = C / (A + B + C)
## Y utiliza a, b y c como las probabilidades de que sean muestreados
## los elementos correspondiente del vector a muestrear.
 
 
####################################################
## Semillas
 
## Los números aleatorios generados por computadora, no son verdaderamente
## aleatorios. Los algoritmos internos necesarios son determinísticos
## en último término, pero producen secuencias de valores que aproximan
## propiedades estadísticas aceptablemente similares a las esperadas
## por la distribución uniforme.
## Ver:
# http://www.random.org/;
# http://en.wikipedia.org/wiki/Pseudorandom_number_generator).
?Random
?RNG/
## También existe un paquete llamado "random", para la generación de
## números verdaderamente aleatorios (aunque en muchos casos no son
## necesarios).
 
## Las simulaciones que usan métodos Monte Carlo pueden ser sensibles a
## los problemas. Para el uso de R existe al menos un libro relacionado:
#  http://books.google.com/books/about/Introducing_Monte_Carlo_Methods_with_R.html?id=WIjMyiEiHCsC
 
## Las secuencias generadas por estos algoritmos dependen del valor
## inicial llamado "semilla". Si ingresamos una semilla manualmente
# podemos repetir el resultado de cualquiera de las funciones basadas
## en generadores de números aleatorios (como rnorm, runif y sample):
set.seed(1) ## se puede usar cualquier otro número
sample(10)
# 3  4  5  7  2  8  9  6 10  1
## si ejecutamos de vuelta el comando, va a ser distinto
sample(10)
## a menos que volvamos a "plantar" la misma "semilla":
set.seed(1)
sample(10)
 
## Se obtendrán resultados análogos si usamos rnorm, runif u otros
## generadores de números aleatorios.
 
## Esto puede ser útil si queremos que nuestros resultados sean
## reproducibles *exactamente* en otro momento o lugar.

Created by Pretty R at inside-R.org

sábado, 23 de julio de 2011

break / next

En programación existen ciertas herramientas fundamentales y universales denominadas colectivamente Estructuras de Control. En esta categoría se encuentran los "loops" (lazos, ciclos) y "condicionales".
En esta entrada se va a mostrar el uso del "break" y "next" en R

### Comando "break"
?"break"
## Interrumpe un loop.
## Ejemplo:
for(i in 1:10) {
  print(i)
  break
}
i
## Imprime sólo el primer valor
 
## Combinado con if, sirve para interrumpir en el momento que ocurre
## cierta situación:
## Ejemplo:
## Quiero interrumpirlo cuando i llegue al 7
for(i in 1:10) {
  print(i)
  if(i >= 7) break
}
 
### Comando "next"
?"next"
## Interrumpe la iteración del loop y pasa a la siguiente
for(i in 1:10) {
  if(i == 4) next
  print(i)
}
## Se saltea el número 4
 
for(i in 1:10) {
  if(i %% 2 == 0) next
  print(i)
}
## Se saltea los pares... (¿qué hace el %%?)

if / else


En programación existen ciertas herramientas fundamentales y universales denominadas colectivamente Estructuras de Control. En esta categoría se encuentran los "loops" (lazos, ciclos) y "condicionales".
En esta entrada se va a mostrar el uso del "if" y "else" en R

### Condicionales: if
## Estructura básica:
if(condición) comando
## Se ejecuta el comando si la condición es TRUE
 
## Entonces:
if(3 > 2) print('  :-D  ')
if(3 < 2) print('  :-(  ')
 
## Al igual que vimos antes, se pueden usar las { } para insertar varios
## comandos juntos:
if(3 > 2) {
  print('Es verdadero!')
  print(' :-D ')
}
 
### Usando if / else
### La combinación de if + else sirve para lograr este flujo:
                        / Sí -----> comando1
Se cumple la condición? |
                        \ No -----> comando2

## La estructura básica de if / else es:
if(condición) comando1 else comando2
## Cuando se cumple una condición dada, se ejecuta el comando1, en caso
## contrario, se ejecuta el comando2
## Estructura básica, pero usando {}:
if(condición) {
  comando1 ## Si la condición se cumple (TRUE)
} else {
  comando2 ## Si la condición no se cumple (FALSE)
}
 
## Ejemplo: combinacion de un loop con if/else
## Clasificación de números del 1 al 10 en 2 grupos
clasif <- numeric(10) # Acá voy a guardar el resultado
for(i in 1:10) {
  if(i > 5) { ## Se si se cumple, entonces es "grande"
    texto <- paste(i, 'grande')
    print(texto)
    clasif[i] <- texto ## *1
  } else {    ## Si no se cumple, entonces es "chico"
    texto <- paste(i, 'chico')
    print(texto)
    clasif[i] <- texto ## *2
  }
  # clasif[i] <- texto ## Este comando puede sustituir las lineas *1 y *2
}
clasif
 
## Ídem, pero con ifelse
## esquema conceptual (no correr):
ifelse(condición, comando1, comando2)
x <- 1:10
clasif <- ifelse(x > 5, 'grande', 'chico')
clasif <- paste(x, clasif)
clasif
## Para este tipo de casos sencillos el ifelse es mucho más rápido y simple
## que usar un for + if + else
 
## Hay casos sin embargo, en que usar loops combinados con if y else es
## la mejor opción. Nótese que es necesario hacer un "anidamiento" de
## bloques de comando...
## Ejemplo: Clasificación en tres grupos
clasif <- numeric(10)
for(i in 1:10) {
  if(i > 5) {
    texto <- paste(i, 'grande')
    print(texto)
    clasif[i] <- texto
  } else {
    if(i > 2) {
      texto <- paste(i, 'mediano')
      print(texto)
      clasif[i] <- texto
    } else {
      texto <- paste(i, 'chico')
      print(texto)
      clasif[i] <- texto
    }
  }
}
clasif
 
## Esquemáticamente:
         / Sí -----> "grande"
(i > 5)? |                    / Sí -----> "mediano"
         \ No -----> (i > 2)? |
                              \ No -----> "chico"

Ejemplos extra


## Ejemplo: usar if + break para detener un loop
## Primero necesito la función para evaluar la serie de Leibniz:
leib <- function(n) {
  x <- 0:n
  imp <- 2 * x + 1
  inv <- 1 / imp
  elementos <- inv * (- 1) ^ x
  sum(elementos)
}
 
### Para alcanzar una presición predeterminada
n <- 5     # el n inicial es 5
error <- 0.001 # esta es la presición que queremos
out   <- NULL # en out vamos a ir pegando lo generado **
for(i in 1:100) {
  x   <- leib(n) * 4 # x es el valor de pi estimado
  dif <- abs(x - pi) # dif es |x - pi|
  out <- c(out, dif) # le agrego la diferencia a out
  if(dif < error) { # Si cumple la condición que interesa...
    break          # ... lo interrumpimos con un break
  }
  n <- n + 1
}
plot(out)
 
## ** Usar NULL es una forma de generar un objeto vacío, al cual luego
## se le van a concatenar elementos. En este caso, en cada iteración se
## le concatena un valor nuevo (dif).
 
## Ejemplo: usar un while para lograr lo mismo, pero sin molestarnos
## en poner un número de iteraciones.
#### Para este caso el while es una opción más conveniente
#### ya que a diferencia del for no tiene un límite de iteraciones
#### predeterminada.
#### Para ejecutar lo mismo haciendo un while:
n <- 5
error <- 0.001
dif <- Inf  # dif tiene que estar definido antes de empezar el while
            # porque es necesario para evaluar la condición entre paréntesis
out <- NULL 
while(dif > error) {  # Cáda vez que empieza el loop, se evalúa 
  x <- leib(n) * 4    # la condición entre paréntesis.
  dif <- abs(x - pi) # Si no es verdadera, se corta el loop.
  out <- c(out, dif)
  n <- n + 1        # si no pongo esto hace un loop infinito!
}
plot(out)
 
#### Este proceso se puede implementar en una función.
#### Para esto necesitamos algunos argumentos (n, error y show), aunque
#### show lo ponemos simplemente para que haga opcionalmente un gráfico...
estimaPi <- function(n=5, error=0.01, show=TRUE) {
  dif <- Inf
  while(dif > error) {
    x <- leib(n) * 4
    dif <- abs(x - pi)
    out <- c(out, dif)
    n <- n + 1 # si no pongo esto hace un loop infinito!
  }
  if(show) plot(out# hace el plot si show == TRUE
 
list(estimacion=x, nfinal=n - 1, error=error, dif=out)
}
## Hasta hoar sólo me definió la función, aún no la corrí, para eso:
pi1 <- estimaPi()
str(pi1) # pi1 es una lista
pi1$     # usando la tecla "tab" se pueden ver los objetos contenidos en pi1
abs(pi1$estimacion - pi) ## ¿Qué tan distinto es?
 
## Ejercicio opcional:
## Modificar la función estimaPi para que haga un plot de los
## sucesivos estimados del valor pi.

while

En programación existen ciertas herramientas fundamentales y universales denominadas colectivamente Estructuras de Control. En esta categoría se encuentran los "loops" (lazos, ciclos) y "condicionales".
En esta entrada se va a mostrar el uso del "while" en R

## Estructura básica:
while(condición) comando
## Se ejecuta el "comando" repetidamente, siempre que la "condición" sea TRUE
 
#############################################
## Paréntesis: sobre operadores lógicos
## La "condición" es el resultado de una operación lógica, así que debe ser
## TRUE o FALSE
## Algunas operaciones lógicas sencillas:
3 > 2  # TRUE  (¿tres es mayor que dos?)
3 < 2  # FALSE (¿tres es menor que dos?)
3 == 2 # FALSE (¿tres es igual a dos?)
## Con el ! se cambia un TRUE por FALSE y viceversa (equivale a multiplicar
## por -1):
!(3 > 2)
!(3 == 2)
 
## Se pueden juntar varias condiciones con los operadores AND y OR
 
## AND:
condición1 && condición2 ## ¿Se cumplen ambas condiciones?
## Da TRUE sí y sólo sí condición1 es TRUE y condición2 es TRUE
## Ej.:
3 > 2 && 4/5 < 7/8  ## TRUE
3 > 2 && 4/5 == 7/8 ## FALSE
 
## OR:
condición1 || condición2 ## ¿Se cumple al menos una condición?
## Da TRUE si al menos una de las dos condiciones es TRUE
3 > 2 || 4/5 < 7/8  ## TRUE
3 > 2 || 4/5 == 7/8 ## TRUE
1 > 2 || 4/5 == 7/8 ## FALSE
## Final del paréntesis.
#############################################
 
## Retomando el uso de while, veamos un ejemplo:
x <- 1 ## valor inicial...
while(x < 11) { ## La condición es que x sea menor a 11
  print(x)
  x <- x + 1 ## Aumento el valor de x de a 1 por iteración
}
## imprime los enteros del 1 al 10
## En este caso es igual a un
for(i in 1:10) print(i)
## pero es más útil si no sabemos de antemano cuántas iteraciones
## necesitamos realizar...
 
## Ejemplo:
## Esta serie tiende a cero cuando n tiende a infinito
## f( n ) = f(n - 1) / n
## f( 0 ) = 5
## Va a tomar los valores:
## f(0)=5, f(1)=5, f(2)=5/2, f(3)=(5/2)/3, f(4)=((5/2)/3)/4, ...= 5
## Para calcular los sucesivos valores de la serie, llegando hasta un valor
## menor a 0.001 (o sea: f( n ) < 0.001)
f <- 5 ## Valor inicial
n <- 0
while(f > 0.001) {
  n <- n + 1
  f <- f / n
  print(f)
}
f ## ¿A qué número llegó? (¿qué tan cerca del cero?)
n ## ¿Cuántas iteraciones fueron necesarias?
## La cantidad de iteraciones depende del valor inicial de f
## (probar con 100, 1000)

Created by Pretty R at inside-R.org