Nuestros vecinos portugueses acaban de abrir un foro para sus usuarios de R un poco al estilo de nuestro r-help-es (¿todavía no estás dado de alta en él?).

Espero que les sirva de base para organizar una comunidad vibrante de usuarios. Y que algún día podamos organizar unas Jornadas Ibéricas de Usuarios de R.

Cuando tenía 18 años, pensaba, llegué a aprender todo lo que había que saber sobre factorización de matrices. Incluida la inutilidad de Jordan. El otro día, con un ciento y pico por ciento más de años, he descubierto una clase entera de factorizaciones que aquellos planes de estudios viejunos no contemplaban y que, ¡carajo!, aparte de útiles engarzan con otras ideas la mar de interesantes.

Se trata de factorizaciones positivas de matrices igualmente positivas.

Dentro de unos días os copiaré aquí unas líneas de un artículo del 83 que bien pudiera haber sido escrito el mes pasado. Pero hoy no voy a ir tan lejos. Me quedo con uno del 89 que recomiendo que hojeéis: Clinical vs Actuarial Judgement.

No, no vais a aprender en él nada que no sepáis. Os podrá parecer viejuno el uso de clínico o actuarial para denotar conceptos que ahora conocemos por otros nombres. Pero el tema, desgraciadamente, sigue siendo reciente. Todavía hay clínicos —más bien, todavía siguen siendo mayoría— que desconfían de los cálculos probabilísticos proporcionados por las máquinas (actuariales en la terminología del artículo) más que de su propia intuición.

Anuncio algo que no he conseguido hacer: agrupar grafos por topología. Pero no me he quedado lejos. Y espero que si alguien tiene alguna idea al respecto, nos lo haga saber al resto en la coda.

Contexto (disfrazado). Hay usuarios que tienen correos electrónicos. La relación esperada es de uno a uno. Pero la realidad es, como siempre, mucho más compleja: hay usuarios que tienen varios correos y correos compartidos por varios usuarios.

Hoy he escrito

last.date <- max(Filter(function(x) format(x, "%m") == "03",
	all.filled.data$Date))

y he echado mucho de menos Scala.

Más sobre Scala:

• Si yo fuera rey, todos los niños aprenderían Scala.
• Al tipo que inventó Scala le gustan tanto o más los oneliners que a mí.
• Todavía me llevo mal con el compilador.
• La gente viene a Suiza y aprende el alemán malhablado de aquí; yo, ya véis, Scala. Soy así de sociable.
• Detrás de Scala vendrá Spark.

Me gustaría no tener que hacer más t-tests en la vida, pero no va a ser el caso.

El problema al que me refiero le surgió a alguien en una galaxia lejana y, de alguna manera, me salpicó y me involucró. Es, simplificándolo mucho, el siguiente.

Tiene una muestra $X = x_1, \dots, x_n$ y quiere ver si la media es o no cero. ¿Solución de libro? El t-test. Pero le salen cosas raras e inesperadas. De ahí lo del salpicón.

Avisé en mi entrada del otro día: no me preguntéis por qué (imponer restricciones en un problema de mínimos cuadrados).

Pero cuanto más pienso sobre ello, menos claro lo tengo. ¿Por qué restricciones?

Primero, el contexto. O el casi contexto. Porque no es exactamente así. Pero sí parecido. Supongamos que queremos predecir algo y construimos, p.e., 4 modelos. Se nos ocurre (y hay buenas razones para ello) combinar los predictores.

Uno puede pensar en usar la media de las predicciones. O la mediana. O tratar de usar un peso revelado por los datos.

Estoy sin tiempo, así que os suelto el código y me largo a casa a no cenar. Es así:

library(parallel)

cl <- makeCluster(8)

# solo si hay aleatorización
# clusterSetRNGStream(cl, 123)

clusterEvalQ(cl,
{
	# las librerías necesarias tienen que cargarse
	# en cada esclavo
	library(rpart)

	# en la práctica, hay que cargar los datos
	# (¿desde fichero?) en cada esclavo
	my.data <- iris

	# lo mismo con las funciones necesarias
	foo <- function(x, dat){
		train <- 1:nrow(dat) %% 10 != 1
		mod <- rpart(Species ~ ., data = dat[train,])
		res <- predict(mod, dat[!train,])
	}
})

res <- parSapply(cl, 0:9,
	function(x) foo(x, my.data), simplify = F)

Sí, había restricciones. No me preguntéis por qué, pero los coeficientes tenían que ser positivos y sumar uno. Es decir, buscaba la combinación convexa de cuatro vectores que más se aproximase a y en alguna métrica razonable. Y lo resolví así:

# prepare constrained optimization

y <- dat.clean$actual
x <- t(dat.clean[,2:5])

# target function: L2 first, then other metrics

L2 <- function(coef){
  sum(abs((y - colSums(x * coef)))^1.5)
}

# restrictions: coefs > 0, sum(coefs) ~ 1

ui <- rbind(diag(4), c(-1,-1,-1,-1), c(1,1,1,1))
ci <- c(0,0,0,0,-1.000001,0.999999)

theta <- rep(0.25, 4)

best.coef <- constrOptim(theta, L2,
  grad = NULL, ui = ui, ci = ci)

coefs <- best.coef$par

Objetos aparte de x e y, hay:

A los estadísticos se nos acusa en ocasiones de contestar preguntas tontas en las que nadie está interesado.

(Nota: de alguna manera conseguí el artículo al que se refiere el enlace anterior; pero ahora no veo que exista ninguna copia libre por ahí. Si alguien la consigue, por el bien del resto de los lectores, que me avise o que lo haga saber en los comentarios).

A lo que iba. Muchos estadísticos tienen el cerebro reprogramado para tratar de no cometer los llamados errores de tipo I y errores de tipo II (y para ello tratan de estimar una cosa de dudosa utilidad, $P(D|H)$, donde $D$ son los datos y $H$ es cierta hipótesis (que, generalmente, a nadie interesa y que es más difícil de plantear correctamente de lo que parecería).