0. SCRIPTS DE R

1. EL LENGUAJE DE PROGRAMACIÓN DE R

1.1. Generalidades

  • Es un lenguaje más, orientado a objetos, e interpretado (no compilado), pero diseñado por expertos en informática y estadística.
  • Constantes: numéricas (3, 0.1, -2.18,…), de cadena de texto ("a", 'dos palabras',…), lógicas (TRUE y FALSE) y otras erspeciales (Inf, NA, NaN, NULL,…)
  • Operadores lógicos (&, |, !, no confundir con otros lenguajes)
  • Operador de asignación (= o bien <-)
  • Comparaciones (==, !=, <, <=)
  • Clases de objetos:
    • Vectores: por ejemplo c(1, 5, 6), c('edad', 'altura', 'grado'), 1:10
    • Matrices: por ejemplo matrix(data=1:10, ncol=2, byrow=TRUE)
    • Hojas de datos: por ejemplo data.frame(numeros=1:3, letras=c('a', 'b', 'c'))
    • Listas: por ejemplo list(comp1=1:5, comp2=matrix(1:9, 3), comp3="una lista")

1.2. Estructuras de flujo

  • Ciclos for: la variable contador recorre cada componente del vector vector y realiza el cuerpo del bucle. La variable contador es un nombre uevo mientras que el vector vector debe estar ya definido.
for(contador in vector) {
  # cuerpo del bucle
}
  • Condicional if: se pone una expresión lógica (normalmente una comparación o una función cuyo resultado es TRUE o FALSE) de modo que si se evalúa como TRUE, entonces se procesa el bloque entre llaves, y se evalúa como FALSE, entonces se salta el bloque entre llaves. Opcionalmente se puede complementar con un else y otro bloque a procesar, cuando la evaluación de la condición es FALSE
if(condicion) {
  # bloque a procesar si se cumple
  # la condición
}


if(condicion) {
  # bloque a procesar si se cumple
  # la condición
} else {
  # bloque a procesar si NO se cumple
  # la condición
}
  • Condicional while: se evalúa la condicionLogica y si resulta TRUE, entonces se procesa el cuerpo entre llaves, y se vuelve a evaluar la condicionLogica (indefinidamente hasta que resulte FALSE, en cuyo caso se salta el bloque entre llaves).
while( condicionLogica ) {
  # cuerpo del while
}

1.3. Manejo del error

  • Las funciones están preparadas para ciertos argumentos, de modo que no se puede operar con otro tipo de argumentos. Por ejemplo:
mean(c('a', 'b', 'c'))
## Warning in mean.default(c("a", "b", "c")): argument is not numeric or
## logical: returning NA
## [1] NA
aaa
## Error in eval(expr, envir, enclos): objeto 'aaa' no encontrado
  • Una forma de manejar el error, e informar al usuario, es realizar comprobaciones sobre los argumentos con if, y devolver algún mensaje por pantalla. Algunas comprobaciones sobre los datos son:
    • is.null(): devuelve TRUE si el argumento tiene el valor NULL (es decir, si no se ha inicializado).
    • is.na(): devuelve TRUE si el argumento tiene el valor NA.
    • is.numeric(): devuelve TRUE si el argumento está formado por constantes numéricas (vector o matriz).
    • is.matrix(): devuelve TRUE si el argumento es una matriz.
    • is.integer(): devuelve TRUE si el argumento es vector de números enteros.
    • Usar las comparaciones ==, >=, etc.
  • La forma de escapar de un script con un mensaje de error es:
    • stop('Mensaje de error')
  • Y la forma de avisar de algun problema (sin que llegue a ser un error, como, por ejemplo, usar un procedimiento que requiere de unas condiciones, que no se están cumpliendo) es:
    • warning('Mensaje de aviso')

1.4. Definir nuevas funciones en R

  • Las funciones se crean con la estructura:
nombreFuncion = function(argumento1, argumento2, ...) {
  # cuerpo de la función
  # return(nombreObjetoADevolver) opcional
}
  • Las características son:
    • argumento1, argumento2:
      • Son los argumentos de la función, que la función va a manejar en su cuerpo. Puede no haber argumentos
      • Se les puede dar un valor por defecto en la misma definición de la función, por ejemplo miFun = function(yo="Pablo", tu,...) {
    • ...: son optativos, y pueden recoger argumentos de otras funciones que se van a usar en el cuerpo de esta función
    • El cuerpo de la función contiene las instrucciones que queremos que realice cada vez que sea invocada.
    • return(nombreObjetoADevolver) se incluye cuando se desea que la función cree un objeto en memoria.

1.5. Cómo huir del bucle for: la función apply()

  • Los bucles for de R son extremadamente lentos.
  • Muchas funciones de R actúan componente a componente y no hace falta recorrer los vectores con un bucle for().
  • Ejemplo: calcular la suma de los invesos de los 1000 primeros números naturales
# MÉTODO 1:
suma = 0  # inicializa la suma a 0
for(i in 1:1000) {
  suma = suma + 1/i # añade cada sumando
}
suma

# MÉTODO 2: 
inversos = numeric(1000)  # inicializa el vector de inversos
for(i in 1:1000) {
  inversos[i] = 1/i # pone el inverso de cada uno
}
sum(inversos)

# MÉTODO 3: 
sum( 1/(1:1000) )  # la división ya es vectorial
  • La función apply() no es mucho más rápida que un bucle, pero ayuda a simplificar el código.
    • apply(X, MARGIN, FUN) aplica la función FUN sobre todos los vectores fila (si MARGIN=1) o columna (si MARGIN=2) de la matriz X.
  • Ejemplo: simular 1000000 de muestras de tamaño 100 del modelo uniforme en el intervalo [0,1], y calcular las medias de todas las muestras y almacenarlas en un vector llamado medias
# MÉTODO 1:
set.seed(123) # semilla para poder replicar la simulación 
medias = numeric(1e6) # inicializa vector de 1e6 ceros
for(i in 1:1e6) {
  v = runif(n=100, min=0, max=1) # simula muestra de tamaño 100
  medias[i] = mean(v) # su media va al vector de medias
}
summary(medias)

# MÉTODO 2: 
set.seed(123) # semilla para poder replicar la simulación
m = matrix(data=runif(n=100*1e6, min=0, max=1), # simula TODO
           ncol=1e6) # lo coloca en columnas de matriz
medias = apply(X=m, MARGIN=2, FUN="mean") # aplica "mean" a cada columna
summary(medias)

1.6. Funciones útiles para realizar gráficas

  • plot(x, y, ...) abre una ventana gráfica para realizar un gráfica. Un nuevo plot() machaca el anterior.
  • points(x, y, ...) añade puntos a una gráfica existente.
  • Argumentos útiles de la función plot():
    • x: coordenadas \(x\) de los puntos a representar
    • y: coordenadas \(y\) de los puntos a representar
    • type: tipo de gráfica.
      • "p": de puntos aislados
      • "l": los unes con línea
      • "h": como un histograma (verticales del eje X hasta cada punto)
      • "s": une con escalones (primero horizontal luego vertical)
      • "S": une con escalones (primero vertical luego horizontal)
      • "n": no dibuja
    • pch: forma del punto (número del 1 al 25)
    • col: color (número del 1 al 8, o cadena de texto con el nombre en inglés, ver ``)
    • lty: sólo si se elige type="l", es el tipo de línea (números del 1 al 6)
    • lwd: sólo si se elige type="l", es el grosor de la línea (números del 1 en adelante)
    • main: título principal (superior)
    • sub: título secundario (inferior)
    • xlab: etiqueta eje X
    • ylab: etiqueta eje Y
  • Ejemplo:
# DIBUJO DE LA NORMAL TIPIFICADA
x = seq(from=-4, to=4, by=0.01) # recorremos eje x
y = dnorm(x=x) # calculamos las f(x)
plot(x, y, type="l", main="Campana de Gauss", 
     xlab="x", ylab="Densidad de la normal", 
     col="red", lwd=3)
points(x=c(-4,4), y=c(0,0), type="l") # dibujamos eje X con 2 puntos
# vamos a lanzar 1000 puntos al azar en el recuadro
xx = runif(n=1000, min=-4, max=4) # sus coordenadas X
yy = runif(n=1000, min=0, max=1) # sus coordenadas Y
points(xx, yy) # los pintamos como puntos
dentro = yy < dnorm(xx) # condición "estar dentro de la gráfica"
color= rep("green", 1000) # vector de colores constante
color[dentro] = "red" # cambio a rojo si están dentro
points(xx, yy, col=color) # pinto de nuevo con color

  • Matrices de gráficas: la función layout() con una matriz como argumento, con entradas indicando dónde va cada “plot”, sirve para generar una ventana gráfica múltiple.
  • Ejemplo:
layout(matrix(data=c(1,1,1,2,3,3,2,4,5), ncol=3, byrow=TRUE))
stripchart(runif(100))
boxplot(rexp(100))
plot(1:10, sin(1:10), type='l')
hist(x=rnorm(50))
plot(1:10, 1/(1:10), type='l')

1.7. Otras funciones interesantes

  • cumsum(): devuelve las sumas parciales del vector argumento.
  • diff(): devuelve las diferencias entre cada par de componentes consecutivas del vector argumento
  • sign(): devuelve un vector con los signos (-1, 0 ó 1) de las componentes del vector argumento

2. EJERCICIOS BLOQUE PROGRAMACIÓN

  1. Define una función que calcule la suma de ciertos términos de la serie geométrica \(1 + r + r^2 + r^3 + \cdots\):
    • Argumentos:
      • r: la razón, que por defecto valga \(0.5\).
      • first: el primer exponente, que por defecto sea \(0\)
      • last: el último expoenente, que por defecto sea Inf (que corresponde a la serie completa).
    • Cuerpo:
      • Debe verificar que r es un valor numérico
      • Debe verificar que first es un valor entero no negativo
      • Debe verificar que last es un valor entero mayor o igual que first, o bien es Inf
      • Debe retornar el valor de la suma \(\texttt{r}^\texttt{first} + \texttt{r}^{\texttt{first}+1} + \cdots + \texttt{r}^\texttt{last}\)
      • Si last es Inf debes acudir a una fórmula que no involucre usar infinitos.
# Escribe aquí tu definición de función
# y compruebala para ver que cumple con todas
# las condiciones
  1. TEOREMA DEL LÍMITE CENTRAL: Simula 20 muestras de tamaño 300 del modelo uniforme en el intervalo [0,1]. Define una nueva muestra, media, también de tamaño 300, que sea la media aritmética término a término de las 20 muestras (es decir, cada valor de media es la media de los 20 datos que ocupan la misma posición en las 20 muestras). Representa una matriz de 3x3 histogramas, donde los 8 primeros histogramas correspondan a las 8 primeras muestras, y el último histograma sea de la variable media.
# Escribe aquí tu solución
  1. ENFOQUE FRECUENCIALISTA DE LA PROBABILIDAD: Simula un experimento que sea “éxito” con probabilidad 0.3. Haz 1000 simulaciones. Calcula un vector que calcule, a cada componente, la frecuencia relativa de “éxitos” acumulados hasta dicha componente. Representa una gráfica de dicho vector a lo largo del número de experimentos (eje X, números del 1 al 1000). Añade la recta \(y=0.3\) para comprobar que la frecuencia relativa se aproxima a dicha asíntota. Repetir 3 veces, representando con líneas de distinto tipo.
# Escribe aquí tu solución
# x = vector con las 1000 simulaciones (de 0s y 1s)
# frecRelAcum = vector con la frec. relat. acum. de 1s
# plot de x=1:1000, y=frecRelAcum (con tipo "línea")
# añadir recta y=0.3 con abline(h=0.3)
# repetir con un nuevo x y frecRelAcum
# añadir nueva gráfica con points()
# repetir con un nuevo x y frecRelAcum
# añadir nueva gráfica con points()
  1. CONTRASTE DE HIPÓTESIS SIMPLES: Define la función que resuelve el contraste de hipótesis simples \[ \left\{ \begin{array}{ll} H_0: & f_X = f_0 \\ H_1: & f_X = f_1 \end{array} \right. \] que minimiza una combinación lineal de las probabilidades de error tipo I y II. Si la combinación es \(a \cdot \alpha(\delta) + b \cdot \beta(\delta)\), y la muestra es \(\vec{x} = (x_1, x_2, \ldots, x_n)\), entonces el procedimiento óptimo rechaza \(H_0\) si y sólo si \(a \cdot f_0(\vec{x}) - b \cdot f_1(\vec{x}) < 0\).
# Escribe aquí tu solución
# función contrasteHipotSimples()
# argumentos los que dice el enunciado
# calcular el valor a*f0(muestra) - b*f1(muestra)
# crear una lista con el valor y otra con la frase
# aceptar H0 o rechazar H0 según el caso
# devolver la lista
  1. Escribe una función num2signs() que tome un vector de datos numéricos (argumento x) y devuelva un vector de “signos” (que puede ser numérico, por ejemplo, +1 y -1). Que el usuario elija (argumento type) en uno de los dos criterios posibles: (1) comparando cada dato con la mediana de todos los datos (type="median"), (2) comparando cada dato con el anterior (type="increment"). El usuario debe poder elegir el criterio, con el argumento type, cuyo valor por defecto sea "median". Las comparaciones que den lugar a 0 (ni +1 ni -1) se excluirán del vector de signos devuelto.
    • Ayuda: prueba la funcion diff() sobre un vector numérico.
    • Comprueba el funcionamiento con pequeños ejemplos.
# Escribe aquí tu solución
  1. Escribe una función nruns() que, para un vector de datos de dos categorías (argumento x, por ejemplo, de números +1 y -1), devuelva el número de rachas existentes en el vector. Recuerda que una racha es una cadena máxima de valores iguales, situados entre valores del otro tipo, o extremos de la muestra (ojo, una racha puede ser un solo signo). Si el vector tiene más de dos categorías, que devuelva un mensaje de error indicando que “el vector tiene más de dos signos!”.
    • Ayuda 1: la función unique() devuelve los valores “únicos” del vector argumento, sin repetirlos más que una vez.
    • Ayuda 2: puedes usar la función diff() sobre un vector cualquiera de signos +1 y -1, y comparar ambos vectores para ver cómo se pueden contar las rachas que hay (dificultad media-alta).
    • Comprueba el funcionamiento con pequeños ejemplos.
# Escribe aquí tu solución
  1. Define en R la función de probabilidad de la variable aleatoria \(R\) = “número de rachas”, para un proceso de dos signos, en el que se han observado \(n_1\) signos de un tipo y \(n_2\) signos del otro tipo (ver fórmula en teoría). La función tendrá la forma druns(x, n1, n2), donde n1 y n2 son los argumentos de la cantidad de signos de cada tipo y x puede ser cualquier valor posible de \(R\), de modo que druns() devuelve la probabilidad de cualquier número de rachas posible que se indique en x. Comprueba el funcionamiento con pequeños ejemplos.

La función de probabilidad del “número de rachas” es:

# Escribe aquí tu solución

3. CONTRASTES DE HIPÓTESIS NO PARAMÉTRICOS:

3.1. Contraste de bondad de ajuste de Pearson simple

  • Contrasta si una muestra viene de un modelo concreto (de variable discreta o continua, expresado como tabla de probabilidades) o no.
  • chisq.test(x,p): contraste sobre la bondad de ajuste de una muestra a una tabla de probabilidades
    • Argumentos:
      • x: vector con las frecuencias absolutas de cada categoría de la muestra
      • p: vector con las probabilidades de cada categoría de la muestra
    • Devuelve:
      • Valor del estadístico de contraste
      • \(p\)-valor
      • Opcionalmente, un warning si la aproximación no está garantizada

3.2. Contraste de bondad de ajuste de KOLMOGOROV-SMIRNOV

  • Contrasta si una muestra viene de un modelo concreto (de variable continua, expresado por la función \(F\) de distribución acumulada) o no.
  • ks.test(x, y)
    • Argumentos:
      • x: vector con datos de la muestra.
      • y: nombre (entrecomillado) de la “presunta” función \(F\) (según esté programada en R).
    • Devuelve:
      • Valor del estadístico de contraste
      • \(p\)-valor
      • Warning, si hay datos repetidos en la muestra

3.3. Contrasta de normalidad de LILLIEFORS

  • Contrasta si una muestra viene de “algún” modelo normal.
  • lillie.test(x) (instalar y/o cargar package ‘nortest’)
    • Argumentos:
      • x: vector con datos de la muestra.
    • Devuelve:
      • Valor del estadístico de contraste
      • \(p\)-valor

3.4. CONTRASTES DE ALEATORIEDAD

Asumimos que tenemos una variable aleatoria \(X\) de la que se ha extraido una muestra. Todos los datos se han obtenido de \(X\), pero no sabemos si esos datos son independientes entre sí.

Son contrastes para decidir entre:

  • creer que los datos de la muestra son independientes unos de otros, o bien
  • descartar esa creencia, porque parece que unos datos influyen sobre otros.

3.5. CONTRASTE DE RACHAS

Una racha es algo que se puede definir de varias formas: una serie de datos por encima o por debajo de la mediana, una serie de datos en ascenso o descenso, etc.

Cuando hay verdadera aleatoriedad, hay rachas de todo tipo, pero suele haber pocas rachas muy cortas y muy largas.

Los datos numéricos se tranforman en signos (según se comparan para formar las rachas), y cada cadena máxima de signos iguales es una racha, y su longitud es el número de signos que tiene.

  • MUESTRA: 1 variable, datos numéricos, que se transforman en signos
  • CONTRASTE: verificar si las rachas que se han formado son creíbles o no bajo la hipótesis de aleatoriedad (independencia)

\[ \left\{ \begin{array}{ll} H_0: & \text{datos aleatorios, independientes entre sí} \\ H_1: & \text{datos con cierta dependencia entre sí, no aleatorios} \end{array} \right. \]

  • FUNCIÓN: runs.test(x) del package ‘tseries’.
  • ARGUMENTOS:
    • x: vector con los “signos”.
  • DEVUELVE: una lista con valores, entre los que destaca el p.value, base de la decisión del contraste.
  • ¡ATENCIÓN!: usa aproximación normal incluso cuando no se debe (sólo conveniente para muestras con más de 20 signos de cada tipo)

3.6. CONTRASTE DE CUANTILES: signos

Si \(X\) expresa una variable aleatoria, su cuantil de orden \(p_0\) se representa por \(X_{p_0}\).

  • MUESTRA: 1 variable, datos numéricos
  • CONTRASTE: verificar si el valor \(x_0\) es creíble como cuantil del orden \(p_0\) para la muestra observada

\[ \left\{ \begin{array}{ll} H_0: & X_{p_0} = x_0 \\ H_1: & X_{p_0} (\neq, <, >) x_0 \end{array} \right. \]

  • FUNCIÓN: binom.test(x, n, p, alternative)
  • ARGUMENTOS:
    • x: número de signos + (que serían datos inferiores o iguales a \(x_0\)).
    • n: número total de datos.
    • p: orden del cuantil, es decir, \(p_0\), por defecto \(0.5\).
    • alternative: dirección de \(H_1\) (two.sided por defecto para \(\neq\), y atención a less para \(>\) o greater para \(<\)).
  • DEVUELVE: una lista con valores, entre los que destaca el p.value, base de la decisión del contraste.

4. EJERCICIOS BLOQUE CONTRASTES

Se deben contestar respondiendo a lo que se pregunta con palabras relativas al enunciado (prohibido mencionar \(H_0\) o \(H_1\) en las respuestas finales).

Carga el espacio de trabajo mt1021-1415-labo-s1-data.RData. En él están definidas ciertas variables para resolver los ejercicios

4.1. Simula 50 datos de un dado imperfecto definido por la tabla

  X        1    2    3    4    5    6
  f(X)  0.15 0.13 0.20 0.15 0.17 0.20

y luego contrasta si esa muestra es compatible con un dado perfecto, comentando el resultado, refiriendo al nivel de significación. Usa alguna semilla para la simulación.

# Escribe aquí tu solución

Escribe aquí tu comentario.

4.2. Averigua si los datos de la variable ‘x72’ son compatibles con el modelo de Poisson de media 2.5, comentando el resultado, refiriendo al nivel de significación.

# Escribe aquí tu solución

Escribe aquí tu comentario.

4.3. Simula 30 datos del modelo uniforme en el intervalo (0,1) y luego contrasta si esa muestra es compatible con el modelo normal de media 0.5 y varianza 0.25, comentando el resultado, refiriendo al nivel de significación. Usa alguna semilla para la simulación.

# Escribe aquí tu solución

Escribe aquí tu comentario.

4.4. Contrasta si los datos de la variable ‘x74’ son compatibles con el modelo normal de media 5 y varianza 1, comentando el resultado, refiriendo al nivel de significación.

# Escribe aquí tu solución

Escribe aquí tu comentario.

4.5. En caso de incompatibilidad en el ejercicio anterior, ¿se podría admitir al menos que los datos de la variable ‘x74’ son compatibles con el modelo normal? Comenta el resultado, refiriendo al nivel de significación.

# Escribe aquí tu solución

Escribe aquí tu comentario.

4.6. Sospechando sobre la falta de independencia entre las observaciones de una muestra, se comparan los datos con la mediana de los mismos, dando lugar a la cadena de signos de la variable ‘x76’. Realiza un contraste que arroje luz sobre este asunto. Comenta el resultado, refiriendo al nivel de significación.

# Escribe aquí tu solución

Escribe aquí tu comentario.

4.7. La mediana de la distribución de salarios en España se suponía de 650 EUR. Se sospecha que con la crisis ha disminuido. Se muestrea la población resultando los datos de la variable ‘x77’. Comenta el resultado, refiriendo al nivel de significación.

# Escribe aquí tu solución

Escribe aquí tu comentario.