Me quedé el otro día en el modelo probabilístico de los grafos aleatorios exponenciales. Quedaba una última parte y al ensayar su redacción me di cuenta de que me había metido en un huerto: la cosa es mucho más vasta de lo que a primera vista parecía.

Así que me limitaré a repasar lo más básico tratando de no meter demasiado la pata.

Tradicionalmente, se utilizaba para estimar los parámetros de un grafo la llamada técnica de la función de seudo-verosimilitud. Se ve que uno puede escribir

donde es una posible arista del grafo, es el grafo original sin la arista y es una función —la función difererencia— que depende del tipo de configuración. Un poco más en cristiano, que la razón de probabilidades para que exista un cierto vértice puede modelarse como una ecuación lineal en los coeficientes y eso permite estimarlos usando algo similar a regresiones logísticas.

Si vale el símil, sería equivalente a plantear un modelo logístico para predecir la existencia o no existencia de un determinado vértice.

Pero parece que esta técnica está mandada a recoger, i.e., desaconsejada, y suelen utilizarse en su lugar técnicas de Montecarlo. En particular, usando el paquete ergm de R. Como me superan los detalles de este asunto entero, me limitaré a recorrer con cierta minuciosidad el ejemplo que aparece en la función central del paquete, ergm.

Para ello, primero cargamos el paquete e importamos un conjunto de datos, flo.

library(ergm)
data(flo)

flo es una matriz de incidencia: contiene ceros y unos y estos últimos indican que una determinada pareja de familias florentinas mantuvieron en su día vínculos matrimoniales. A partir de dicha matriz se puede crear una red:

flomarriage <- network(flo, directed = FALSE)
flomarriage

Y también añadir atributos a sus nodos (en este caso, la riqueza relativa de las familias):

flomarriage %v% "wealth" <- c(10,36,27,146,55,44,20,8,42,103,48,49,10,48,32,3)
flomarriage

Obviamente, las redes pueden representarse gráficamente haciendo

plot( flomarriage )
plot(flomarriage, vertex.cex=flomarriage %v% "wealth" / 20)

para obtener (en la segunda sentencia) algo así como

La parte interesante llega ahora: plantear un modelo que, por ejemplo, indique si las familias tienen propensión a enlazarse cuando la diferencia de riqueza entre ellas es grande. Por supuesto, controlando por el número de enlaces totales que hay en el modelo, que es el significado del térmimo edges:

gest <- ergm(flomarriage ~ edges + absdiff("wealth"))
summary(gest)

La salida es

==========================
Summary of model fit
==========================

Formula: flomarriage ~ edges + absdiff("wealth")

Iterations: 20

Monte Carlo MLE Results:
Estimate Std. Error MCMC % p-value
edges -1.457666 0.354532 NA absdiff.wealth -0.004176 0.007387 NA 0.573
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Null Deviance: 166.355 on 120 degrees of freedom
Residual Deviance: 107.798 on 118 degrees of freedom
Deviance: 58.557 on 2 degrees of freedom

AIC: 111.8 BIC: 117.37

en la que diría yo que el coeficiente (negativo) de edges indica que la densidad del grafo es relativamente pequeña (que es un grafo con pocos vértices, vamos) y que descontado ese efecto, la diferencia de riqueza entre las familias no parece tener mayor efecto.

Los términos que aparecen a la derecha de la fórmula del modelo representan configuraciones y el paquete implementa muchas de ellas: kstar(n) para el número de estrellas de n vértices, transitive, degree, cycle,… y así hasta treinta o cuarenta de ellas.

Los lectores interesados en el asunto pueden encontrar gratificante el ejercicio consistente en relacionar el coeficiente obtenido para el parámetro edges al modelar redes creadas con sentencias del tipo

mi.red <- network(25, directed=FALSE, density=0.1)

haciendo variar el parámetro density.

Ayer dejamos abierto el problema de la inferencia en grafos. La idea fundamental es la de suponer que un grafo determinado no es tanto un grafo en sí como una realización de un proceso aleatorio de generación de aristas entre un determinado número de nodos.

El planteamiento es análogo al que se hace con las series temporales: no es tan importante la serie en sí como el hecho de que pueda probarse que obedece a un modelo autorregresivo, ARIMA, etc.

Una serie amplia de modelos de grafos aleatorios caen dentro de los llamados modelos exponenciales. Fijado el número de vértices, la probabilidad de observar el grafo (puede suponerse que) es

donde es una constante normalizadora, es un parámetro (similar al de las regresiones) y indica si la configuración A está o no presente en el grafo.

Una configuración es… un patrón. El concepto se comprende mejor a través de ejemplos. Y el más simple es el de una arista (que podría representar un vínculo de amistad entre Juan y Pedro). O un triángulo (Juan, Pedro y Tomás son mutuamente amigos). O cualquiera de los que aparecen en el siguiente gráfico:

El modelo probabilístico para un grafo en el que cada arista tiene una probabilidad dada de ocurrir sería, por lo tanto,

donde indica si existe o no la arista entre los nodos i y j. Si cada arista tiene la misma probabilidad, quedaría, todavía de manera más simple,

donde es el número de aristas. Si se quiere estudiar si el grafo tiene tendencia a organizarse creando triángulos (es decir, que el amigo de mi amigo tiene cierta propensión a ser también mi amigo), pueden plantearse modelos del tipo

donde es el número de triángulos que aparecen en el grafo.

Una vez planteado el marco probabilístico, queda plantearse la estimación de los parámetros y la inferencia, que permite encontrar resupuestas a preguntas del tipo:

• ¿muestra mi grafo una tendencia significativa a las relaciones recíprocas?
• ¿tiende mi grafo de una manera significativa a crear triángulos, i.e., relaciones transitivas?

Queda pendiente para una futura entrada.

Sea un colegio y sus alumnos. Sea el indicador de que el alumno i es amigo del alumno j. Con eso tenemos montado un grafo (o, si se prefiere, una red social).

Muchos análisis que se hacen sobre este tipo de redes son meramente descriptivos pero, ¿es posible la inferencia sobre este tipo de conjunto de datos?

Por ejemplo, en el grafo que describo más arriba, cabría preguntarse si hay reciprocidad, es decir, si es mucho mayor que . O dicho de otro modo, si el que Juan sea amigo de Pedro incrementa notablemente la probabilidad de que Pedro también se considere amigo de Juan.

¿De qué modo puede medirse este tipo de características en una red? ¿Qué tipo de parámetro puede estimarse (preferiblemente con sus intervalos de confianza) que muestre que en el grafo existe una predisposición a la reciprodidad?

Aparte de este tipo de relaciones, existen otras muchas que los expertos buscan en sus conjuntos de datos. Por ejemplo, algunas de las siguientes:

A este tipo de preguntas aspiran a dar respueta los llamados modelos exponenciales para grafos aleatorios, cuyos rudimentos expondremos mañana.

Una de las principales objeciones que se le pueden hacer a mi entrada de ayer es que puede estar confundiendo la causa con efecto: puede que parte de la radialidad de la red que obtuve tenga que ver con el tamaño desproporcionado de Madrid que, a su vez, podría haber sido causado por la radialidad de la red tradicional de las comunicaciones españolas.

Así que enviemos una partida de pescado en malas condiciones a Mercamadrid, convidemos a toda la provincia, veámosla fenecer víctima de contumaces diarreas y rehagamos la simulación suponiendo que

nodos.alt <- nodos
nodos.alt$pop[nodos.alt$prov == "Madrid"] <- 0

¿Qué forma tendría ahora la red? Ejecutando

res  <- do.call(rbind, apply(aristas, 1, function(x) peso.tramos(x[1], x[2], g2, nodos.alt)))
peso <- tapply(res$pop / (res$distancia)^(1), res$tramo, sum)
 
col <- peso
col[col < median(col)] <- 0
col <- rgb(0,0,0, 255 * col/max(col), maxColorValue=255)
plot.grafo(g2, nodos, col = col)
 
g3 <- delete.edges(g2,edges=E(g2)[peso < median(peso)])
col3 <- col[peso >= median(peso)]
plot.grafo(g3, nodos, col = col3)
 
E(g3)$weight <- peso[peso >= median(peso)]
 
centralidad <- data.frame(nodo = V(g3)$name, centralidad = alpha.centrality(g3) )
centralidad <- centralidad[order(-centralidad$centralidad),]
centralidad

se obtiene

Esta vez la península se parte en dos reeditando una suerte de Hispania Tarraconensis en la que las capitales más centrales son Tarragona, Valencia, Barcelona, Castellón, Jaén y Lérida.

Si en lugar de esta versión tan extrema su ponemos que Madrid tiene una poblacion promedio, es decir,

nodos.alt <- nodos
nodos.alt$pop[nodos.alt$prov == "Madrid"] <- median( nodos$pop )

se obtiene una configuración prácticamente similar:

Y, finalmente, si toda la población está distribuida uniformemente en las provincias, es decir,

nodos.alt <- nodos
nodos.alt$pop <- mean( nodos$pop )

las cosas cambian de manera bastante sorprendente:

La vieja Castilla serviría de nexo de comunicaciones, siendo las provincias con mayor centralidad Palencia, Burgos, Segovia, Guadalajara, Toledo, Soria y Madrid. Tal vez porque la densidad de capitales de provincia (al suponer la población igual en todas las provincias) favorece esa zona en que se agrupan de manera un poco más compacta.

¿Será que consigo que alguien se anime a introducir correcciones por superficie (suponiendo una distribución uniforme de la población por kilómetro cuadrado), PIB, acceso a puertos de mar o salidas a Portugal o Francia?

Me propuse hace un tiempo combinar lo que aprendí creando rutas callejeras por Zaragoza con una entrada que escribí sobre la estructura radial de las vías de transporte de España. El problema que me planteo es si tiene sentido que la red de carreteras Española tenga estructura radial habida cuenta de la geometría peninsular bajo ciertas hipótesis, siempre discutibles y mejorables, de partida.

Así que, en primer lugar, cargué los paquetes de R necesarios, un fichero que creé que contenía las capitales de provincia, su latitud, su longitud y la población de las respectivas provincias y fabriqué una red de carreteras muy ineficiente que unía todos los nodos entre sí:

library(maptools)
library(pxR)
library(igraph)
library( geosphere)
 
## datos: provincias y población
 
nodos <- read.table( "http://www.datanalytics.com/uploads/prov_pop_lat_lon.txt",
             sep = ",", dec = ",", header = T, encoding = "latin1")
 
## aristas
 
distancia <- function(p1, p2, nodos){
  alfa  <- c(nodos$lon[nodos$prov==p1], nodos$lat[nodos$prov==p1] )
  omega <- c(nodos$lon[nodos$prov==p2], nodos$lat[nodos$prov==p2] )
  distCosine( alfa, omega ) / 1000	# kms.
}
 
aristas <- expand.grid(nodos$prov, nodos$prov)
colnames(aristas) <- c("desde", "hasta")
aristas <- aristas[ as.character(aristas$desde) < as.character(aristas$hasta), ]
 
aristas$weight <- apply(aristas,1, function(x) distancia(x[1], x[2], nodos))
 
 
## grafo
 
grafo <- graph.data.frame(aristas, directed = F)
 
plot.grafo <- function(g, nodos, col = "black", text = F){
  tmp <- get.edges(g, V(g))
  vertices <- data.frame( (V(g)[get.edges(g, E(g))[,1]])$name,
                          (V(g)[get.edges(g, E(g))[,2]])$name, col = col )
 
  plot(nodos$lon, nodos$lat, xlab ="", ylab = "", xaxt = "n", yaxt="n")
  if( text )
    text(nodos$lon, nodos$lat,nodos$prov)
 
  apply(vertices, 1, function(x){
    x1 <- nodos$lon[nodos$prov == x[1]]
    y1 <- nodos$lat[nodos$prov == x[1]]
    x2 <- nodos$lon[nodos$prov == x[2]]
    y2 <- nodos$lat[nodos$prov == x[2]]
    lines( c(x1,x2), c(y1,y2), col = x[3])
  })
}
 
plot.grafo(grafo, nodos)	# pequeño caos

El resultado es este pequeño caos:

Por simplificar, eliminé todas las autovías que unían capitales de provincia cuando pudiera encontrar una ruta alternativa cuya longitud no excediese a la original por un factor de 1.2 haciendo

exceso.edge <- function(g, e){
  a <- shortest.paths(g)
  b <- shortest.paths(delete.edges(g, e))
  max( incr <- b / a, na.rm = T )
}
 
tmp <- sapply(E(grafo), function(e) exceso.edge(grafo,e))
g2  <- delete.edges(grafo, E(grafo)[tmp < 1.2])
 
plot.grafo(g2, nodos)

para obtener

Finalmente, simulé trayectos entre provincias con este criterio: una persona viaja de A a B con una probabilidad directamente proporcional al producto de las poblaciones de dichas provincias e inversamente proporcional a la distancia (en línea recta) entre ellas. La regla del producto de la población de las provincias es compatible con una muestra aleatoria de parejas de personas sobre la población total modificada en segunda instancia por la distancia entre ellas. Así que haciendo

peso.tramos <- function(a, b, g, nodos){
  data.frame(
    tramo = as.numeric(E(g2, path = get.shortest.paths(g2, from=a, to = b)[[1]])),
    pop = nodos$pop[nodos$prov == a] / 1e6 * nodos$pop[nodos$prov == b] / 1e6,
    distancia = distancia(a,b,nodos)
  )
}
 
res  <- do.call(rbind, apply(aristas, 1, function(x) peso.tramos(x[1], x[2], g2, nodos)))
peso <- tapply(res$pop / (res$distancia)^(1), res$tramo, sum)
 
## un primer gráfico
 
col <- peso
col[col < median(col)] <- 0
col <- rgb(0,0,0, 255 * col/max(col), maxColorValue=255)
plot.grafo(g2, nodos, col = col)

obtuve

En este mapa sólo se han representado la mitad de los tramos de mayor importancia (de acuerdo con el criterio arriba especificado) y en el resto se ha modulado la intensidad en función también de ese criterio.

¿Es una estructura radial? Podemos recurrir de nuevo a la teoría de grafos para medir la centralidad de los nodos:

g3 <- delete.edges(g2,edges=E(g2)[peso < median(peso)])
col3 <- col[peso >= median(peso)]
plot.grafo(g3, nodos, col = col3)
 
E(g3)$weight <- peso[peso >= median(peso)]
 
centralidad <- data.frame(nodo = V(g3)$name, centralidad = alpha.centrality(g3) )
centralidad <- centralidad[order(-centralidad$centralidad),]
centralidad

El resultado da preeminencia a Madrid y otras capitales de su entorno:

La cuestión es: ¿está Madrid en el centro a causa de su población? ¿Es esta población de Madrid grande entre otras cosas, gracias a la estructura radial de las comunicaciones? En una nueva entrega sobre este asunto volveré a analizar el problema con hipótesis de partida distintas.

Óscar Perpiñán me puso el otro día al tanto del paquete osmar de R, que

proporciona la infraestructura para acceder a datos de OpenStreetMap a través de diferentes fuentes, trabajar con ellos con R de una manera unificada y aprovechando la infraestructura que proporcionan otros paquetes como, por ejemplo, sp e igraph.

Hoy voy a ilustrar el uso de este paquete adaptando un ejemplo de sus autores para encontrar la ruta óptima entre dos puntos de Zaragoza, la mercería Bell y el colegio La Salle Montemolín, ambos lugares muy vinculados a mi prehistoria. Comenzaré cargando los paquetes necesarios y los datos de OpenStreetMap correspondientes a Zaragoza:

library( igraph )
library( osmar )
 
api <- osmsource_api(url = "http://api.openstreetmap.org/api/0.6/")
box    <- corner_bbox( -0.90, 41.60, -0.85, 41.69 )
zgz <- get_osm(box, source = api)

En segundo lugar, voy a usar las funciones auxiliares del paquete para extraer las calles del objeto zgz y representarlas geométricamente:

hways_zgz <- subset(zgz, way_ids = find(zgz, way(tags(k == "highway"))))
hways <- find(hways_zgz, way(tags(k == "name")))
hways <- find_down(zgz, way(hways))
hways_zgz <- subset(zgz, ids = hways)
 
plot_ways(hways_zgz, col = "gray" )
title("Calles de Zaragoza")

El resultado es

A continuación, selecciono los puntos incial y final de mi ruta:

id <- find(zgz, node(tags(v == "Moda infantil Bell")))[1]
hway_start_node <- find_nearest_node(zgz, id, way(tags(k == "highway")))
hway_start <- subset(zgz, node(hway_start_node))
 
id <- find(zgz, node(tags(v == "La Salle Montemolín")))[1]
hway_end_node <- find_nearest_node(zgz, id, way(tags(k == "highway")))
hway_end <- subset(zgz, node(hway_end_node))

Puedo representarlos mediante

plot_nodes(hway_start, add = TRUE, col = "red", pch = 19, cex = 1)
plot_nodes(hway_end, add = TRUE, col = "blue", pch = 19, cex = 1)

para obtener

Finalmente, utilizo la infraestructura proporcionada por el paquete igraph (para el análisis de redes sociales) para calcular la ruta más corta entre ambos puntos haciendo

gr_zgz <- as_igraph(hways_zgz)
route <- get.shortest.paths(gr_zgz,
                            from = as.character(hway_start_node),
                            to = as.character(hway_end_node))[[1]]
 
route_nodes <- as.numeric(V(gr_zgz)[route]$name)
route_ids <- find_up(hways_zgz, node(route_nodes))
route_zgz <- subset(hways_zgz, ids = route_ids)
route_zgz

y puedo finalmente representarla en el mapa mediante

plot_ways(route_zgz, add = TRUE, col = "black", lwd = 2)

para obtener el resultado final

Es incluso posible (véase esto) obtener una lista de las calles y nodos que hay que seguir para ir de uno de los puntos al otro.

Interesante, ¿eh?

Notas:

• No he conseguido calcular rutas en zonas más amplias: descargar los mapas de carreteras/calles de toda España, por ejemplo, no es tan sencillo como descargar los de zonas más pequeñas (una ciudad, por ejemplo). El API por defecto de OpenStreetMap, el que he usado, sólo permite descargar zonas que ocupen menos de la cuarta parte de un grado cuadrado. Existen procedimientos para extracciones más amplias, pero no las he probado.
• Para otras ubicaciones distintas de las que he probado (y que existen como tales en el objeto zgz como, por ejemplo, el campus de la Universidad de Zaragoza) el código que he mostrado falla de manera que no he podido todavía explicar.

He leído recientemente el artículo Social Ties and their Relevance to Churn in Mobile
Telecom Networks porque ilustra una técnica muy de moda: el análisis de redes sociales (SNA) en en ámbito de las telecomunicaciones y, en particular, la construcción de indicadores tempranos de baja (churn) de clientes de telefonía móvil. Más aún, permite rediseñar estrategias basadas en los resultados para retener clientes: al clasificarlos mejor usando técnicas de SNA, pueden diseñarse estrategias activas para aquellos que no sólo tienen una mayor predisposición a darse de baja sino, además, a arrastrar con ellos a parte de su entorno social.

El artículo, en resumen, introduce dos indicadores. El primero, p(k), es más ilustrativo que práctico: se trata de la probabilidad de que un cliente que tiene k conexiones —una conexión es alguien con quien el cliente ha hablado durante un determinado periodo— que se han dado de baja previamente se dé él mismo de baja. El gráfico siguiente muestra cómo p(k) es una función creciente de k. Sin embargo, el indicador puede no ser particularmente útil dado que, estoy seguro, el número de clientes para los que k > 1 es, casi seguro, muy pequeño.

En la segunda parte los autores construyen un modelo de propagación. Les interesa no sólo contar —y construir, de paso, probabilidades de corte frecuentista— sino explicar la dinámica y aprovecharla para construir modelos más útiles. La idea es la siguiente: un cliente que se da de baja transmite una señal a aquellos con los que se comunica. La señal puede ser del tipo esta compañía es malísima, me voy a ir a esta otra. No se sabe realmente cómo es la influencia, pero los autores la aproximan de la siguiente manera:

1. Asignan a cada cliente que se da de baja en un periodo determinado un cierto nivel de energía.
2. Un porcentaje de este nivel de energía se transmite de ellos a sus contactos en función de ciertos criterios (a mayor nivel de contacto, mayor flujo de energía). Este criterio preserva la energía: la energía total del sistema antes y después de la redistribución es la misma.
3. Los contactos que tienen un nivel de energía mayor que cero lo transmiten recursivamente a los suyos.
4. El proceso se itera hasta que alcanza un equilibrio razonable.

Al final, a muchos clientes (técnicamente, a los que pertenecen a la unión de las componentes conexas que contienen a las bajas) se les habrá asociado un nivel de energía. Y este nivel de energía es, según los autores, un indicador temprano de baja de alto valor predictivo.

¿Será?

Me he entretenido en los últimos tiempos tratando de responder una pregunta que, sin inquietarme, no dejaba de despertar mi curiosidad.

En la escuela nos enseñaron a definir palabras. Una de las primeras reglas de aquel juego era que el término definido no podía usarse en la definición: casa no se puede utilizar para definir casa. Los niños lo entendíamos. Sin embargo, los mayores hacían trampa: en el DRAE, construir se define en términos de edificar y edificar, en términos de construir.

Efectivamente, cójase el diccionario. El DRAE, por antonomasia. Búsquese una palabra. Cualquiera. En su definición aparecen otras. Búsquense estas a su vez. Y continúese recursivamente. Pueden pasar dos cosas:

• Volver a tropezar con la palabra original.
• No volver nunca a tropezar con ella.

Supongo que esas palabras que aparecen en los ciclos tienen una importancia léxica distinta de las del resto. Uno podría llamarlas palabras axiomáticas, palabras cuyo significado debería conocer el hablante antes de consultar la herramienta que define, es decir, el diccionario.

Quizás uno pueda contemplar la dicusión anterior de manera euclídea. Cada definición vendría a ser un teorema de la geometría euclidiana. La demostración de un teorema puede remitir a otros teoremas previos. Pero no indefinidamente: existen cinco postulados que se dan por buenos sin demostración, que se suponen ciertos de antemano.

Igualmente, en el diccionario, uno podría preguntarse cuáles son esos términos que se supone debieran darse por sabidos y que un diccionario euclídeo debiera abstenerse de definir. O, al menos, marcar explícitamente como tales.

Para conseguir mi objetivo he hecho lo siguiente:

1. Descargar la lista de palabras definidas en el DRAE, disponibles aquí y aquí.
2. Consultar (programáticamente, por supuesto) en el DRAE cada una de ellas.
3. Buscar la raíz de los términos que aparecen en la definición usando mi lematizador.
4. Crear una tabla con tres columnas:
   • lema
   • raíz del término que aparece en la definición
   • número de veces que aparece en la definición

Luego he analizado este conjunto de datos utilizando métodos de análisis de redes sociales. En efecto, considero que las palabras del diccionario, unas 88000, forman una red social en la que A es amiga de B si A aparece en la definición de B.

Y eso me permite responder una serie de preguntas. Por ejemplo, ¿cuántas palabras carecen de amigos? Es decir, ¿cuántas no aparecen en la definición de ninguna otra? Pues de las 87654 palabras del DRAE, son, exactamente, 51506, es decir, un 58.76 %. Incluyen desde a-, aarónico, aaronita, aba, ababa y ababillarse hasta zurumbo, zurupeto, zutanejo, zutuhil, zuzar y zuzo.

Las restantes 36148 palabras se usan en las definiciones de otras y van desde a, ababol, abacá, abacal, abacería, y ábaco hasta zurribanda, zurriburri, zurrón, zutano, zutujil y zuzón.

El siguiente paso del análisis consiste en eliminar aquellos términos que sólo entran en la definición de términos sin amigos o, recursivamente, en la definición de términos eliminados en el paso anterior. Por ejemplo, si A es un término sin amigos y B es un término que se usa en la definición de A y no en ningún otro, lo filtraría en este paso. Tras este filtrado, quedan 24683 términos, un 28.15 % de los términos originales. Seguro que Euclides pensaría que demasiados.

Los 11465 términos que se caen van desde ababol, abacal, abada, abajar, abakuá y aballar hasta zurriago, zurribanda, zurriburri, zutano, zutujil y zuzón. Encuentro en la lista términos como tridimensional, tropecientos, sobreexplotar, rutherfordio, presuntamente y perversión junto a otros como zangolotear, podrigorio, segueta, tolmera o tósigo, de cuya existencia a cabo de tener noticia.

De entre los restantes 24683 términos encontramos dos tipos. Por un lado, 197 familias aisladas de términos que son amigos entre sí, pero que no son amigos de otros términos. Por ejemplo, forman parte de estas familias parejas como violonchelista y violonchelo o triplas como tabulador, tabuladora y tabular.

Pero existe una familia de términos que comprende 24331 de ellos, un 27.75 % del total, que forma un clúster completo y que se extiende desde a, abacá, abacería, ábaco, abad y abadejo hasta zurdo, zuro, zurra, zurrador, zurrar y zurrón.

Quiero dejar constancia de que mis números son aproximados. Es posible que haya errores y, en efecto, he detectado algunos en el lematizador. Por ejemplo, este ha asignado a “nota” (musical) la raíz “notar” (verbo) en alguna ocasión. Etc.

No obstante, pienso que el número de términos en las definiciones es excesivo. ¿Debería la Academia esforzarse por reducir su número, por tratar de que la lista de palabras axiomáticas fuese más corta? Puede. Según algunos expertos, el número de palabras (en inglés) que utiliza activamente un hablante medio ronda las 20.000 y conoce (pasivamente) unas 40.000. Y estas cifras estarían dentro de los órdenes de magnitud que indico para el DRAE.

¿Cuál será, me pregunto, la opinión de mis lectores?

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.