¿Cómo instalar Ocaml para Windows con MSYS2?

En primer lugar instala MSYS2, para lo cual puedes seguir las instrucciones de instalación que hay en la página del programa http://www.msys2.org/, las intrucciones son muy claras.

Posteriormente confirmaremos si el paquete existe con el commando

$ pacman -Ss ocaml

Lo que imprimirá en pantalla:

Como puedes obervar, en el caso de OCaml, existen dos repositorios mingw32 y mingw64, los sufijos determinan la versión de 32 o 64bits respectivamente. De igual forma el prefijo mingw-w64-i686- es para la versión de 32-bits y; mingw-w64-x86_64- es para la versión de 64-bits.

Nota además que en este ejemplo, la versión 4.04.0-1 ya se encuentra instalada.

Una vez elegida la versión (32/64-bits), procedemos a ejecutar el siguiente comando:

pacman -S mingw-w64-i686_64-ocaml

para la versión de 64-bits o:

pacman -S mingw-w64-x86_64-ocaml

para la versión de 32-bits. La traza final de la instalación es la siguiente:

Y ahora verificamos que OCaml efectivamente esta instalado.

MyRandomPicker: Sortea ganadores de consursos

MyRandomPicker es una herramienta gratuita para efectuar ordenamientos aleatorios. De tal forma que ayuda en la elección fiable de los ganadores de algún concurso o rifa. Dado que esta desarrollado con base de código R, es necesario que tenga instalado el proyecto R (https://cran.r-project.org/).

Análisis de la Red Social Youtube con tuber

'tuber' es un paquete creado para enlazar con las API de YouTube desde R, es muy fácil de usar y en mi experiencia tiene un rendimiento eficiente. Para aprender más sobre las API de YouTube, dirígete a https://developers.google.com/youtube/v3/

Mapear archivo KML

Con datos de la Secretaría de Energía (SENER) se mapeará una capa de información en formato KML (del acrónimo en inglés Keyhole Markup Language). Los datos corresponden a las zonas en proceso de electrificación.

El código para extraer las coordenadas geográficas, la altitud así como el identificador (ID), es un aporte de Fabio Veronesi (https://www.r-bloggers.com/extract-coordinates-and-other-data-from-kml-in-r/), ligeras modificaciones fueron efectuadas. No obstante como veremos, los datos de la SENER sólo contienen latitud y longitud, la forma de hacerlo es la siguiente.

kml.text <- readLines('SENER.kml')
re <- "<coordinates> *([^<]+?) *<\\/coordinates>"  
coords <- grep(re,kml.text)  
   
re2 <- "src_id:"  
SCR.ID <- grep(re2,kml.text)  
   
re3 <- "<tr><td><b>Name:</b><td>"  
Name <- grep(re3,kml.text)  
   
kml.coordinates <- matrix(0,length(coords),4,dimnames=list(c(),c("ID","LON","LAT","ELEV")))  
kml.names <- matrix(0,length(coords),1)  
   
for(i in 1:length(coords)){  
    sub.coords <- coords[i]  
    temp1 <- gsub("<coordinates>"," ",kml.text[sub.coords])  
    temp2 <- gsub("</coordinates>"," ",temp1)  
    coordinates <- as.numeric(unlist(strsplit(temp2,",")))  
   
    sub.ID <- SCR.ID[i]  
    ID <- as.numeric(gsub("<tr><td><b>src_id:</b><td>"," ",kml.text[sub.ID]))  
   
    sub.Name <- Name[i]  
    NAME <- gsub(paste("<tr><td><b>Name:</b><td>"),"",kml.text[sub.Name])  
   
    kml.coordinates[i,] <- matrix(c(ID,coordinates),ncol=4)  
    kml.names[i,] <- matrix(c(NAME),ncol=1)  
}

Con esta parte del código extraemos la información del archivo KML.

> head(kml.coordinates)
     ID       LON      LAT ELEV
[1,] NA -92.13972 17.18167    0
[2,] NA -91.99667 17.18222    0
[3,] NA -91.20667 16.37694    0
[4,] NA -91.74639 17.03611    0
[5,] NA -92.11972 17.09611    0
[6,] NA -92.80972 17.40472    0

Como vemos, la información proporcionada por la dependencia se limita a las coordenadas geográficas. Por lo que, en este ejemplo, dejaremos la matriz sólo con las columnas que contienen información.

kml.coordinates <- kml.coordinates[,2:3]

El siguiente paso será plotear estas coordenadas, el código empleado es el siguiente:

library(RgoogleMaps)
z0 <- 5                     # zoom
y0 <- 21.525818             # Latitud zero
x0 <- -100.522399           # longitud zero


N <- length(kml.coordinates)
color <- c()                # Inicializo vector
color[1:N] <- 'gray20'

PuebMap  <- GetMap(center= c(y0,x0), zoom=z0, 
             GRAYSCALE = F, size=c(640,640),
             destfile = file.path(tempdir(),"sener.png"),
             maptype = 'terrain', SCALE = 2,)

PlotOnStaticMap(PuebMap)

{PlotOnStaticMap(PuebMap, lat = kml.coordinates[,2], 
                lon = kml.coordinates[,1], cex=1, pch=19,
                destfile = "sener.png", 
                col= color ,
                add=TRUE)};

Esto produce el siguiente mapa:

Esta representación nos permite observar que los procesos de electrificación se están llevando a cabo principalmente en las áreas montañosas.

Análisis Pokemon

La base de datos fue reunida de diversas fuentes por Alberto Barradas (https://www.kaggle.com/abcsds), la cual incluye 721 pokemon con algunas de sus estadísticas básicas:

> pokemon <- as.data.frame(read.csv('Pokemon.csv', header=T))
> str(pokemon)
'data.frame':   800 obs. of  13 variables:
 $ X.        : int  1 2 3 3 4 5 6 6 6 7 ...
 $ Name      : Factor w/ 800 levels "Abomasnow","AbomasnowMega Abomasnow",..: 81 330 746 747 103 104 100 101 102 666 ...
 $ Type.1    : Factor w/ 18 levels "Bug","Dark","Dragon",..: 10 10 10 10 7 7 7 7 7 18 ...
 $ Type.2    : Factor w/ 19 levels "","Bug","Dark",..: 15 15 15 15 1 1 9 4 9 1 ...
 $ Total     : int  318 405 525 625 309 405 534 634 634 314 ...
 $ HP        : int  45 60 80 80 39 58 78 78 78 44 ...
 $ Attack    : int  49 62 82 100 52 64 84 130 104 48 ...
 $ Defense   : int  49 63 83 123 43 58 78 111 78 65 ...
 $ Sp..Atk   : int  65 80 100 122 60 80 109 130 159 50 ...
 $ Sp..Def   : int  65 80 100 120 50 65 85 85 115 64 ...
 $ Speed     : int  45 60 80 80 65 80 100 100 100 43 ...
 $ Generation: Ord.factor w/ 6 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Legendary : Factor w/ 2 levels "False","True": 1 1 1 1 1 1 1 1 1 1 ...

• X.: ID para cada pokemon
• Name: Nombre de cada pokemon
• Type.1: Cada pokemon tiene un tipo, el cual determina vulnerabilidades o fortalezas a los ataques.
• Type.2; Algunos pokemones son tipo dual y tienen un segundo "tipo".
• Total: Suma de las estadísticas que se listan a continuación, da una idea general sobre que tan fuerte es el pokemon.
• HP: Hit points, o salud, define que tanto daño puede soportar un pokemon después de pelear
• Attack: Modificador de base para ataques normales.
• Defense: Resistencia al daño base contra ataques normales.
• SP Atk: Ataques especiales, el modificador de base para ataques especiales como rafaga de fuego, etc.
• SP Def: Resistencia al daño base contra ataques especiales.
• Speed: Determina que pokemon ataca primero cada ronda.
• Generation: La generación a la que pertenece el pokemon.
• Legendary: Variable booleana que indica TRUE cuando el pokemon es catalogado como legendario.

Promedio "Total" por tipo (1) de pokemon y generación

pkmn <- melt(with(pokemon, tapply(Total, list(Generation,Type.1), mean)))

pkmn <- pkmn[,c(2,1,3)]

ina  <- which(is.na(pkmn$value))

pkmn <- pkmn[-ina,]

colnames(pkmn) <- c('Type','Generation','Total')

freq <- count(pokemon, c('Type.1','Generation'))

pkmn$Freq <- freq[,3]

rownames(pkmn) <- seq(1,nrow(pkmn))

pkmn

p <- ggplot(pkmn) +

    geom_point(aes(x      = as.character(Type),

                   y      = as.character(Generation),

                   size   = Freq, 

                   color  = Total)) +

    scale_x_discrete(name = 'Type') +

    scale_y_discrete(name = 'Generation') +

    scale_colour_gradient(name='Avg Total', low='yellow', high='red3') +

    theme_bw() + 

    theme(axis.text.x=element_text(angle = 90, hjust = 0, vjust = 0.5)) + 

    ggtitle('Pokemons by Generation, type \n and Total Score')

p

Los pokemon tipo "Dragón" tienen en promedio los mejores puntajes, por otro lado los pokemones tipo "Agua" son los más abundantes.

Tipos de pokemon y puntaje "Total"

pkmn.1 <- melt(with(pokemon, tapply(Total, Type.1, mean )))

colnames(pkmn.1) <- c('Type.1','Total')

pkmn.1$freq <- melt(table(pokemon$Type.1))[,2]

mypalette <- rev(rainbow(20))

treemap(pkmn.1, index=c('Type.1'), vSize='freq', vColor='Total', 

        palette= mypalette,type='value', title.legend='Avg Total', 

        title='Pokemons by type and Total')

Clasificando pokemones legendarios

Convertimos la variable "Generation" a escala ordinal

pokemon$Generation <- as.ordered(pokemon$Generation)

Dividimos la base de datos original, en una muestra de entrenamiento con el 60% de los datos y el resto será la muestra de evaluación del modelo:

set.seed(1111)

id.train <- sample(seq(1,nrow(pokemon)), round(nrow(pokemon)*0.6))

train <- pokemon[id.train,]        # Muestra entrenamiento

test  <- pokemon[-id.train,]       # Muestra evaluación

El modelo sólo considerará las variables numéricas y la variable ordinal que definimos previamente.

> form <- as.formula(paste('Legendary', paste(names(train[c(5:12)]), 

+                    collapse=' + '), sep=' ~ '))

> form

Legendary ~ Total + HP + Attack + Defense + Sp..Atk + Sp..Def + 

    Speed + Generation

Regresión logística

> log <- glm(form, data=train, family = binomial())

> summary(log)

Call:

glm(formula = form, family = binomial(), data = train)

Deviance Residuals: 

     Min        1Q    Median        3Q       Max  

-2.41554  -0.14484  -0.03142  -0.00470   2.42037  

Coefficients: (1 not defined because of singularities)

               Estimate Std. Error z value Pr(>|z|)    

(Intercept)  -22.619056   3.338589  -6.775 1.24e-11 ***

Total          0.047109   0.011872   3.968 7.25e-05 ***

HP            -0.021168   0.015201  -1.393   0.1638    

Attack        -0.026050   0.015030  -1.733   0.0831 .  

Defense       -0.008286   0.011818  -0.701   0.4832    

Sp..Atk       -0.005889   0.014302  -0.412   0.6805    

Sp..Def       -0.002507   0.013974  -0.179   0.8576    

Speed                NA         NA      NA       NA    

Generation.L   1.660818   0.768322   2.162   0.0306 *  

Generation.Q  -0.438145   0.718330  -0.610   0.5419    

Generation.C   0.113477   0.711620   0.159   0.8733    

Generation^4  -0.715514   0.640430  -1.117   0.2639    

Generation^5  -0.565000   0.570013  -0.991   0.3216    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 284.85  on 479  degrees of freedom

Residual deviance: 112.57  on 468  degrees of freedom

AIC: 136.57

Number of Fisher Scoring iterations: 8

La fortaleza general del pokemon es el mejor predictor de los pokemones legendarios, seguido del ataque, y en menor medida la resistencia al daño. De igual forma, la primera generación de pokemones, se distingue significativamente de la 6a generación en cuanto a la cantidad de legendarios.

> cm.l <- table(train$Legendary, ifelse(fitted(log)>0.5,1,0));cm.l

       

          0   1

  False 430   8

  True   17  25

> cr.l <- sum(diag(cm.l)/sum(cm.l));cr.l             # Razón de exactitud

[1] 0.9479167

> rsq<- (1 - (log$deviance/log$null.deviance));rsq   #Pseudo R cuadrada

[1] 0.6047947

La matriz de confusión sobre la muestra de entrenamiento, y la razón de "exactitud" indican que es un buen modelo, no así el valor de la pseudo-R cuadrada. En cuanto al pronóstico sobre la muestra de evaluación:

> cm.l <- table(test$Legendary, ifelse(test$score.l>0.5,1,0));cm.l

       

          0   1

  False 292   5

  True   10  13

> cr.l <- sum(diag(cm.l)/sum(cm.l));cr.l

[1] 0.953125

>    pred.l      <- prediction(test$score.l,test$Legendary)

>    perf.l      <- performance(pred.l,"tpr","fpr")

>    max(attr(perf.l,'y.values')[[1]]-attr(perf.l,'x.values')[[1]])

[1] 0.9360269

> auc.l <- paste('Logistic=',round(performance(pred.l,'auc')@y.values[[1]],3))

> auc.l

[1] "Logistic= 0.984"

De acuerdo con la matriz de confusión, la razón de "exactitud", el estadístico de Kolmogórov-Smirnov y el Área Bajo la Curva, se puede apreciar que es un buen modelo.

A pesar de lo bueno que resulta ser el modelo de clasificación de máxima entropia, o de regresión logística, en el siguiente gráfico se muestran los resultados del desempeño predictivo de otros dos modelos de clasificación: Splines de Regresión Adaptativa Multivariante (MARS por sus siglas en inglés) y AdaBoost con el algoritmo SAMME.

clr <- palette(rev(rainbow(3)));clr <- palette(rev(rainbow(3)))

plot(perf.l, col=clr[1], lty=6, lwd=2,        # LOG

     xlab='Tasa Falso Positivo',ylab='Tasa Verdadero Positivo')

plot(perf.m, add=T, col=clr[2],lty=5, lwd=2)  # MARS

plot(perf.a, add=T, col=clr[3],lty=4, lwd=2)  # ADABOOST

abline(a=0, b= 1)

legend(0.20,0.80,                                 

      c(auc.l,auc.m,auc.a),                       

      col=clr, lwd=3, lty= c(6,5,4))              

grid(col='gray85')

GoogleMaps + Bubble plot en R

En esta entrada muestro cómo graficar referencias geográficas empleando GoogleMaps y R. Para empezar necesitamos referencias geográficas que deseemos plotear, en mi caso graficaré unidades económicas de la Comisión Federal de Electricidad de los estados de Tlaxcala y Puebla, los datos los obtuve del Directorio Estadístico Nacional de Unidades Económicas (DENUE) del INEGI.

Instalamos y cargamos la librería

install.packages('RgoogleMaps')
library(RgoogleMaps)

Posterior a ello, definimos algunos parámetros como el zoom de nuestro mapa, la latitud y la longitud centro, y adicional a ello una distancia (Euclidiana) que será el radio de círculos concéntricos a partir de nuestro centro.

z0 <- 8                     # zoom
y0 <- 19.040442             # Latitud zero
x0 <- -98.250160            # longitud zero
k  <- 0.4                   # Define el salto del circulo
                            # Saltos de 11 Km por cada décima

Colocaremos en el centro del mapa (o cualquier coordenada), un marcador que nos indique alguna referencia adicional.

markers = paste0("&markers=color:blue|label:I|19.040442,-98.250160")

Acto seguido, obtener nuestro mapa base:

PuebMap  <- GetMap(center= c(y0,x0), zoom=z0, GRAYSCALE = TRUE,
             size=c(640,640),markers=markers,
             destfile = file.path(tempdir(),"interm.png"),
             maptype = 'mapmaker-hybrid')

Enseguida, importamos nuestra base de datos de las unidades económicas, la cual deberá tener una estructura parecida a esto:

> head(DF)
       id      lat       lon        x1 x2
1 3535207 18.87717 -97.72655   87.5724  3
2 3398690 18.49201 -97.40435 1167.6320 40
3 3263793 19.04774 -98.04557   87.5724  3
4 3467483 19.12102 -98.24399 1167.6320 40
5 3316191 18.95277 -97.91986   87.5724  3
6 3457112 19.06409 -98.23185   87.5724  3

x1 y x2 representan variables cuantitativas a representar mediante un gráfico de burbujas. Lo siguiente es elaborar el gráfico sobre nuestro mapa base.

tpm <- bubbleMap(DF,coords=c("lon","lat"),

                 crs=sp::CRS("+proj=longlat +datum=WGS84"),

                 max.radius=500,legendLoc="bottomright",

                 filename = file.path(tempdir(),"interm.png"),

                 map=PuebMap,zcol='zinc',

                 key.entries = c(100,250,500,750,2500,6000));

Lucirá de la siguiente manera:

Lo siguiente que haremos será graficar n círculos concentricos de radio n*k en nuestro mapa, con las siguientes líneas de comando:

r1 <- k*1

r2 <- k*2

r3 <- k*3

n  <- seq(0,359,0.05)

a  <- seq(1,360,0.05)

theta <- (360*n)/a

     x1 <- x0 + r1*cos(theta)

     y1 <- y0 + r1*sin(theta)

     x2 <- x0 + r2*cos(theta)

     y2 <- y0 + r2*sin(theta)

     x3 <- x0 + r3*cos(theta)

     y3 <- y0 + r3*sin(theta)

circle <- data.frame(c(x1,x2,x3),c(y1,y2,y3))

N <-length(x1)                        # Tamaño del vector

color <- c()                          # Inicializo vector

color[1:N]           <- 'gold'

color[(N+1):(2*N)]   <- 'darkorange'

color[(2*N+1):(3*N)] <- 'red'

{PlotOnStaticMap(PuebMap, lat = circle[,2], 

                lon = circle[,1], cex=0.5, pch=19,

                destfile = "interm.png", 

                col= color ,add=TRUE)};

Quedará algo como esto:

Podemos definir un nuevo nivel de zoom:

Por último y con fines de análisis geoespacial, para definir en que radio se encuentra cada una de las Unidades Económicas (UE), empleamos las siguientes líneas:

D <- dim(DF)[1]; D
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 1 o no
DF$PotR1 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r1^2){
      DF[j,'PotR1'] <- 'dentro'
      } else DF[j,'PotR1'] <- 'fuera'
}
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 2 o no
DF$PotR2 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r2^2){
      DF[j,'PotR2'] <- 'dentro'
      } else DF[j,'PotR2'] <- 'fuera'
}
# Este vector me indica si la Unidad Económca se encuentra en el área del radio 3 o no
DF$PotR3 <- NULL
for(j in 1:D){
    if (((DF[j,2] - y0)^2 + (DF[j,3] - x0)^2) <= r3^2){
      DF[j,'PotR3'] <- 'dentro'
      } else DF[j,'PotR3'] <- 'fuera'
}

Finalmente, en este ejemplo, se analiza a las variables cuantitativas en función de su pertenencia o no al área definida por cada radio, tal como el número de posibles empleados en cada radio:

sumER_1 <- 0
for (j in 1:D){
if (DF[j,'PotR1']=='dentro'){
    sumER_1 <- DF[j,'pers']+ sumER_1
    }
}
sumER_1                       # Potencial Radio 1

sumER_2 <- 0
for (j in 1:D){
if (DF[j,'PotR2']=='dentro'){
    sumER_2 <- DF[j,'pers']+ sumER_2
    }
}
sumER_2 <- sumER_2 - sumER_1             
sumER_2                       # Potencial Radio 2

sumER_3 <- 0
for (j in 1:D){
if (DF[j,'PotR3']=='dentro'){
    sumER_3 <- DF[j,'pers']+ sumER_3
    }
}
sumER_3 <- sumER_3 - sumER_2 - sumER_1   
sumER_3                       # Potencial Radio 3


sumER_x <- sum(DF['pers']) - sumER_1 - sumER_2 - sumER_3
sumER_x                       # Potencial fuera de cualquier radio

Información que finalmente representamos en una tabla:

> potencial <- data.frame(rbind(c(sumER_1, sumER_2, sumER_3, sumER_x)

+                               )

+                         )

> colnames(potencial) <- c(' Radio 1', ' Radio 2', ' Radio 3', ' Radio 4')

> rownames(potencial) <- 'Empleados'

> round(potencial)

           Radio 1  Radio 2  Radio 3  Radio 4

Empleados     1539      372      398      128

Lo que nos indica que existen aproximadamente 1539 empleados en el radio 1, que enfatizemos es de 44 Km a partir de la coordenada central que hemos fijado.

Visualización de interacciones de Twitter con R (1 de 2)

En esta entrada muestro la visualización de interacciones de Twitter, con ayuda de R y tomando sólo como ejemplo el hashtag "#PRI", una vez que han pasado las elecciones.

Los paquetes necesarios son "graphTweets", "networkD3" y "igraph". Se emplea la información de 1,000 tweets del 5 de junio de 2016. Nuestra red lucirá así:

 Y si analizamos la formación de comunidades, nuestra red se verá de la siguiente forma:

Existen otros paquetes de R con los cuales también podemos graficar las interacciones como "rgl" o "tcltk".