Parcial II

Curso IMSER 2012

Instrucciones:

La carpeta comprimida “parcial-II.zip” tiene un script de R (“parcial.R”) en el cual usted deberá guardar todos los comandos del ejercicio, siguiendo la demarcación que se muestra en el archivo mismo.

La hoja de cálculo llamada 'datos.xls' contiene tres variables muestreadas. La primer variable son los tamaños corporales (en Kg) de 40 competidores del último Mundal de Magic, 20 hombres y 20 mujeres. La segunda variable es el sexo de los mismos, codificados como 1 para mujeres y 2 para varones. La tercer variable es la altura de los participantes, en metros.

Nota: los ejercicios del parcial son dependientes de los anteriores en el sentido de que utilizan objetos creados, pero no implica que no se puedan tratar de resolver independientemente.

Una vez terminado el parcial usted deberá comprimir la carpeta y subirla al EVA en la página correspondiente. Nótese que debe contener al menos dos archivos:

• “parcial.R”
  
• Un archivo de texto plano (.txt o .csv) con los datos de “datos.xls”, según lo indicado en la parte a.

1. Modelo Discreto para la funcion del Crecimiento Exponencial

La función exponencial busca simular un crecimiento exponencial tal como lo describiera Thomas R. Malthus, pero en su versión discreta (en el que la población, i.e. cantidad de individuos, en un tiempo t + 1 alcanza un valor equivalente al tamaño poblacional del tiempo anterior t multiplicado por una constante lambda). Para esto se eligió utilizar un loop for (nótese igualmente que existe una solución cerrada que permitiría, en R, evitar esto).

Considere ahora las distintas funciones candidatas y determine cuáles de ellas NO SON CORRECTAS (es decir, que tienen errores de diseño del loop interno).

Para ver las opciones vea el archivo .txt en este link.

Bajo condiciones muy poco realistas o de corta duración, las poblaciones naturales tienen el potencial de crecer exponencialmente en el medio donde viven. El crecimiento de una población esta determinado por el aumento de su número de individuos, el cual es constante en términos relativos y suele representarse por el parámetro \( \lambda \). Considerando los \( N_t \) como los tamaños poblacionales (número de individuos) para intervalos de tiempo \( t \) discretos, entonces lambda se obtiene así:

\[ \frac{N_t}{N_{t - 1}} = \lambda \; \text{(ec. 1.1)} \]

donde \( \lambda \) (lambda) es el parámetro poblacional de crecimiento exponencial discreto y su valor determina crecimiento (\( \lambda > 1 \)), decrecimiento (\( \lambda < 1 \)) o estásis (\( \lambda = 1 \)).

A partir de la ecuación 1.1 es posible encontrar una solución de \( N \) para todos los \( t \) posibles, a partir de un valor inicial \( N_0 \), asumiendo que la tasa de crecimiento \( \lambda \) es constante. Dicha solución es:

\[ N_t = N_0 \cdot \lambda ^ t \; \text{(ec. 1.2)} \]

En los siguientes ejercicios usted tomará los datos de la población humana global para estimar el parámetro de crecimiento anual y comparar gráficamente la predicción de este modelo con los datos reales.

1.a. Estimar los \( \lambda_t \)

El archivo poblacion.humana.rda contiene una data.frame con los datos de la población humana para todos los años entre 1950 y 2012. Con estos datos, calcule para cada par de años consecutivos el valor de \( \lambda \) usando la primer ecuación indicada. Nótese que se obtendrán un total de \( n - 1 \) valores de lambda, siendo \( n \) el total de datos (en este caso 63). Entonces para evitar confusiones, llamaremos de aquí en más \( \lambda_t \) al cociente \( N_t / N_{t - 1} \) (y no al cociente \( N_{t + 1} / N_t \)).

El vector de todos estos $\lambda$s se llamará lambda.

Nota: este ejercicio se puede realizar de varias maneras, utilice la que le resulte más cómoda.

1.b. (*) Función para calcular los \( \lambda_t \)

Aprovechando el código creado para el ejercicio anterior diseñe una función llamada calcLambda que tome como argumento datos de población (un vector numérico cualquiera) y devuelva un vector con las tasas de crecimiento para los datos iniciales, siguendo la misma convención indicada en el ejercicio anterior.

Recuerde que calcLambda debe funcionar correctamente para cualquier vector numérico que se use como argumento, sin importar su nombre, longitud o valor de los datos.

1.c. Estadísticos

A partir del vector lambda obtener los siguientes objetos:

• lprom: el valor promedio
• lmed: la mediana
• lq: vector con los percentiles 2.5 y 97.5 de la distribución de valores. Vea la función quantile para hacer la estimación.

Finalmente deberá guardar todos estos (4) valores en un vector llamado lest.

1.d. Modelo

Utilizando la equación 1.2 y los cuatro valores de lambda obtenidos construya los siguientes vectores:

• Nprom: tamaños poblacionales esperados para un crecimiento exponencial con una tasa de crecimiento constante igual a lprom.
• Nmed: tamaños poblacionales esperados para un crecimiento exponencial con una tasa de crecimiento constante igual a lmed.
• Nq1: tamaños poblacionales esperados para un crecimiento exponencial con una tasa de crecimiento constante igual a lq1.
• Nq2: tamaños poblacionales esperados para un crecimiento exponencial con una tasa de crecimiento constante igual a lq2.

Tome a consideración los siguientes detalles: (1) la población inicial (\( N_0 \)) debe ser la correspondiente a 1950, (2) los cuatro vectores tendrán la misma longitud y esta será igual a la cantidad de datos de población (63) (de forma que haya un valor distinto correspondiente a cada año de la serie original) y (3) el ejercicio se puede realizar tanto con y sin el uso de loops, aprovechando las capacidades de vectorización de R.

1.e. (*) Gráfico

Para este ejercicio se comenzará graficando los datos originales con el siguiente comando:

plot(poblacion.humana, xlab = "Año", ylab = "N", pch = 19, col = '#696969', cex = 1.5)

Posteriormente usted deberá agregar 4 curvas al gráfico utilizando un loop: las correspondientes a los vectores creados en el ejercicio anterior (1.d). El gráfico final será similar al siguiente:

Nótese que para agregar las líneas se utilizó el comando lines y que para hacerlo es necesario especificar en todos los casos los valores de ambos ejes (en el caso particular del plot inicial, al ser una data.frame de 2 columnas R no necesita que esto sea especificado).

Sugerencias: (1) juntar los cuatro vectores en una matriz o data.frame, (2) crear vectores con los valores de los parámetros que se quieran usar para las diferentes curvas, como col o lty (la elección de colores y tipos de línea son opcionales, el único requisito importante es que se pueda distinguir Nprom de Nmed).

2. Caminata del borracho

Como se mostró en las lecciones del curso, la “caminata del borracho” o “caminata aleatoria” es un proceso iterativo en el que una partícula se mueve aleatoriamente. En el caso más sencillo el movimiento es en una sóla dimensión, tal como se mostró en la lección mencionada.

2.a. (*) Loop

En este ejercicio usted deberá simular una caminata aleatoria similar, pero con ciertas diferencias. En este caso la partícula o “el borracho” se va a mover bajo las mismas reglas (empieza en el 0 y el movimiento es determinado por valores aleatorios tomados de una distribución normal \( N(0, 1) \)), pero en vez de detenerse cada vez que “choca” con los límites impuestos, simplemente va a “rebotar”. Entonces, las reglas seran:

1. La partícula empieza en la posición inicial 0.
2. El movimiento se determina de la siguiente manera: \( P_{t + 1} = P_t + Z_t \) en donde \( Z_t \) es una variable aleatoria con distribución \( N(0, 1) \).
3. Si la partícula sale de los límites -10 y +10, entonces es forzada a volver a entrar en la posición -9.5 o +9.5 respectivamente.
4. El número de iteraciones de la caminata es fijado de antemano en 500.

Hay total libertad de utilizar los nombres que se quiera para los distintos objetos.

La imagen que se muestra a continuación se obtuvo fijando los números aleatorios con set.seed(10) antes de correr la simulación y graficando con el comando:

plot(x, xlab = "Iteración", ylab = "Posición", type = "o", pch = 20, cex = 0.75)

2.b. (*) Función

Utilizando el código ya generado y modificándolo en donde sea pertinente, escriba una función que ejecute la simulación de caminata aleatoria llamada caminata. Dicha función aceptará los argumentos it y lim: el número de iteraciones y límite/borde de la caminata respectivamente. Nótese que ahora la partícula al salir de los límites deberá ser llevada a las posiciones -lim + 0.5 y lim - 0.5 según el caso.

El código y la imagen a continuación sirven de referencia para hacer el ejercicio:

# Corrida de la función 'caminata', con it = 250 y lim = 5:
set.seed(10)
plot(caminata(250, 5), xlab = "Iteración", ylab = "Posición", type = "o", 
    pch = 20)

3. Ómnibus

Para simular la cantidad de pasajeros de un ómnibus urbano se ha creado el código que aparece sobre el final del ejercicio (y que también se encuentra en el script parcial-II.R). El criterio es el siguiente: el bus recorre 25 paradas, empezando el trayecto sin pasajeros. En cada parada se subirán de 0 a 6 personas, pero debido a que existe un máximo estipulado de personas, a partir del momento en que se alcanza el total de 44 pasajeros el vehículo deja de parar.

(*) 3.a.

Completar el código: las líneas en blanco que se encuentran dentro de los límites del código indican en dónde debe cambiarse. El resto de las líneas están correctas.

3.b.

Modifique el código de la parte anterior de forma tal que haga lo mismo, pero utilizando un loop while.

(*) 3.c.

Modifique el código reparado en la parte 3.a para crear una función que ejecute la misma simulación, en la que el número de paradas y capacidad máxima del bus sean los argumentos de la misma.

3.d.

Hacer una variante del código (de cualquiera de las partes anteriores) en la que además de subir personas, a partir de la parada 10 se bajen entre 0 y 4 pasajeros por parada.

Código fuente:

paradas <- 25
pasajeros <- 0

registro[1] <- pasajeros
for (i in 1:paradas) {

  if (pasajeros >= 44) {

    registro[i:paradas] <- 44
    cat('Bus lleno!\n')
    break

  registro[i] <- pasajeros
  cat('Parada', i, 'hay', pasajeros, 'pasajeros\n')
}
plot(registro, xlab='Parada', ylab='No. de pasajeros')