Modulo 2: Escribiendo funciones en R

1
Modulo 2 Escribiendo funciones en R

• Por que queremos escribir funciones?
• Como escribimos funciones? Sintaxis basica
• Argumento de la funcion
• Salida de la funcion
• Ejemplos ( y lo que aprendimos de R al escribir
  esas funciones)
• Intervalo de confianza para la desviacion media
  poblacional.
• Intervalo de confianza para el cociente de dos
  desviaciones medias
• Calculo del periodograma
• Graficando el periodograma acumulado
• Estimando el espectro del proceso en base a las
  autocorrelaciones de la serie

2
Por que escribimos funciones?

• Las dos motivaciones principales son
  probablemente
• Tenemos que realizar un calculo laborioso muchas
  veces para diferentes conjuntos de datos. Es mas
  facil escribir una funcion y aplicarla luego que
  volver a realizar una vez y otra todos los
  calculos
• Algo que tenemos que calcular todavia no esta
  incorporado al software porque es de aparicion
  reciente. Generalmente toma un tiempo entre que
  una nueva formula aparece en un articulo y es
  incorporada al software comercial, a veces nunca
  es incorporada.
• Ejemplo recientemente una alumna en un programa
  de verano estaba comparando la variabilidad en el
  nivel de expresion de 100 genes de interes con la
  variabilidad de controles (ribozomas) en varios
  micro-arreglos. Minitab tiene el test de Levene
  para comparar varianzas pero yo queria que
  tambien compare las desviaciones medias, tenia la
  formula pero eso no esta todavia en ningun
  paquete, era mas rapido escribir la funcion y
  ensenarle a usarla que hacer los calculos
  requeridos paso por paso para los 24 experimentos
  que ella tenia.

3
Sintaxis basica para escribir funciones

• nombrelt-function(x,y)
• Comandos
• Salida

• Ustedes deciden el nombre de la funcion
• Dentro del parentesis (x,y) pueden escribir
  cuantos argumentos crean necesarios y darles el
  nombre que quieran
• Los comandos van separados por punto y coma
• Todo el cuerpo de la funcion va encerrado entre
  llaves
• El nombre que va al final (salida en el
  ejemplo) es el resultado que ustedes quieran que
  aparezca en pantalla. Los calculos intermedios no
  se imprimen
• Si queremos incluir comentario lo empezamos con
• Se acostumbra a poner al comienzo una
  descripcion de la funcion empezando esa linea con
  

4
Ejemplo 1 calculando un intervalo de confianza
para la desviacion media poblacional

• La desviacion media con respecto a la mediana
  (MAD) en una muestra se define como el promedio
  de las distancias de los datos con respecto a la
  mediana
• En el articulo Confidence Intervals for Mean
  Absolute Deviations de Bonett Seier (2003, The
  American Statistical Association, Vol 57 4) se
  dedujo el siguiente intervalo de confianza para
  la desviacion media poblacional

5
Intervalo de confianza para la desviacion media
poblacional

• La idea es construir un intervalo de confianza
  para el logaritmo de MAD y luego calcular el
  exponente de los extremos del intervalo.

• Podemos definir una funcion citau

6
Argumento de la funcion

• Para calcular la funcion citau tenemos que
  indicar a que datos (x) tenemos que aplicarla y
  con cuanta confianza queremos el intervalo (z
  seria el valor critico).
• Por ello citau es funcion de x y z.
• citault-function(x,z)

7
Los calculos a ejecutar

• mdmedian(x)
• mmean(x)
• vvar(x)
• taumean(abs(x-md))
• del(m-md)/tau
• nlength(x)
• cn/(n-1)
• gamv/(tau2)
• varlnt(delgam-1)/n
• setausqrt(varlnt)
• forlowlog(ctau)-zsetau
• foruplog(ctau)zsetau
• lowerendexp(forlow)
• upperendexp(forup)

• Utilizamos lo siguiente
• length de observaciones
• mean median
• median mediana
• Var varianza
• log (logaritmo natural)
• exp(a) ea
• sqrt (raiz cuadrada)
• / -

8
La salida

• Tenemos que crear un objeto con los resultados
  que querramos como salida. En el ejemplo del
  intervalo de confianza quisieramos el extremo
  inferior y superior del intervalo y el estimador
  puntual. Creamos un objeto que llamaremos ci
• cilt- c(lowerend, tau,upperend)
• Podriamos ponerle nombre a cada elemento de la
  salida introduciendo un vector de nombres como
• alt-c(extremo-inferior, estimacion-puntual,ext
  remo-superior)
• names(ci)a
• La ultima expresion de la funcion debe ser el
  nombre del objeto donde hemos almacenado el
  resultado que queremos que aparezca en pantalla,
  en este caso
• ci

9
Ejecutando la funcion

• Imaginemos que tenemos un conjunto de 100
  observaciones correspondientes a una muestra
  aleatoria de una poblacion. Los datos estan en el
  archivo WTAM1G. Primero leemos esos datos con
• Xlt-scan(awtam1g.dat)
• Para calcular el intervalo de 95 de confianza,
  z1.96
• Una vez que hemos copiado y pegado la funcion
  citau en R, simplemente escribimos
• citau(x,1.96) y obtenemos
• 0.6204765 0.7347000 0.8876145
• Si quisieramos el intervalo con 90 de confianza
  escribiriamos
• citau(x,1.645)

10
Ejemplo 2 intervalo de confianza para el
cociente de dos desviaciones medias poblacionales

• Dos poblaciones se podrian comparar en cuanto a
  dispersion en terminos de sus varianzas (test de
  Levene o de Barlett) o en terminos de sus
  desviaciones medias. En el articulo Confidence
  Intervals for Mean Absolute Deviations de Bonett
  Seier (2003, The American Statistical
  Association, Vol 57 4) se dedujo la siguiente
  formula para comparar dos poblaciones
  construyendo un intervalo de confianza para el
  cociente de las desviaciones medias
  poblacionales. Si el intervalo (con (1-alfa)100
  confianza) no cubre el valor 1 significa que
  rechazariamos la hipotesis nula de que las dos
  poblaciones tienen la misma desviacion media.

11

• La formula para el intervalo de confianza para el
  cociente de dos desviaciones medias es
• En la siguiente pagina (en dos columnas) esta la
  funcion rataus que implementa la formula. Las
  funciones y operaciones en R utilizadas son
  similares a las utilizadas en el ejemplo
  anterior. Recomendamos leerla relacionando las
  lineas del programa con la formula. En este caso
  tenemos dos conjuntos de datos x e y, por ello
  el argumento de la funcion es x,y,z (el valor
  critico).
• ratault-function(x,y,z)

12

• ratauslt-function(x,y,z)
• md1median(x)
• md2median(y)
• m1mean(x)
• m2mean(y)
• v1var(x)
• v2var(y)
• tau1mean(abs(x-md1))
• tau2mean(abs(y-md2))
• del1(m1-md1)/tau1
• del2(m2-md2)/tau2
• n1length(x)
• n2length(y)
• c1n1/(n1-1)
• c2n2/(n2-1)

• gam1v1/(tau12)
• gam2v2/(tau22)
• varlnt1(del1gam1-1)/n1
• varlnt2(del2gam2-1)/n2
• forratlog((c1tau1)/(c2tau2))
• seratsqrt(varlnt1varlnt2)
• forlowforrat-zserat
• forupforratzserat
• ratauexp(forrat)
• lowerendexp(forlow)
• upperendexp(forup)
• cilt- c(exinf,ratau,exsup) ci

13
Ejemplo 3. Graficando el periodograma de una
serie de tiempo.

• Escribimos la funcion perioplot(x) para calcular
  y graficar el periodograma de una serie de tiempo
  x. La transformada finita de Fourier fft es
  indispensable para hacer este calculo. En la
  pagina siguiente aparece el grafico obtenido
  aplicando la funcion perioplot a la serie de la
  temperatura del mar vista previamente.
• funcion para calcular el periodograma
• perioplotlt-function(x)
• adjxx-mean(x) substracts
  the mean of the series
• tffft(adjx) calculates
  finite Fourier transform
• nflength(tf) n2nf/21 decides
  the number of frequencies
• pritflt-tfc(1n2) takes the
  elements of the Fourier transform
• intensitylt-(abs(pritf2))/nf calculates
  the ordinates of periodogram
• nyquist1/2 pfreqlt-seq(0,nf/2,by1)
  preparation for frequencies
• freqlt-pfreq/(length(pfreq)-1)nyquist
  calculates frequencies
• plot(freq,intensity,type"l")

14
Continuacion de la pagina anterior (Ejemplo 3
graficando el periodograma)

• tempmarlt-scan(atempsea.dat)
• perioplot(tempmar)
• produce el siguiente grafico (primero tenemos que
  copiar y pegar la funcion perioplot de la pagina
  anterior en R)

15
Ejemplo 4. Uso de la funcion cusum de R para el
periodogram acumulado.

• Esta funcion se parece mucho a la anterior pero
  se utiliza tambien la funcion cumsum que produce
  la suma acumulada de los elementos de un vector.
• periocumlt-function(x)
• adjxx-mean(x)
• tffft(adjx) nflength(tf) n2nf/21
• pritflt-tfc(1n2)
• intensitylt-(abs(pritf2))/nf
• nyquist1/2 pfreqlt-seq(0,nf/2,by1)
• freqlt-pfreq/(length(pfreq)-1)nyquist
• cumintlt-cumsum(intensity)/(max(cumsum(intensity)))
  
• plot(freq,cumint,type"l")

16
Continuacion del Ejemplo 4 el periodograma
acumulado.

• Asumiendo que ya hemos leido anteriormente los
  datos de la serie de temperatura del mar,
• escribimos
• periocum(tempmar) para obtener el grafico.

17
Continuacion del Ejemplo 4

• Escribimos la funcion periocum de las paginas
  anteriores con fines didacticos pero R tiene ya
  incorporada una funcion para calcular el
  periodograma acumulado
• Si escribimos
• cpgram(tempmar) obtenemos el grafico de la
  derecha. Ese grafico incluye las bandas para el
  test de ruido blanco.

18
Ejemplo 5. Estimando el espectro de una series de
tiempo

• La idea de incluir esta funcion es mostrar como
  en R puede usarse operaciones con matrices en
  lugar de los ciclos iterativos (loops) del tipo
  do..while etc. que usamos en otros lenguajes.
  Por lo general se recomienda evitar los loops.
• La formula que queremos calcular en este caso es
  

19
Notar que estamos utilizando dos operadores que
no habiamos visto antes outer y . El lector
puede usar el help en R para ver detalles
sobre ellos.

• Esta funcion estima el espectro de una serie y
  tiene como entrada las autocorrelaciones rj de
  la serie que han sido puestas en au. Como
  dijimos el proposito de incluirla en estas notas
  es enfatizar que una suma se puede calcular en
  base a operaciones de matrices en lugar de un
  loop o ciclo iterativo.
• estspeclt-function(au)
• Mlt-length(au) counts how many
  autocorrelations (M) we read
• jlt-seq(1,M,by1) creates subindeces j1M
• lam0.5(1cos(jpi/M)) calculates Tukeys
  weights
• wlt-seq(0,pi,bypi/50) calculates angular
  frequencies
• lact(lam)au multiplies each weight by the
  corresponding autocorrelation
• flt-function(j,w) cos(jw)
• zlt-outer(j,w,f)
  calculates cos j w for all values of w and j
• szlt-lacz obtains the
  sum of weights correlations cos jw
• hlt-(1/(2pi))(12sz) calculates
  h(w)
• plot(w,h ,type"l",mainEstimated spectral
  density )

20
Ejemplo 5 -continuacion

• Leemos las primeras 30 auto- correlaciones de la
  temperatura del mar
• ault-c(0.827,0.537,0.241,0.003,-0.155,-0.223,
  -0.196, -0.087,0.075,0.253,,0.379,0.402,0.306,0.12
  4,-0.084,-0.256,-0.380,-0.429,-0.382,-0.247,-0.055
  ,0.141,0.288,0.352,0.289,0.123,-0.076,-0.261,-0.39
  0,-0.444)
• Luego escribimos
• estspec(au)
• para obtener el grafico de la derecha
• Nota R tiene incorporada la funcion spec.pgram
  para estimar el espectro utilizando un metodo
  diferente de estimacion.

21
Ahora practiquemos el escribir y ejecutar
funciones.

• Ejercicio 1. a) El resumen de 5 numeros de un
  conjunto de datos esta compuesto por el minimo,
  primer cuartil,mediana, tercer cuartil y maximo.
• Podemos usar funciones existentes en R para
  escribir nuevas funciones y podemos ponerles el
  nombre que deseamos.
• Usar las funciones max(), min(), median(), y
  quantile( , ) para escribir una nueva funcion que
  calcule el resumen de 5 numeros del conjunto de
  datos x. Llamar a la nueva funcion cinco.
• b) Para comprobar si su funcion esta bien
  escrita, ingresar los siguientes datos
• xlt-c(0,2,9,1,8,7,5,5,6,2,5,1,9,4,8)
• Ejecutar la funcion escribiendo
• cinco(x)
• El resultado debe ser 0.0 2.0 5.0 7.5 9.0
• c)Ya existen funciones en R que calculan el
  resumen de los cinco numeros quantile(x) y
  summary(x) (la ultima tambien da la media).
  Ejecutar esas funciones y comparar los resultados
  con los de la funcion cinco(x).
• Que haria Ud para que la funcion 5 incluya los
  nombres de los resultados tal como lo hacen las
  funciones quantile y summary?

22

• Ejercicio 2.
• Escribir una funcion histolog que calcule el
  logaritmo natural de los datos en x y luego
  grafique el histograma de los datos transformados
  (la salida de la funcion debe ser el histograma).
  
• Ejercicio 3.
• Escribir una funcion (Ud decide el nombre) que
  para observaciones a una variable x, arroje como
  salida el numero de observaciones, la media, la
  mediana y el valor absoluto de la diferencia
  entre media y mediana. Aplique la funcion a un
  conjunto de datos.
• Ejercicio 4.
• Escribir una funcion (Ud decide el nombre) que
  para observaciones a dos variables x e y , arroje
  como salida la media de cada variable y la
  correlacion entre las dos.
• Aplique la funcion a observaciones de dos
  variables.