No sé si seguir leyendo libros. Sus autores los llenan de letras. Y es un lujo poder disponer del tiempo de leerlas todas.

Uno de esos libros llenos de letras es Linked, de Barabasi. Es un libro estupendo y recomendable. Pero podría ocupar 20 páginas si el autor fuese un poco más escueto y no se empeñase de llenarlo todo de anécdotas y colores.

Su primer capítulo trata sobre las redes sociales aleatorias, también conocidas como redes de Poisson o de Erdös-Rényi. Una de tales redes aleatorias es una colección de n nodos y enlaces entre ellos de manera que la probabilidad de que dos nodos x e y al azar estén unidos es p.

Dos propiedades matemáticas de ellas cita Barabasi. Hoy discutiré esas dos y una más.

La primera es que el grado de los nodos sigue una distribución de Poisson —y de ahí el nombre—. Efectivamente, la probabilidad de que un nodo tenga grado k tiene una distribución binomial B(n,p) y siempre que z = np sea pequeño, puede aproximarse por una distribución de Poisson de parámetro z.

El código

n.nodos <- 1000

build.network <- function( n.nodos, n.enlaces ){
  rs <- matrix( sample( n.nodos, 2 * n.enlaces, replace = T),
                nrow = n.enlaces, ncol = 2  )
  rs <- subset( rs,  ! rs[,1] == rs[,2] )

  if( nrow( rs ) < n.enlaces )
    return( rbind( rs , build.network( n.nodos, n.enlaces - nrow( rs ) ) ) )

  rs
}

link.x.node <- function( z, n.nodos = 1000 ){
  p <- z / n.nodos

  n.enlaces.potenciales <- n.nodos * ( n.nodos - 1 ) / 2

  # cuantos enlaces existen?
  n.enlaces <- rbinom( 1, n.enlaces.potenciales, p )
  kk <- build.network( n.nodos, n.enlaces )
  kk <- table( table( kk ) )
  cbind( z, as.numeric( names( kk )), kk )
}

node.dist <- function( z, n.iter ){
  res <- replicate( n.iter, link.x.node(z, n.nodos = n.nodos) )
  res <- mapply( cbind, res, 1:n.iter )
  res <- data.frame( do.call( rbind, res ) )
  colnames( res ) <- c( "z", "n.enlaces", "count", "iter" )

  boxplot( res$count ~ res$n.enlaces )
  tmp <- unique( res$n.enlaces )
  lines( tmp, n.nodos * dpois( tmp, lambda = z ), col = "red" )

  res
}

res <- node.dist( 5, 100 )

crea 100 redes aleatorias con n.nodos nodos y z = 5. Luego cuenta cuántos nodos tienen grado k y, finalmente, crea una gráfica similar a

que es un diagrama de cajas del grado sobre el que se superpone el número esperado de nodos de grado k si se acepta que su distribución es de Poisson. Como puede apreciarse, coinciden. Pero ya lo sabíamos porque conocíamos de viejo la relación entre las distribuciones de Poisson y la binomial.

La segunda propiedad de este tipo de redes es más enjundiosa. Tiene que ver con lo siguiente: en una red aleatoria, ¿existirán varias componentes inconexas? ¿o todos los nodos acabarán unidos a todos los nodos?

Supongamos que existe una componente gigante G. ¿Cuál será la probabilidad q de que alguno de los nodos x no pertenezca a ella? Si es un nodo con k vecinos , entonces

y dizque (yo no tengo aún claro el motivo: no acabo de ver por qué son sucesos independientes)

Si un nodo no pertenece a G, tampoco lo harán ninguno de sus vecinos. De ahí q sea igual a

y que, en general,

De ahí se obtiene la ecuación y si es la probabilidad de que un nodo pertenezca a G, se obtiene .

¿Tiene soluciones esta ecuación entre 0 y 1? ¿Para qué valores de z? Existe una raíz trivial s = 0: en tal caso el término de la izquierda y el de la derecha valen 0. En 1, el término de la izquierda vale 1 y el de la derecha, menos de uno.

La siguiente figura ilustra la situación:

En ella se representa el término de la izquierda,

curve( I, 0, 1, xlab = "", ylab = "" )

como una línea negra sólida y luego el término de la derecha con los valores z = 2

curve( 1 - exp( - 2 * x), add = T, lty = 2 )

y z = 0.5

curve( 1 - exp( - 0.5 * x), add = T, lty = 2 )

con líneas punteadas. Como podrá comprobar quién aún se derivan funciones, la derivada del término de la derecha en 0 es mayor o menor que uno dependiendo de si z es mayor o menor que 1 y eso obliga a que haya o no raíz de la ecuación en (0,1).

En resumen, cuando z < 1 no existe componente gigante y la red es una especie de sopa de letras. Pero si z > 1, aparece una supercomponente que contiene un porcentaje de los nodos dado por la solución de la ecuación anterior. Raíz que puede obtenerse con R así:

foo.optim <- function( x, z = 3 ) ( x - 1 + exp( -z * x) )**2
optimize( foo.optim, z = 2, interval = c(0,1))$minimum
optimize( foo.optim, z = 1.5, interval = c(0,1))$minimum

Para terminar esta larga entrada, vamos a ver si podemos fiarnos o no de Erdös. Porque una cosa es que este hombre nos diga que el tamaño de la componente grande es s y otra distinta, que lo sea.

Para eso vamos a crear muchas redes aleatorias con distintos valores de z y calcular el tamaño de la componente más grande para después compararlo con el calculado teóricamente. Y efectivamente, corriendo

library( igraph )
 
prop.giant.component <- function( z, n.nodos = 1000 ){
  p <- z / n.nodos
  n.enlaces.potenciales <- n.nodos * ( n.nodos - 1 ) / 2
  n.enlaces <- round( n.enlaces.potenciales * p )
  my.graph <- build.network( n.nodos, n.enlaces )
  my.graph <- graph.data.frame( my.graph )
 
  kk <- max( clusters( my.graph )$csize )
  kk / n.nodos
}
 
# prop.giant.component( 2 )
 
foo <- function( z, n.iter = 100 )
  replicate( n.iter, prop.giant.component( z ) )
 
zetas <- seq( 0.1, 5, by = 0.2)
 
res <- sapply( zetas, foo )
 
kk <- data.frame( res = as.numeric( t( res ) ), zetas = zetas )
 
plot( kk$zetas, kk$res, xlab = "", ylab = "" )
 
zetas <- unique( kk$zetas )
zetas <- zetas[ zetas > 1 ]
 
foo.optim <- function( x, z = 3 ) ( x - 1 + exp( -z * x) )**2
foo <- function( zeta ) optimize( foo.optim, z = zeta, interval = c(0,1))$minimum
theorical.sizes <- sapply( zetas, foo )
 
lines( zetas, theorical.sizes, col = "red")

se obtiene

con lo que nos damos por satisfechos y aguardamos a disponer de otro ratillo de asueto para sacarle la miga al capítulo II.

Incluso a los que conocemos el mercado desde dentro, la lectura de artículos como éste nos descubre un asombroso brave new world. Tanto los nuevos métodos con que dizque se afrontan los problemas más pedestres (como la detección de fraude, la retención de los mejores clientes, etc.) como la misma naturaleza de las áreas en las que se aplican (lucha antiterrorista, predicción de motines, elecciones sangrientas, actos de represión,… ¡e incluso el lanzamiento de cohetes por parte de Hizbolá!) parecen anunciar que ya tocamos lo que Asimov llamaba psicohistoria con la yema de los dedos.

No obstante, parece conveniente distinguir las promesas, tan pomposamente descritos en algunos medios, de las realidades, mucho menos noticiables y que apenas se dan a conocer tímidamente en letra chiquita. Por eso, recojo hoy dos noticias recientes con la que alimentar el sano escepticismo.

El primero es el del fin del sistema de códigos de alerta del departamento de interior de los EE.UU., que aspiraba a predecir la probabilidad de un ataque terrorista. La CNN, que hace lo que debe hacer un medio de comunicación, es decir, aportar contexto y análisis a lo escueto de la nota, afirma (mi traducción):

Los ataques fueron más frecuentes cuando el nivel era amarillo (riesgo significativo) que cuando era naranja (riesgo alto). La única vez en que fue rojo (riesgo severo), no pasó nada. Nunca ha sido azul o verde, los niveles que indican una menor probabilidad de riesgo.

La segunda noticia al respecto sobre los límites de este tipo de proyectos, en Wired.

Hace un tiempo comencé una serie de entradas, que serán finalmente tres, sobre los “mejores” paquetes de R. Esta va a ser la segunda entrega. Siento haber tardado tanto en realizarla: quienes me conocen saben que ocioso no he permanecido. De mis actividades de este periodo daré cumplida cuenta en entradas subsiguientes.

Tengo que añadir también como preámbulo que ha sido una conversación sobre análisis de redes sociales con un ex-compañero muy ducho en apropiarse de contraseñas ajenas la que me ha empujado finalmente a ahondar este estudio que tenía, junto a tantos, postergado en una esquina de mi disco duro.

En la primera entrada hice un ránking de los participantes más activos en la lista de correo r-help mediante técnicas análogas a las que usa Google para asignar pesos —o pageranks— a páginas. En ésta voy a analizar la evolución de estos pesos año a año.

En primer lugar, ahí va el número de personas que han respondido correos en la la lista en cada año (nótese que el 2010 está truncado en abril):

Luego, una vez calculados los pesos de cada una de ellas, muestro el podio en cada uno de los años:

Es constatable cómo Ripley, ”campeón absoluto” de manera global, parece haberse alejado de r-help (aunque sigue muy activo en la más elitista r-devel). También cómo Venables e Ihaka, pioneros en el mundo de R, dejaron pronto paso a otros entusiastas del lenguaje.

El único español que ha ocupado los lugares de privilegio es alobo, al que no conozco, pero que usaba una cuenta del CSIC a primeros de siglo.

Finalmente, presento una gráfica que muestra la evolución histórica del peso agregado de los respondedores por decil: se aprecia un claro proceso de concentración: una élite de expertos en R —creciente también en número, todo hay que decirlo— está concentrando a los líderes (en la terminología utilizada por los expertos en redes sociales) de la lista de correo de R.

(Los lectores más inquietos habrán advertido cómo ciertas series se cortan cuando no lo debieran en el gráfico anterior. Se trata de un artefacto causado por empates en pesos.)

Hace no mucho leí un articulillo de SAS sobre el impacto de ciertas marcas en determinadas redes sociales. Como este tema, así como sus posibles aplicaciones, siempre me ha intrigado, llevado de la curiosidad y del aburrimiento, decidí realizar un estudio análogo.

El artículo de SAS utiliza como materia prima resúmenes de publicaciones científicas que tratan de determinados medicamentos. A los autores les interesa conocer de qué marca de medicamentos escribe cada autor ponderando a éstos últimos en función de su impacto. El impacto lo miden a través de su peso en la red de colaboraciones científicas: tiene alto impacto un autor que ha escrito muchos artículos en colaboración con otros autores que también han escrito muchos artículos.

En defintiva, algo no muy distinto del famoso PageRank de Google: una página tiene un peso que se calcula ponderando el peso de las páginas que apuntan a ella mediante un algoritmo que, al menos a primera vista, parece recursivo.

En esta entrada voy a describir la primera fase de mi pequeño experimento. Éste consiste en asignar pesos a los distintos paquetes de R en función de su importancia en la lista de correo r-help. En la primera fase he asignado a los participantes en dicho foro un peso que mide su nivel de impacto. He aquí cómo:

Descarga del histórico de correos de r-help

Es sencillo: para descargar todos los correos por mes desde abril de 1997, basta con ejecutar (¿tengo que decir que uso Linux?)

wget -nd -r -l1 --accept gz https://stat.ethz.ch/pipermail/r-help/

Sacas la ropa de la lavadora, tiendes y, a la vuelta, ejecutas

zcat *.gz > all_mails

Preprocesamiento del fichero de correos

El fichero creado en el paso anterior ocupa 417MB (a fecha de hoy) y tiene formato mbox. Podría procesarse como un fichero de texto normal, pero Python, como en tantas ocasiones, acude en tu rescate. El módulo mailbox permite manipular este tipo de ficheros y en mi caso, con poco esfuerzo, crear dos tablas:

• Una contiene el ID de cada correo junto con el correo electrónico de su autor.
• Otra, la relación entre cada mensaje y el mensaje al que hace referencia, es decir, el mensaje del que es respuesta. En ella, obviamente, sólo se guardan los mensajes que son respuesta a otros mensajes anteriores.

La creación de estos ficheros, aunque conceptualmente simple, se complica por las excepciones, distintas configuraciones, etc. de los distintos sistemas de correo electrónico. Pero no son problemas que 20 líneas de código en Python no puedan solventar.

Asignación del “r-help-rank” a los participantes de la lista

Dada la naturaleza de la lista, decidí no asignar pesos sino por respuestas a mensajes, no por iniciar una conversación. Dos participantes están conectados si participan en una misma discusión. El peso, r-help-rank, de un participante lo construí como la suma de los pesos de los participantes que intervienen en una conversación común. Más concretamente:

1. En un paso inicial, asigno un peso idéntico a cada uno de los 7177 participantes (más bien, respondedores) en la lista.
2. Para cada participante, identifico las conversaciones en las que ha tomado parte y sumo los pesos de los restantes participantes en ellas.
3. Normalizo  los pesos (para que tengan suma 1) e itero el paso anterior hasta que el vector de pesos se estabiliza.

El algoritmo es similar al del cálculo por el método de las potencias del primer valor propio de la matriz cuya coordenada (i,j) cuenta el número de veces en que los miembros i y j de la lista participan en una misma conversación.

And the winner is… Ripley!

La tabla siguiente contiene los pesos de los 20 participantes con el r-help-rank más alto (sólo muestro parte del correo electrónico):

Que conste que “gb” en la posición vigésima no soy yo. De hecho, yo aparezco en la posición 452 de la lista. Si alguien tiene curiosidad por encontrarse en ella, que me deje un mensaje. No voy a publicar los resultados completos por no hacer públicos correos electrónicos de terceros sin permiso.

Caveats

Son muchas las objeciones que se le pueden realizar a este estudio. Soy consciente de muchas de ellas y es posible que, en función de la disponibilidad de tiempo, me disponga a resolver las que buenamente pueda. En tanto, estoy abierto (y de hecho, expectante) a atender comentarios y sugerencias de mis lectores.