Ejercicios de TIplE

Índice

1 Días 1 y 2: Ubuntu

  1. Crea una carpeta llamada informes dentro de tu carpeta personal. 
    cd
mkdir informes
  2. Mediante el procesador de texto de LibreOffice, crea un documento llamado informe-ubuntu (LibreOffice le añadirá la extensión .odt por omisión) y guárdalo en la carpeta informes.

    El documento debe contener la ecuación: \(\Gamma(\alpha) = \psi_\beta.\Delta^\omega\)

    Prueba a incluir también

    • fórmulas más complejas
    • una tabla
    • una lista numerada
    • una lista no numerada
    • un título
    • una cabecera y un pie de página
    • un texto en cursiva, otro en negrita y otro monoespaciado
    • diversos colores de texto y colores de fondo
  3. Prueba una sustitución global (por ejemplo, cambia todas las comas por punto y coma) mediante Edit / Find and Replace.
  4. Exporta el documento a PDF (se creará un documento llamado informe-ubuntu.pdf en la carpeta informes).
  5. Crea una carpeta llamada « .secreta » dentro de tu carpeta personal. 
    mkdir ~/.secreta
  6. Mueve el fichero editable informe-ubuntu.odt a .secreta.
  7. Localiza en el buscador el editor de imágenes GIMP. Arrastra su icono hasta fijarlo en el panel de la izquierda (lanzador). Luego, pulsa sobre el icono con el botón derecho y elige «No mantener en el lanzador» para eliminar el icono.
  8. Descarga el escudo de la universidad de la dirección http://www.uniovi.es/documents/31582/23504333/Logo+Universidad+de+Oviedo+centrado.png/c976d431-0532-4725-919f-da660abba00f?t=1468927003744 (o busca «escudo uniovi» y elige uno en formato PNG).
  9. Abre dicho escudo con Gimp y modifícalo (por ejemplo, con la herramienta de selección rectangular selecciona alguna parte y elimínala con Editar / Cortar; o ve a Imagen / Autorecortar imagen).
  10. Guarda el fichero en la carpeta informes. Se guardará con la extensión .xcf que guarda todos los elementos del proceso de edición de la imagen (capas, textos…).
  11. Exporta el fichero a la carpeta .secreta, con formato .jpg. Abre desde Nautilus, con el Visor de Imágenes, el nuevo fichero creado. Puedes comprobar que ya no es trasparente y (si no escogiste calidad 100%) que, al ampliarlo, se ven algunos defectos).
  12. En el documento informe-ubuntu.odt inserta el escudo de la universidad en la cabecera, con un tamaño más o menos adecuado.
  13. Vuelve a exportar el documento a formato PDF a la carpeta informes.
  14. Carga el archivo paises.ods de http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos y explora las definiciones de referencias relativas y absolutas en http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/apuntes/libreoffice-calc/practica1-217.odt.
  15. Calcula el total mundial de población. 
    5.189.503
  16. Calcula la extensión de cada país en quilómetros cuadrados. 
    extensión <- 1000*población / densidad
  17. Abre una nueva hoja en paises.ods y crea la siguiente tabla:
    tasa de I.V.A. = 16 %
    PRECIO sin I.V.A. PRECIO con I.V.A.  
    200    
    115    
    2200    
    30    
  18. Crea fórmulas para calcular los precios con IVA, haciendo uso del valor almacenado en la celda con el 16. Cambia la tasa de 16 a 21, de forma que todos los precios con IVA se actualicen automáticamente.
  19. Calcula el P.I.B. percápita por país y el global. 
    P.P.C. <- P.I.B.*1e8 / (población*1000) # $ EE.UU. / habitante
PPCglobal <- sum(P.I.B.*1e8) / sum(población*1000)
  20. Calcula el consumo de calorías medio global. 
    weighted.mean (calorías, población, na.rm=TRUE) # discutible

2 Días 3 al 7: GNU bash y GNU coreutils

  1. Localiza dónde se encuentra el fichero llamado «sources.list» y muestra su contenido. 
    locate sources.list                # rápido pero si es demasiado reciente no lo encuentra
find / -name sources.list          # lento pero seguro
  2. Busca en carleos2.epv.uniovi.es todos los ficheros .pdf del directorio /tmp/ROK 
    find /tmp/ROK -name '*.pdf'
  3. Crea en tu carpeta personal la siguiente estructura de directorios: 
                  +-----+         
              |  ~  |
              +-----+
                | | 
        +-------+ +------+
        |                |
        v                v
   +---------+      +--------+
   |Descargas|      |Imágenes|
   +---------+      +--------+
     |     |          |     |
     V     V          V     V
+------+ +-----+    +---+ +---+
|libros| |pelis|    |JPG| |PNG|        
+------+ +-----+    +---+ +---+

    mkdir ~/Descargas
mkdir ~/Descargas/libros
mkdir ~/Descargas/pelis
mkdir ~/Imágenes/JPG
mkdir ~/Imágenes/PNG
# o todo a la vez:
mkdir -p ~/Descargas/{libros,pelis} ~/Imágenes/{JPG,PNG}
  4. Consigue una serie de imágenes en formato JPG (por ejemplo, las de la final en Badajoz de la Incubadora de Sondeos y Experimentos) y colócalas en la carpeta «JPG». 
    cd ~/Imágenes/JPG
wget -r --no-parent http://bellman.ciencias.uniovi.es/incubadora/edicion4/fotos/badajoz-asturianos/
find bellman.ciencias.uniovi.es -iname '*.jpg' -exec cp {} . \;
rm -rf bellman.ciencias.uniovi.es
  5. Mediante un bucle «for», conviértelas al formato PNG, de forma que los ficheros PNG resultantes acaben en la carpeta «PNG». La orden para convertir imagen.jpg en imagen.png es convert imagen.jpg imagen.png. 
    for i in *.JPG; do mv $i $(basename $i .JPG).jpg; done 
for i in *.jpg; do convert $i ../PNG/$(basename $i .jpg).png; done
## otra forma de hacerlo, sin bucles "for":
rename 'y/\.JPG$/.jpg/' *.JPG
cp *.jpg ../PNG
cd ../PNG
mogrify -format png *.jpg
  6. La orden wc permite contar renglones, palabras y caracteres de un fichero de texto. En la carpeta datos tienes un fichero llamado vacas.csv. Cuenta cuántos registros contiene ayudándote de wc. Luego carga el fichero en Libreoffice Calc y comprueba el resultado. 
    wget http://bellman.ciencias.uniovi.es/~carleos/manadine/datos/vacas.csv
wc -l vacas.csv
## son 50.000 registros más un renglón de cabecera
  7. En la misma carpeta está el fichero pisa2012.csv.gz. (También aquí.) El sufijo .gz indica que está comprimido. Puedes descomprimirlo con el comando gunzip pisa2012.csv.gz, con lo que se convierte en un fichero de texto normal y corriente.
    • ¿Cuántos registros contiene?
    • ¿Cuántas variables? Ten en cuenta las órdenes head y tr; pueden ser útiles (considera las opciones -n 1 para head y -c y -d para tr). 
    wget http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/pisa2012.csv.gz
gunzip pisa2012.csv.gz
wc -l
## 271324 renglones (un registro menos)
head -n1 pisa2012.csv | tr -cd ',' | wc -c
## da 614 comas, luego hay 615 campos o variables

    Si no quieres descomprimir un fichero físicamente, sino sólo usarlo al vuelo, puedes usar zcat. Por ejemplo: 

    zcat pisa2012.csv.gz | wc

    p <- read.csv (gzfile ("pisa2012.csv.gz"))
nrow (p)
ncol (p)
dim (p)
  8. Compara los resultados de las siguientes órdenes e intenta averiguar la razón: 
    time zcat pisa2012.csv.gz | wc
time zcat pisa2012.csv.gz | wc -l

    # wc -l no pierde el tiempo contando palabras ni caracteres
  9. Intenta averiguar qué hace esta orden: 
    zcat pisa2012.csv.gz | cut -d',' -f1 | uniq -c | sort -n

    ## "cut" se queda con el primer campo, CNT = país
## "uniq" cuenta cuántos registros hay de cada país
## "sort" ordena para que los más numerosos estén al final

    p <- read.csv (gzfile ("pisa2012.csv.gz"))
table (p$CNT)
sort (summary (p$CNT))
  10. Ídem: 
    zcat piso2012.csv.gz | head -n 1 > /tmp/muestra.csv
zcat piso2012.csv.gz | tail -n +2 | shuf -n 100 >> /tmp/muestra.csv

    ## crea un fichero con cien registros extraídos al azar del principal;
## los cien registros están precedidos por una cabecera con los nombres de los campos;

    p1 <- p [sample (nrow(p), 100), ]
  11. Obtén un fichero /tmp/muestra-es.csv que contenga los registros de PISA relacionados con España (CNT=ESP). 
    : zcat piso2012.csv.gz | head -n 1 > /tmp/muestra.csv
: zcat piso2012.csv.gz | grep ^ESP >> /tmp/muestra.csv

    p2 <- p [p$CNT == "ESP", ]
p2 <- p [which (p$CNT == "ESP"), ] # por si hubiera algu'n CNT ausente
  12. Averigua qué hace esta instrucción explorando cada tubo: 
    cat /usr/share/common-licenses/GPL|tr ' ' '\n'|sort|uniq -c|sort -n
    • Mete en el tubo una instrucción para convertir mayúsculas a minúsculas. 
      # por ejemplo: tr '[A-Z]' '[a-z]'
cat /usr/share/common-licenses/GPL|tr ' ' '\n'|tr '[A-Z]' '[a-z]'|sort|uniq -c|sort -n
    • Ídem para eliminar todos los símbolos que no sean letras minúsculas. 
      # por ejemplo: tr -cd '[a-z]\n'
cat /usr/share/common-licenses/GPL|tr ' ' '\n'|tr '[A-Z]' '[a-z]'|tr -cd '[a-z]\n'|sort|uniq -c|sort -n
    • Intenta hacer lo mismo con los ficheros texto1 … texto4 sitos en ftp://bellman.ciencias.uniovi.es. Para quitar acentos, tendrás que meter en el tubo la orden iconv ó recode. 
      # con "iconv" o con "recode":
cat texto*|tr ' ' '\n'|tr '[A-Z]' '[a-z]'|tr -cd '[a-z]\n'|iconv -t ascii//translit|sort|uniq -c|sort -n
cat texto*|tr ' ' '\n'|tr '[A-Z]' '[a-z]'|tr -cd '[a-z]\n'|recode ..flat|sort|uniq -c|sort -n
  13. Obtengamos el código fuente de un paquete de Ubuntu, por ejemplo así: 
    apt-get source r-cran-rcmdr

    Metámonos en el directorio con las fuentes: 

    ls -lrt
cd rcmdr-2.0-3

    ¿Cuántos ficheros hay? 

    find . | wc                              # juntos ficheros y directorios;
find . -ls | cut -c 13 | sort | uniq -c  # por separado;

    ¿Cuántos son de código fuente R? 

    ## Dos formas:
find . -name '*.R' | wc
find . | grep '\.R$' | wc

    ¿Cuántos renglones supone el código fuente de R? 

    ## Dos formas:
find . -name '*.R' -exec cat {} \; | wc
find . -name '*.R' | xargs cat | wc

    ¿Cuántas definiciones de funciones de R hay? Una definición debe contener algo como « <- function ( » 

    find . -name '*.R' | xargs cat | grep '<-.*function.*(' | wc
  14. Obtén las frecuencias de los sexos en el fichero vacas.csv, en la terminal y en R 
    cut -d';' -f10 vacas.csv|sort|uniq -c

    summary(read.csv2("vacas.csv")$sexo)
  15. Obtén los municipios más frecuentes en el fichero vacas.csv, en la terminal y en R. 
    cut -d';' -f1 vacas.csv|sort|uniq -c|sort -n|tail
R -e 'head(summary(read.csv2("vacas.csv")$nomMunicipio))'
  16. Halla en qué se diferencian los ficheros vacas.csv y vacas0.csv mediante el comando diff (pista: se ha usado 
    cat vacas.csv | sed 's/;;/;0;/g' | sed 's/;;/;0;/g' > vacas0.csv

    para crear el fichero vacas0.csv; ¿por qué se repetirá el sed?).

  17. A partir de la información de esta página y este directorio intenta averiguar qué hace esta orden: 
    for i in $(seq 1987 2008); do wget http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/vuelos/$i.csv.bz2; done
  18. En la misma carpeta que vacas.csv está el fichero vuelos.csv.gz. Su contenido está descrito en http://stat-computing.org/dataexpo/2009/the-data.html
    • ¿Cuántos registros contiene? (intenta hacerlo lo más rápido posible) 
    zcat vuelos.csv.gz | wc -l  # salen 123.534.970 renglones
    • ¿Cuántas variables? 
    zcat vuelos.csv.gz | head -n1 | tr -cd ',' | wc -c # sale 28, así que 29 campos
    • ¿Cuántos registros se refieren al aeropuerto internacional de la capital de Iowa? (Pista: halla el código de su aeropuerto.) 
    zcat vuelos.csv.gz | grep -w DSM | wc -l # salen 446.224
  19. Visita la página del I.N.E.
    • Pincha a la izquierda en el menú de navegación.
    • Pincha en "Productos y servicios".
    • Pincha en "Información estadística".
    • Pincha en "Ficheros de microdatos".
    • Pincha en "Censos de población".
    • Pincha en "Censos 2011".
    • Pincha en "Ficheros de microdatos: Nacional".
    • Guarda el archivo y descomprímelo.
    • En el apartado "Personas y hogares", pincha en "diseño de registro y valores válidos de las variables".
    • ¿Cuántas personas de cada sexo aparecen registradas? 
    cat MicrodatosCP_NV_per_nacional_3VAR.txt | cut -c 43 | sort | uniq -c

3 Días 8 al 22: R

  1. Ejecuta R en un terminal, crea una variable llamada altura que valga 170 y sal de R guardando el espacio de trabajo. Comprueba la naturaleza de los nuevos ficheros que aparecen en el directorio actual: 
    ls -la                     # deben aparecer .Rhistory y .Rdata;
file .*                    # .Rhistory debe identificarse como texto y .RData como un fichero comprimido;
less .Rhistory             # deben aparecer los comandos realizados;
zcat .RData | less         # no se ve gran cosa; se trata de un fichero binario;
zcat .RData | strings      # las cadenas contenidas en el fichero;
zcat .RData | hexdump -C   # el contenido del fichero octeto por octeto;
  2. Carga en R el fichero vacas.csv mediante el comando read.csv2 o mediante Rcmdr > Datos > Importar > Texto. Llamaremos v a la variable que contiene los datos. 
    v <- read.csv2 ("vacas.csv", na.strings="0")
  3. Representa la variable peso mediante un histograma. 
    hist(v$peso)
  4. Representa el peso en función del sexo mediante cajas. 
    boxplot (peso ~ sexo, v)
  5. Halla la diferencia en media y mediana entre sexos. 
    summary (v$peso [v$sexo == "H"])
summary (v$peso [v$sexo == "M"])
by (v$peso, v$sexo, summary)
  6. Representa gráficamente las frecuencias de los municipios. 
    pie     (table (v$nomMunicipio))   # tarta
barplot (table (v$nomMunicipio))   # barras
  7. (Ejecución por lotes) Prueba a usar R en un tubo de terminal: 
    echo '1-2*(1-pnorm(1.96))' | R --vanilla
  8. Crea mediante un editor de textos un pequeño programa en R. Por ejemplo: 
    alfa    <- 0.05
cuantil <- qnorm (1 - alfa/2)
cat ("El cuantil gausiano bilateral asociado a alfa=", alfa, "es", cuantil, ".\n")

    Llama al fichero ejemplo.R. Entonces puedes ejecutarlo de alguna de estas maneras: 

    cat ejemplo.R | R --vanilla
R --vanilla < ejemplo.R          # equivale a la anterior;
R -e 'source("ejemplo.R")'       # no hace eco de cada orden;
R CMD BATCH ejemplo.R            # la salida se guarda en ejemplo.Rout
  9. Sea x el número de incidencias semanales registrados en cierto control. Los valores de x obtenidos hasta ahora son: x <- c (3, 7, 5, 9, 12, 3, 14, 2).
    1. ¿Cómo seleccionar el primer valor de x? 
      x[1]
    2. ¿Cómo seleccionar el tercero? 
      x[3]
    3. ¿Cómo obtener un vector que contenga los dos primeros valores? 
      x[1:2]
    4. ¿Cómo obtener un vector que contenga todos los valores menos el tercero? 
      x[-3]
    5. ¿Qué posiciones ocupan los valores menores que 10? (Pista: which). 
      which (x < 10)
    6. ¿Cómo obtener un vector que contenga los valores menores que diez? 
      x[x<10]
x[which(x<10)]  # funciona aunque x contenga NA
    7. Sustituye los elementos de x menores que diez por un cero. 
      x[x<10] <- 0
  10. Crea un vector llamado precipitacion en R y que contenga los siguientes datos (pulgadas anuales de precipitación sobre Mineápolis - San Pablo):
    0'77 1'74 0'81 1'20 1'95 1'20
    0'47 1'43 3'37 2'20 3'00 3'09
    1'51 2'10 0'52 1'62 1'31 0'32
    0'59 0'81 2'81 1'87 1'18 1'35
    4'75 2'48 0'96 1'89 0'90 2'05 
    precipitacion <- c (0.77 , 1.74 , 0.81 , 1.20 , 1.95 , 1.20 ,
                    0.47 , 1.43 , 3.37 , 2.20 , 3.00 , 3.09 ,
                    1.51 , 2.10 , 0.52 , 1.62 , 1.31 , 0.32 ,
                    0.59 , 0.81 , 2.81 , 1.87 , 1.18 , 1.35 ,
                    4.75 , 2.48 , 0.96 , 1.89 , 0.90 , 2.05 )
## puedes obtener lo anterior usando un editor de texto y remplazando...
## ...apóstrofos por puntos
## ...tabuladores por comas
## otra opción es
## 1. escribir en la consola
##     precipitacion <- scan(dec="'")
## 2. pulsar Intro
## 3. pegar en la consola el texto con los datos seleccionados
## 4. pulsar Intro dos veces
  11. Trasforma precipitacion para que esté expresado en milímetros (litros por metro cuadrado). 
    precipitacion <- precipitacion * 25.4
  12. Halla la media, mediana, desviación típica y recorrido intercuartílico de precipitacion. 
    mean (precipitacion)
median (precipitacion)
sd (precipitacion)
IQR (precipitacion)
  13. Ídem de los veinte últimos años. 
    mean (precipitacion[11:30]) # etc.
  14. Ídem de todos los años salvo el antepenúltimo. 
    mean (precipitacion[-28])
  15. Ídem de todos los años secos (un año es seco si llueve menos de cincuenta milímetros). 
    mean (precipitacion[precipitacion<50])
  16. Halla la media, mediana, desviación típica y recorrido intercuartílico de la variable número de hijos:
    número de hijos 0 1 2 3 4 5 6 7
    número de familias 259 486 329 49 40 3 0 1

    (Aparte: puedes usar scan() para importar en R cómodamente datos separados por espacios (por ejemplo, al copiarlos desde HTML ó PDF y pegarlos en el terminal.) 

    weighted.mean (nh, nf)  # una posibilidad, sólo para la media;
muestra <- rep (nh, nf) # otra posibilidad: creamos las muestra completa;
mean (muestra)          # etc.
  17. Crea un dataframe que contenga sólo las variables número de hijos y número de coches, pero que represente a la muestra completa:
    número de hijos 0 0 1 1 1 2 2 2
    número de coches 0 1 2 0 1 0 1 2
    número de familias 18 31 14 3 10 9 5 21 
    nhijos <- c(0,0,1,1,1,2,2,2); ncoches <- c(0,1,2,0,1,0,1,2); nfam <- c(18,31,14,3,10,9,5,21)
## opción similar a la del ejercicio anterior
nhijos.muestra  <- rep (nhijos,  nfam)
ncoches.muestra <- rep (ncoches, nfam)
datos <- data.frame (num.hijos = nhijos.muestra, num.coches = noches.muestra)
## otra versión más elaborada y generalizable
## (sirve para casos con el dataframe ya construido)
datos.aux <- data.frame (nhijos, ncoches, nfam)
datos <- datos.aux [rep (1:nrow(datos.aux), datos.aux$nfam)), -3]
  18. En un experimento se examinó el efecto de la densidad plantar sobre la producción de tomate. Se plantaron tres variedades de tomate a cuatro diferentes densidades. Cada combinación se replicó tres veces. He aquí los resultados:
    Densidad Variedad 1 Variedad 2 Variedad 3
    1 9'2 12'4 5'0 8'9 9'2 6'0 16'3 15'2 9'4
    2 12'4 14'5 8'6 12'7 14'0 12'3 18'2 18'0 16'9
    3 12'9 16'4 12'1 14'6 16'0 14'7 20'8 20'6 18'7
    4 10'9 14'3 9'2 12'6 13'0 13'0 18'3 16'0 13'0

    Creamos la tabla de producciones así: 

    prod <- c ( 9.2, 12.4,  5.0,  8.9,  9.2,  6.0, 16.3, 15.2,  9.4,
           12.4, 14.5,  8.6, 12.7, 14.0, 12.3, 18.2, 18.0, 16.9, 
           12.9, 16.4, 12.1, 14.6, 16.0, 14.7, 20.8, 20.6, 18.7, 
           10.9, 14.3,  9.2, 12.6, 13.0, 13.0, 18.3, 16.0, 13.0)

    ¿Cómo usarías rep para crear vectores que representen la variedad y la densidad?

    (Pista: rep(1:nd,rep(nr*nv,nd)), rep(1:nd,each=nr*nv); rep(rep(1:nv,rep(nr,nv)),nd), rep(1:nv,each=nr,nd)) 

    densidad <- rep (1:4, each=9)
variedad <- rep (1:3, each=3, 4)
  19. Halla medias y desviaciones típicas de las prod asociadas a la variedad 1. 
    mean (prod[variedad==1])
  20. Halla medianas y recorridos intercuartílicos de las prod asociadas a la variedad 2 y a la densidad 3. 
    median (prod [variedad==2 & densidad==3])
  21. Crea una matriz con los datos de prod, densidad y variedad, por columnas. 
    m <- cbind (prod, densidad, variedad)
  22. De la anterior matriz, crea una submatriz que contenga sólo las filas asociadas a la variedad 2 y a la densidad 3. 
    m1 <- m [variedad==2 & densidad==3, ]
  23. De la anterior matriz, crea una submatriz columna que contenga sólo la densidad. (Pista: explora la opción drop de [, o bien explora la orden matrix.) 
    m [, "densidad", drop=FALSE]
  24. Haz un histograma de precipitacion. 
    hist (precipitacion, col=2)
  25. Mete el histograma en una variable y analiza sus componentes mediante names, str, ls.str, as.list… 
    h <- hist (precipitacion, col=2)
names (h)
str (h)
ls.str (h)
as.list (h)
  26. Haz un gráfico de cajas de precipitacion. 
    boxplot (precipitacion)
  27. Haz un gráfico que represente la secuencia temporal de precipitacion. 
    plot (precipitacion, type="l")
  28. Haz un gráfico de barras de x <- rbinom (200, 10, .7). Que aparezca la frecuencia sobre cada barra. 
    frecuenc <- table (x)
abscisas <- barplot (frecuenc)
text (abscisas, frecuenc + 0.3, frecuenc, xpd = TRUE)
## El 0.3 es arbitrario y puede ser necesario retocarlo.
  29. Haz un gráfico de la función \(f(x)=2x+x^2\) entre \(-2\) y \(3\). 
    plot (x <- seq(-2,3,.01), 2*x+x^2, type="l")
  30. Representa prod en función de densidad con un gráfico de dispersión (plot) y un gráfico de cajas (boxplot). 
    plot (prod ~ densidad)
boxplot (prod ~ densidad)
  31. Realiza una regresión simple (lm) de la prod en función de densidad. 
    summary (lm (prod ~ densidad))
  32. Representa gráficamente (abline) la recta de regresión junto a los puntos. 
    plot (densidad, prod)
abline (lm (prod ~ densidad), col=2)
  33. Genera un vector que contenga el \(R^2\) múltiple de la regresión y el p-valor asociado a la ordenada del origen. 
    l <- summary (lm (prod ~ densidad))
v <- c (l$r.squared, l$coefficients["(Intercept)","Pr(>|t|)"])
  34. Realiza una regresión múltiple para ajustar una parábola a esos puntos. Consulta el capítulo 11 de Introducción a R para el modelo de la fórmula. 
    p <- lm (prod ~ densidad + I(densidad^2)) $ coef
  35. Representa la parábola (seq, lines) junto a la nube de puntos. 
    abscisas <- seq (min(densidad), max(densidad), length.out = 100)
ordenadas <- p[1] + p[2]*abscisas + p[3]*abscisas**2
plot (prod ~ densidad)
lines (abscisas, ordenadas, col=3)
  36. Representa gráficamente la relación entre alzada y peso para las vacas.csv. 
    v <- read.csv2 ("vacas.csv")
plot (v$peso ~ v$alzada)
v1 <- v [which (v$peso > 50 & v$alzada > 50 & v$alzada < 150), ] # "which" elimina los NA
plot (v1$peso ~ v1$alzada)
  37. Realiza una regresión lineal y otra parabólica de peso frente a alzada y representa las curvas ajustadas junto a la nube de puntos. 
    plot (v1$peso ~ v1$alzada)
l <- lm (peso ~ alzada, v1)
abline (l, col=2)
p <- lm (peso ~ alzada + I(alzada**2), v1) $ coef
x <- seq (min(v1$alzada), max(v1$alzada), length.out=100)
y <- p[1] + p[2]*x + p[3]*x^2
lines (x, y, col=3)
  38. Guarda la representación anterior como png y como pdf. 
    savePlot()       # guarda el gráfico actual como PNG
## empezando de cero:
png()
plot (v1$peso ~ v1$alzada)
abline (l, col=2)
lines (x, y, col=3)
dev.off()
pdf()
plot (v1$peso ~ v1$alzada)
abline (l, col=2)
lines (x, y, col=3)
dev.off()
  39. Crea una lista que contenga las siguientes componentes con estas etiquetas:
    R2
    coeficiente de determinación del ajuste cuadrático
    beta0
    ordenada en el origen
    sigbeta0
    p-valor asociado a \(\beta_0\)
    peso100
    predicción de peso cuando alzada vale 100 (ajuste parabólico) 
    s <- summary (lm (peso ~ alzada + I(alzada**2), v1))
lista <- list (R2       = s$r.squared,
beta0    = s$coef["(Intercept)","Estimate"],
sigbeta0 = s$coef["(Intercept)","Pr(>|t|)"],
peso100  = s$coef["(Intercept)","Estimate"] +
100 * s$coef["alzada","Estimate"] +
10000 * s$coef["alzada","Estimate"])
  40. Calcula la edad de cada animal usando substr. 
    edad <- function (nac, des) {
    dia.nac <- as.numeric (substr (nac, 1, 2))
    mes.nac <- as.numeric (substr (nac, 4, 5))
    ano.nac <- as.numeric (substr (nac, 7, 10))
    dia.des <- as.numeric (substr (nac, 1, 2))
    mes.des <- as.numeric (substr (nac, 4, 5))
    ano.des <- as.numeric (substr (nac, 7, 10))
    dia.des-dia.nac + 30*(mes.des-mes.nac) + 365*(ano.des-ano.nac)
}
edad1 <- edad (v1$f_nacim, v1$f_destete)
  41. Calcula la edad de cada animal usando as.Date. 
    edad2 <- as.Date (v1$f_destete, format="%d/%m/%Y") -
         as.Date (v1$f_nacim,   format="%d/%m/%Y")
  42. Definir una función asim que calcule el coeficiente de asimetría de Físher. 
    asim0 <- function (x)
    mean ((x-mean(x))**3) / sd(x)**3
asim00 <- function (x)
    mean ((x-mean(x,na.rm=TRUE))**3,na.rm=TRUE) / sd(x,na.rm=TRUE)**3
asim <- function (x, ...)
    mean ((x-mean(x,...))**3,...) / sd(x,...)**3
  43. Hallar una función en R que calcule el coeficiente de asimetría de Físher. 
    e1071::skewness
  44. Definir una función en R que
    • dibuje un histograma a partir de un vector numérico dado como primer argumento
    • admita un segundo argumento opcional, el color del histograma (por omisión, rojo)
    • admita un tercer argumento opcional, el número de barras (por omisión, diez); hay que forzar que el histograma tenga exactamente el número de barras pedido
    • cada barra debe mostrar por encima el porcentaje de observaciones que representa 
    Histograma <- function (v, color = 2, nbarras = 10) # "v" es el vector nume'rico
{
## calculamos los extremos horizontales de las barras
cortes <- seq (min(v), max(v), (max(v)-min(v))/nbarras)
## realizamos un histograma auxiliar para calcular "counts"
## no hacen falta los adornos (col, labels) pero si' fijar las barras
aux <- hist (v,
breaks = cortes) # fuerza el nu'mero de barras
## calculamos porcentajes redondeados
## suponemos que no hay datos ausentes; si los hubiera, en vez
## de «length(v)» puede usarse «sum(!is.na(v))» o «sum(aux$counts)»
porcentajes <- round (aux$counts / length(v) * 100)
etiquetas   <- paste (porcentajes, "%")
hist (v,
col    = color,
breaks = cortes,
labels = etiquetas) # tienen que ser cadenas; "paste" lo asegura
}
  45. Representa el peso en función de la edad con un gráfico. Añádele una recta de regresión. Guarda el gráfico como PNG y como PDF. 
    png()
plot (v1$peso ~ v1$edad)
abline (lm (peso~edad, v1), col=2)
dev.off()
pdf()
plot (v1$peso ~ v1$edad)
abline (lm (peso~edad, v1), col=2)
dev.off()
  46. Ídem, pero realiza un gráfico por cada año. 
    ano <- function (fecha) as.numeric (substr (fecha, 7, 10))
anos <- sort (unique (ano (v1$f_nacim)))
for (i in anos) {
    v2 <- v1 [ano(v1$f_nacim) == i, ]
    png()
    plot (v2$peso ~ v2$edad)
    abline (lm (peso~edad, v2), col=2)
    dev.off()
    pdf()
    plot (v2$peso ~ v2$edad)
    abline (lm (peso~edad, v2), col=2)
    dev.off()
}
  47. Fíjate en los datos mtcars. Convierte las variables mgp, disp, wt en unidades métricas.
  48. Haz un análisis descriptivo de cada variable.
  49. Haz un análisis del consumo en función de am.
  50. Haz un análisis del consumo en función de am y vs combinados.
  51. Haz un análisis del consumo en función de wt, disp y cyl.
  52. Por cada nivel de cyl, haz un gráfico del consumo respecto al peso, que incluya una recta y una parábola de regresión. 
    pdf()
for (i in unique(mtcars$cyl))
    {
        m <- mtcars[mtcars$cyl==i,]
        plot (consumo ~ peso, m, main=paste("num. cilindros =",i))
        abline (lm (consumo ~ peso, m))
        x <- 700:1500
        p <- coef (lm (consumo ~ peso + I(peso**2), m))
        y <- p[1] + p[2] * x + p[3] * x**2
        lines (x, y)
    }
dev.off()
  53. Crea un dataframe con las mismas variables que mtcars pero donde todas las variables estén tipificadas (con media 0 y desviación típica 1). 
    mtcars01 <- scale (mtcars)
  54. Ídem con mediana 0 y recorrido intercuartílico 1. 
    mtcars51 <- mtcars
for (j in 1:ncol(mtcars)) # for (j in names(mtcars))
    mtcars51[,j] <- (mtcars[,j] - median(mtcars[,j])) / IQR(mtcars[,j])
## otra forma:
mtcars51 <- apply (mtcars, 2, function (x) (x - median(x)) / IQR(x))
  55. Por cada país de PISA, crea un gráfico png y otro pdf que contenga:
    • un histograma de PERSEV
    • un histograma de OPENPS
    • la nube de puntos entre ambas variables

    Puedes generar los nombres de los ficheros usando paste o paste0. 

    download.file ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/pisa2012.csv.gz",
               "/tmp/pisa.csv.gz")
P <- read.csv (gzfile ("/tmp/pisa.csv.gz"), dec=",", na.strings="9999")
by (P, P$CNT,
    function (p)
    {
        png (paste0 (p$CNT[1], ".png"), width=500, height=1500)
        par (mfrow = c(3,1))
        hist(p$PERSEV); hist(p$OPENPS); plot(p$PERSEV,p$OPENPS)
        dev.off ()
        pdf (paste0 (p$CNT[1], ".pdf"))
        hist(p$PERSEV); hist(p$OPENPS); plot(p$PERSEV,p$OPENPS)
        dev.off ()
    })
  56. Por cada país de PISA, crea un dataframe con sus datos y guárdalo mediante save en un fichero .rda y mediante write.table en un fichero .csv.
  57. Por cada año y cada aeropuerto from de vuelos.csv, haz un gráfico de barras con los diez aeropuertos to más frecuentes.
  58. A partir de mtcars, obtén un dataframe que contenga los coches que
    • pesen más de 1.500 quilos y
    • consuman más de 15 litros por cada cien quilómetros

    o bien

    • su nombre empiece por T pero no contenga Corolla

    (Pistas: rownames, substr, grep.)

  59. Halla la desviación media de las variables de mtcars que sean cuantitativas continuas (digamos que tienen más de cuatro categorías).
  60. Halla la desviación media de las variables de vacas.csv que sean cuantitativas continuas (digamos que tienen más de cuatro categorías) condicionando a los datos de los últimos siete años. El resultado ha de ser en todo caso númerico (no NA).
  61. Carga los datos de paises.ods en R, habiéndolos exportado en LibreOffice a paises.csv dejando las opciones por omisión. 
    p <- read.csv ("paises.csv")
p$DENSIDAD <- as.numeric (gsub (",", ".", p$DENSIDAD))
  62. Basándote en mtcars, crea un dataframe llamado mtcoches donde
    1. El consumo venga en litros por 100 km.
    2. El peso venga en toneladas métricas.
    3. El cubicaje venga en centímetros cúbicos.
    4. Las variables vs y am sean factores con etiquetas adecuadas. 
    mtcoches          <- mtcars
mtcoches$consumo  <- 1 / (mtcoches$mpg * 1.6 / 3.9) * 100
mtcoches$peso     <- mtcoches$wt * 0.453
mtcoches$cubicaje <- mtcoches$disp * 2.54^3
## dos maneras de pasar a factor; una la probamos con VS y la otra con AM;
## la variable VS la trasformamos como carácter y al final la factorizamos:
mtcoches$vs [mtcoches$vs == 0] <- "V"
mtcoches$vs [mtcoches$vs == 1] <- "S"
mtcoches$vs <- factor (mtcoches$vs)
## la variable AS la factorizamos y luego cambiamos sus niveles:
mtcoches$am <- factor (mtcoches$am)
levels (mtcoches$am) <- c ("automático", "manual")
## para lo anterior, ¡ojo con mantener el orden de los niveles originales!
## eliminamos las variables originales:
mtcoches[c("mpg","wt","disp","vs","am")] <- NULL
  63. Guarda ese dataframe en un fichero en formato binario (.rda): save 
    save (mtcoches, file="mtcoches.rda")
  64. Sal de R sin guardar el entorno (por ejemplo, con q("no")), vuelve a entrar y recupera mtcoches a partir del fichero: load 
    load ("mtcoches.rda")
  65. Guarda ese dataframe en un fichero .csv: write.table, write.csv, write.csv2 
    write.csv2 (mtcoches, file="mtcoches.csv")
  66. Sal de R sin guardar el entorno (por ejemplo, con q("no")), revisa el contenido del fichero .csv, vuelve a entrar en R y recupera mtcoches a partir del fichero. 
    mtcoches <- read.csv2 ("mtcoches.csv")
  67. Crea un modelo lineal (llámalo l) que prediga el consumo a partir del peso y del cubicaje. 
    l <- lm (consumo ~ peso + cubicaje, mtcoches)
  68. Añade a mtcoches una columna con el consumo predicho por l. 
    mtcoches$consumo.predicho <- predict (l, mtcoches)
## otra forma:
mtcoches$consumo.predicho <- l$coef["(Intercept)"] +
                             l$coef["peso"] * mtcoches$peso +
                             l$coef["cubicaje"] * mtcoches$cubicaje
## otra forma:
mtcoches$consumo.predicho <- l$coef[1] +
                             l$coef[2] * mtcoches$peso +
                             l$coef[3] * mtcoches$cubicaje
  69. Calcula el error cuadrático medio de dicha predicción. 
    mean ((mtcoches$consumo - mtcoches$consumo.predicho) ^ 2)
  70. Para calcular el coeficiente de determinación en tanto por cien, he probado 
    100 * summary (l) ["r.squared"]

    pero me da error. ¿Por qué? 

    ## porque [] extrae una "sublista de tamaño 1" y me interesa "el elemento";
## puedo hacerlo así, con [[]]
100 * summary (l) [["r.squared"]] 
## o trasformar la sublista en vector, así:
100 * as.numeric (summary (l) ["r.squared"])
  71. Vamos a asignar a cada coche una supuesta fecha de fabricación, aleatoria. Añade a mtcoches una columna con fechas aleatorias desde el 1 de enero de 1970 al 31 de diciembre de 1979. 
    ## una posibilidad:
nr    <- nrow (mtcoches)
días  <- sample (1:28,      nr, replace=TRUE) # el menor mes tiene 28 días
meses <- sample (1:12,      nr, replace=TRUE)
años  <- sample (1970:1979, nr, replace=TRUE)
mtcoches$fecha <- as.Date (paste (años, meses, días, sep="-"))
## otra:
mtcoches$fecha <- sample (seq (as.Date ("1970-1-1"), 
                               as.Date ("1979-12-31"),
                               by="days"),
                          nrow (mtcoches),
                          replace = TRUE)
  72. Añade otra columna con el día de la semana correspondiente a cada fecha: weekdays 
    mtcoches$diasemana <- weekdays (mtcoches$fecha)
  73. Construye otra columna llamada hpu añadiendo a la variable potencia un ruido aleatorio con distribución uniforme entre \(-1\) y \(+1\). 
    mtcoches$hp.u <- mtcoches$hp + runif (nrow(mtcoches), -1, +1)
  74. Construye otra columna hpg añadiendo a la variable potencia un ruido aleatorio con distribución gausiana con media \(0\) y desviación típica \(5\). 
    mtcoches$hpg <- mtcoches$hp + rnorm (nrow(mtcoches), 0, 5)
  75. Genera una muestra de 10.000 valores gausianos con media 1'7 y desviación típica 0'05. Representa un histograma de dicha muestra y añádele una curva de densidad gausiana con los parámetros dichos. 
    muestra <- rnorm (10000, 1.7, 0.05)
h <- hist (muestra)
x <- seq (min(h$breaks), max(h$breaks), length.out=100)
lines (x, dnorm (x, 1.7, 0.05), col=2)
  76. El bucle 
    for (i in iris) print (summary (i))

    escribe los resúmenes descriptivos de cada variable, pero no sus nombres. Modifícalo para que ponga el nombre de cada variable antes del resumen. 

    for (i in 1:ncol(iris)) {cat ("-----------------\nVariable:", names(iris)[i], "\n"); print (summary (iris[[i]]))}
  77. Genera una matriz que contenga la tabla de multiplicar. 
    ## versión 1
tablamul <- matrix (NA, 10, 10)
for (i in 1:10) for (j in 1:10) tablamul[i,j] <- i * j
## versión 2
tablamul <- 1:10 %o% 1:10
## versión 3
tablamul <- outer (1:10, 1:10, "*")
  78. Echa un vistazo al fichero /var/log/dpkg.log. Cárgalo en R (read.csv, …). Al principio de cada renglón está la hora de cada registro. Crea un vector con las horas mediante strptime 1 y calcula su desviación típica y su recorrido. 
    a <- read.table ("/var/log/dpkg.log", sep="\t")  # no hay tabuladores: todo en una columna
horas <- strptime (a$V1, format = "%Y-%m-%d %H:%M:%S")
sd (horas)                                     # en segundos;
diff (range (as.numeric (horas)))              # recorrido en segundos;
diff (range (horas))                           # varía la unidad de medida;
attr (unclass (diff (range (horas))), "units") # para poder determinar la unidad;
  79. Obtén los minutos de cada registro del vector anterior: unclass 
    names (unclass (horas))
unclass (horas) $ min
  80. Definir una función que, dados los coeficientes «a,b,c» de una ecuación de segundo grado (\(a x^2 + b x + c = 0\)), devuelva sus raíces reales (sin parte imaginaria), en caso de haberlas. 
    resolver2grado <- function (a, b, c) {
  d <- b^2 - 4*a*c # discriminante
  if (d > 0) (-b + c(-1,1) * sqrt(d)) / (2 * a)
  else if (d == 0) -b / 2 / a
  else {warning ("sin solución real"); c()}
}
  81. Sea \(f\) la función de densidad de una distribución triangular simétrica tal que los extremos del soporte sean \(-a\) y \(a\).
    1. Definir en R la función f que calcule dicha densidad en un punto x dado como argumento. Por omisión, usar a=2.
    2. Comprobar si f puede aplicarse a vectores, por ejemplo, f(0:5). Si no, crear una versión fv vectorizada de f. 
    ### función de densidad Triangular (-a, a, 0)

## versión 1: con "if" y "else
f1 <- function (x, a=2)  # por omisión, a=2
{
if (x < -a)
0
else if (x <= 0)
(a+x) / a^2
else if (x <= a)
(a-x) / a^2
else 0
}
## versión 2: mediante un paquete
install.packages("EnvStats") # tras buscar en internet "r triangular distrib"
f2 <- function (x, a=2)
EnvStats::dtri(x, -a, a, 0)
## versión 3: multiplicando por booleanos (lógicos)
f3 <- function (x, a=2)
1/a^2 * ((a + x) * (x > -a) * (x <= 0) +
(a - x) * (x >  0) * (x <  a))

## vamos a vectorizar la definición de "f1"
fv1 <- Vectorize(f1)
## así ya se pueden hacer estas cosas:
xx <- seq (-3, 3, 0.1); yy <- fv1(xx); plot (xx, yy, type="l")
plot (fv1, -3, 3)
## otra opción sería usar "ifelse":
fv1 <- function (x, a=2)  
ifelse (x < -a,
0,
ifelse (x <= 0,
(a+x) / a^2,
ifelse (x <= a,
(a-x) / a^2,
0)))
  82. Ídem pero devolviendo siempre dos soluciones complejas. 
    resolver2grado <- function (a, b, c)
  (-b + c(-1,1) * sqrt (b^2 - 4 * a * c + 0i)) / (2 * a)
## o bien:
resolver2grado <- function (a, b, c)
  (-b + c(-1,1) * sqrt (as.complex (b^2 - 4 * a * c))) / (2 * a)
  83. Se define el fectorial del entero \(z\in\mathbf{Z}\) como \[ \mathcal F(z) = \prod_{i=z}^{-1}\,i %= i\cdot(i+1)\cdots(-2)\cdot(-1) \] o bien, recursivamente, como \[ \mathcal F(-1) = -1\qquad \mathcal F(z) = z\cdot\mathcal F(z') \] donde \(z'=z+1\) si \(z<-1\) y \(z'=z-1\) si \(z>-1\).

    Escribe una función en R que calcule el fectorial de su argumento. 

    ## versio'n iterativa
fectorial <- function (n) {
    producto <- 1
    for (i in n:-1) producto <- producto * i
    producto
}
## versio'n recursiva
fectorial <- function (n)
    if (n == -1) -1 else n * fectorial (n + ifelse(n>=0, -1, 1))
## versio'n vectorial
fectorial <- function (n) prod (n:-1)
  84. Dibujar la función de densidad de una distribución ji-cuadrado de 5 grados de libertad.

    Sombrear las dos colas del 5% de probabilidad cada una, es decir, del 0 al percentil 5 y del percentil 95 hasta «infinito». 

    minx <- 0; maxx <- qchisq(0.99,5)
plot (x <- seq(minx,maxx,0.01), dchisq (x, 5), type="l")
lines (x <- seq(minx,qchisq(0.05,5),0.001), dchisq (x, 5), type="h")
lines (x <- seq(qchisq(0.95,5),maxx,0.001), dchisq (x, 5), type="h")
  85. Hacer un diagrama de dispersión de la longitud de pétalo frente a la longitud de sépalo, con los datos del dataframe «iris» (si no estuviera, ejecutar «library(datasets);data(iris)»).

    Que la setosa aparezcan de color rojo; las versicolor, de color azul; y las virginica, de verde. 

    plot (Petal.Length ~ Sepal.Length, iris, col = rep (c("red","blue","green"), each=50))
  86. Definir una variable «m» que contenga en la celda «i;j» el valor de «i» \(\times\) «j». 
    m <- matrix (NA, 10, 10)
for (i in 1:10) for (j in 1:10) m[i,j] <- i * j
# otra opción
for (i in 1:10) m[i,] <- i * 1:10
  87. Definir una variable «m» que contenga por ejemplo en la celda «3;5» una cadena con el texto "3 × 5 = 15". 
    m <- matrix (NA, 10, 10)
for (i in 1:10) for (j in 1:10) m[i,j] <- paste (i, "×", j, "=", i * j)
  88. Definir una variable «d» que sea un dataframe de 1000 filas y que
    • la primera columna sean valores aleatorios de una uniforme U(0;1)
    • la segunda, aleatorios gausianos N(170;15)
    • la tercera, aleatorios exponenciales de esperanza 5
    • la cuarta, aleatorios puasones de esperanza 20
    • la quinta, letras aleatorias del alfabeto inglés (variable «letters») con la misma probabilidad
    • la sexta, vocales (a, e, i, o, u) aleatorias donde la «u» tenga el doble de probabilidad que cada una de las otras 
    unif <- runif (1000, 0, 1)
gaus <- rnorm (1000, 170, 15)
expo <- rexp  (1000, 1/5)
puas <- rpois (1000, 20)
letr <- sample (letters, 1000, replace=TRUE)
voca <- sample (c("a","e","i","o","u"), 1000, replace=TRUE, prob=c(1,1,1,1,2))
d <- data.frame (unif, gaus, expo, puas, letr, voca)
  89. Obtener sendos vectores de nombres de variables cuantitativas y cualitativas de «d». 
    cuantitativas <- cualitativas <- c()
for (i in 1:ncol(d))
  if (is.numeric (d[[i]]))
      cuantitativas <- c (cuantitativas, names(d)[i])
  else cualitativas <- c (cualitativas,  names(d)[i])
# otra forma
cuantitativas <- names (d) [sapply (d, is.factor)]
cualitativas  <- names (d) [sapply (d, is.numeric)]
  90. Definir una lista «l» de longitud 6 tal que su elemento «i»-ésimo sea el resultado de aplicar «summary» a la «i»-ésima variable de «d». 
    # una forma:
l <- list(); for (i in 1:6) l <- c (l, list (summary (d[,i])))
# otra forma:
l <- list(); for (i in 1:6) l[[i]] <- summary (d[,i])
# otra forma:
l <- lapply (d, summary)
  91. Añadir nombres a la lista del ejercicio anterior para identificar cada sumario con su variable. 
    names (l) <- names (d)
# otra forma (resultado distinto)
l <- list(); for (i in 1:6) l[[i]] <- list (names(d)[i], summary(d[,i]))
  92. Considerar los datos «mtcars». Obtener una lista que contenga por la media y desviación típica de «mpg» por cada valor de «cyl». 
    lista <- by (mtcars$mpg, mtcars$cyl, function (x) c (mean (x), sd (x)))
  93. Generar 10.000 valores de medias muestrales de muestras de tamaño 100 de una población gausiana N(170;15). Calcular la desviación típica de las medias muestrales. 
    medias <- apply (replicate (10000, rnorm (100, 170, 15)),
                 2,
                 mean)
sd (medias) # debería estar en torno a 1,5 = 15÷√100
  94. Definir una función que indique si su argumento es un número primo. Usarla para hallar los primos menores que 1000. 
    ## El 1 no es primo.
## El 2 es primo.
## Si n>2 entonces n es primo si es divisible sólo por 1 y por sí mismo.

## Para comprobar esto último, basta probar con divisores entre 2 y la raíz cuadrada de n.

## Un número d es divisor de n si   n %% d == 0   es decir, si el resto de la división es cero.

es.primo <- function (n)
{
    if (n == 1) FALSE
    else if (n == 2) TRUE
    else ! any (sapply (2:sqrt(n), function (d) n %% d == 0))
}

which (sapply (1:1000, es.primo))

## alternativa
install.packages ("primes")
which (sapply (1:1000, primes::is_prime))
  95. Cargar en R todos los datos disponibles en la carpeta datos, a saber:
    • mtcoches 
      load (url ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/mtcoches.rda"))
    • paises 
      paises <- read.csv ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/paises.csv", 
                    dec=",")
paises <- read.csv2 ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/paises.csv", 
                     sep=",") # equivale al anterior
    • pisa2012 
      system ("wget http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/pisa2012.csv.gz")
pisa2012 <- read.csv (gzfile ("pisa2012.csv.gz"), dec=",")
    • vacas 
      ### vacas.csv
vacas <- read.csv ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/vacas.csv", 
                   sep=";")
## sería equivalente porque no contiene números decimales:
vacas <- read.csv2 ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/vacas.csv")      

### vacas.txt
vacas <- read.table ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/vacas.txt", 
                     header = TRUE, 
                     fill = TRUE)
vacas <- read.csv2 ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/vacas.txt", 
                    sep = "") # equivalente al anterior

### vacas0.csv
vacas <- read.csv2 ("http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/vacas0.csv")
    • vuelos http://stat-computing.org/dataexpo/2009/the-data.html 
      cd ~/Descargas
wget http://bellman.ciencias.uniovi.es/~carleos/master/manadine/curso1/TIplE/ejercicios/datos/vuelos.csv.gz # descargar mediante HTTP
scp alumno@carleos2.epv.uniovi.es:/home/manadine/dat/vuelos.csv.gz .                                        # descargar mediante SSH
gunzip vuelos.csv.gz  # descomprimir
head vuelos.csv       # explorar estructura del fichero: es un CSV
wc -l vuelos.csv      # contar nu'mero de renglones = 1 + nu'mero de registros

      en nuestros ordenadores seri'a imposible cargar los datos de vuelos completos; lo haremos a trozos, procesa'ndolos segu'n un objetivo concreto; por ejemplo,

      1. mostrar distribución de vuelos por meses (Month)
      2. mostrar distribución de tiempos de retraso en llegadas (ArrDelay) 
      ## setwd ("datos") 
ntotal <- system ("wc -l vuelos.csv", intern = TRUE) # nu'mero de renglones
ntotal <- as.numeric (strsplit (ntotal, ' ') [[1]] [1])
ntrozo <- 1e6L # la L es necesaria para que aparezca 1000000 al usar en Bash:
system (paste ("split -d -a 3 -l", ntrozo, "vuelos.csv vuelos-trozo."))
vuelos <- read.csv ("vuelos-trozo.000")
campos <- names (vuelos)
meses  <- table (vuelos$Month)          # ejemplo de variable cualitativa
combinar <- function (vector1, vector2) # combinar dos vectores con nombres
    sapply (union (names (vector1), names (vector2)),
            function (i)
                sum (vector1[i], vector2[i], na.rm = TRUE))
## ejemplo de variable cuantitativa
v          <- na.omit (vuelos$ArrDelay)
ncuan      <- length (v)
minimo     <- min (v)
maximo     <- max (v)
suma1      <- sum (v)   # para media
suma2      <- sum (v^2) # para varianza
suma3      <- sum (v^3) # para asimetria
suma4      <- sum (v^4) # para curtosis
histograma <- hist (v, plot = FALSE,
                    breaks = c (-Inf, seq(minimo, maximo, length.out=101), Inf))
cortes     <- histograma$breaks
recuento   <- histograma$counts
## bucle
for (i in 1 : (ntotal %/% ntrozo))
    print (c (iter = i,
              system.time (
              {
                  ## no funciona leer a trozos desde R; cuando i = 9 se produce
                  ##    Process R buĉita at Sun Nov 17 21:47:40 2019
                  ## rm (vuelos); gc ()
                  ## vuelos <- read.csv ("vuelos.csv", nrows = ntrozo, header = FALSE,
                  ##                     skip = 1 + ntrozo * i)
                  vuelos         <- read.csv (paste0 ("vuelos-trozo.",
                                                      gsub(" ","0",
                                                           format (i, width=3))),
                                              header = FALSE)
                  names (vuelos) <- campos
                  meses          <- combinar (meses, table (vuelos$Month))
                  v              <- na.omit (vuelos$ArrDelay)
                  ncuan          <- ncuan + length (v)
                  histograma     <- hist (v, plot = FALSE, breaks = cortes)
                  recuento       <- recuento + histograma$counts
                  minimo     <- min (minimo, v)
                  maximo     <- max (maximo, v)
                  suma1      <- suma1 + sum (v)   # para media
                  suma2      <- suma2 + sum (v^2) # para varianza
                  suma3      <- suma3 + sum (v^3) # para asimetria
                  suma4      <- suma4 + sum (v^4) # para curtosis
              }) ["elapsed"]))
## descriptivos numéricos
media     <- suma1 / ncuan
media2    <- suma2 / ncuan
media3    <- suma3 / ncuan
media4    <- suma4 / ncuan
varianza  <- media2 - media^2
destipica <- sqrt (varianza)
## https://en.wikipedia.org/wiki/Skewness
asimetria <- (media3 - 3*media*varianza - media^3) / destipica^3
curtosis  <- (media4 - 4*media3*media + 6*media2*media^2 - 3*media^4) /
    destipica^4 - 3

### resultados

## distribución de la variable cualitativa (meses)

barplot (meses) # o bien, si no tenemos terminal gráfica:
sapply (round (1000 * prop.table (meses)), function (Ni) paste (rep("#", Ni), collapse=""))

####################################################################################### 10
    ################################################################################### 11
 ###################################################################################### 12
    ################################################################################### 1
           ############################################################################ 2
  ##################################################################################### 3
     ################################################################################## 4
   #################################################################################### 5
    ################################################################################### 6
 ###################################################################################### 7
 ###################################################################################### 8
      ################################################################################# 9

## distribución de la variable cuantitativa (retrasos)

mapply (function (Xi, Ni) paste (paste (rep("#", Ni), collapse=""), "≤", Xi), 
        cortes[-1], 
        round (100 * prop.table (recuento)))
## [55] " ≤ -41.0999999999999"                                                           
## [56] "##### ≤ -17.75"                                                                 
## [57] "############################################################ ≤ 5.60000000000014"
## [58] "####################### ≤ 28.95"                                                
## [59] "###### ≤ 52.3000000000002"                                                      
## [60] "## ≤ 75.6500000000001"                                                          
## [61] "# ≤ 99"                                                                         
## [62] "# ≤ 122.35"                                                                     
## [63] " ≤ 145.7"                                                                       

c (media, destipica, asimetria, curtosis) #  7.049963 30.750806  5.781361 92.913388
momento3 <- asimetria * destipica^3
momento4 <- curtosis * destipica^4 + 3
## hallemos distribución mejor ajustada a esos momentos:
## opción 1 - https://ca.wikipedia.org/wiki/Distribuci%C3%B3_SU_de_Johnson
library (SuppDists)
parámetros <- JohnsonFit (c (media, varianza, momento3, momento4),
                          "use") # da error de ajuste
plot (function (x) dJohnson (x, parámetros), -50, 150) 
## opción 2 - https://es.wikipedia.org/wiki/Expansiones_de_Edgeworth
library(PDQutils)
dapx_gca(0, c (media, varianza, momento3, momento4))

4 Día 23: SED, AWK, …

  1. Del fichero vacas.txt, modifica las fechas para que aparezcan en formato yanqui mes-día-año. 
    cat vacas.txt | sed -r 's#([0-9]{2})/([0-9]{2})/([0-9]{4})#\2-\1-\3#g'
  2. Del fichero vacas.txt, halla con AWK la desviación típica de las alzadas. 
    awk '($4>0&&$4<999){n++;suma+=$4;suma2+=$4^2}END{print sqrt(suma2/n-(suma/n)^2)}' vacas.txt
  3. Intenta cargar en R el fichero vacas.txt. Investiga qué ocurre mediante AWK e intenta corregirlo. 
    awk '{print NF}' vacas.txt | head
awk '(NF==11)' vacas.txt | head
awk '(NF!=11)' vacas.txt | head       # falta el genotipo
awk '(NF==11){print $0}(NF!=11){print $0,"GENOTIPODESCONOCIDO"}' vacas.txt > vacas-reparado.txt
  4. Crea un fichero listado.txt a partir de la salida del comando de Bash siguiente: 
    find /usr/share/doc -ls

    find /usr/share/doc -ls > listado.txt

    La columna segunda contiene la cantidad de bloques ocupados en KiB (quibioctetos).

    La columna sétima contiene el tamaño de los ficheros en octetos.

  5. Mediante cut, obtén un fichero octetos.txt que contenga sólo las dos columnas mencionadas arriba.

    Si no están bien alineadas las columnas en el fichero original, prueba a entubarlo a través del comando column -t 

    cat listado.txt | column -t | cut -c 10-15,45-55 > octetos.txt
  6. Haz lo mismo que en el anterior, pero mediante awk. 
    cat listado.txt | awk '{print $2,$7}' > octetos.txt
  7. Mediante awk, obtén la media y la desviación típica de cada columna. 
    awk '{n++;suma+=$1}END{print suma/n}' octetos.txt
awk '{n++;suma+=$1;suma2+=$1^2}END{print sqrt(suma2/n-(suma/n)^2)}' octetos.txt
awk '{n++;suma+=$2}END{print suma/n}' octetos.txt
awk '{n++;suma+=$2;suma2+=$2^2}END{print sqrt(suma2/n-(suma/n)^2)}' octetos.txt
  8. Obtén el máximo de cada columna. 
    sort -n -k1 octetos.txt | tail -n 1
sort -n -k2 octetos.txt | tail -n 1
  9. Mediante awk, obtén un fichero octetos+.txt que contenga, además, una columna con la diferencia entre octetos ocupados y octetos reales. 
    awk '{print $1,$2,1024*$1-$2}' octetos.txt > octetos+.txt
  10. Trasforma en Perl mediante a2p el último programa AWK: 
    echo '{print $1,$2,1024*$1-$2}' | a2p
  11. Compáralo con este programa en Python: 
    import fileinput,string
for renglon in fileinput.input():
  campos = renglon.split()
  print campos[0], campos[1], \
    string.atoi(campos[0])*1024 - string.atoi(campos[1])

    # - Perl requiere $ para nombrar variables.
# - Perl usa llaves para bloques; Python usa sangrías.
# - Perl lee ficheros y trata textos como números, de fábrica;
#   Python requiere cargar bibliotecas.
  12. Compáralo con este programa en Lisp 
    (defun dividir (cadena)
  (let ((s (make-string-input-stream cadena)))
    (loop for campo = (read s nil)
         while campo
         collect campo)))
(loop for renglon = (read-line *standard-input* nil)
   while renglon
   do (let ((campos (dividir renglon)))
        (format t "~@{~a ~}~%" 
                (car campos) 
                (cadr campos) 
                (- (* 1024 (car campos))
                   (cadr campos)))))

    # - Lisp usa paréntesis para la estructura del código fuente.
# - Lisp trae de fábrica lectura y escritura de ficheros.
# - La división de cadenas se puede hacer mediante una biblioteca,
#   pero es fácil definir una función para ello.

Nota al pie de página:

1

Quizá tengas que usar también Sys.setlocale("LC_TIME", "C") porque el nombre de los meses en syslog está en inglés y puede que la configuración de idioma (locale) de tu ordenador esté en español u otro idioma.

Autor: Carlos Enrique Carleos Artime

Created: 2021-01-19 mar 21:29

Emacs 25.2.2 (Org mode 8.2.10)

Validate