Rango Finito

fotoscódigoobservatorioshermanocerdo temas plots

probabilidades

Relevos

¿Cada cuántos días se muere la persona más vieja del mundo? ¿Cómo estimar la probabilidad de que la persona más vieja del mundo muera hoy si se sabe apenas que es la persona más vieja del mundo y se conoce desde hace cuánto lo es? ¿Cuál sería un modelo razonable para describir este proceso de relevos que marca, podría decirse, el lento desprendimiento colectivo del pasado humano vivo? ¿Cuál es la probabilidad (actual) de que las N (N ≥ 2) personas más viejas del mundo hayan nacido en estrictamente distintos días? ¿Cuántos días en la historia cuentan ahora mismo con exactamente una persona viva nacida en ese día (calculándolos, digamos, como ventanas de 24 horas alrededor del nacimiento de una persona donde no queda nadie más nacido en esas horas)?

Un clasificador de Bayes falso

Lo que sigue es un pequeño complemento a Tres Modelos que estuve conversando con Santiago hoy: una forma de usar Bayes (asignar la clasificación de acuerdo a qué color tiene la mayor probabilidad en cada punto) aunque no se cuente con la distribución de densidad de las poblaciones que se estudian es asumir alguna distribución de densidad a priori basada en datos empíricos y proceder desde ahí. Hay una discusión filosófica al fondo de eso pero no importa la posición que se tenga el método sigue teniendo sentido. Sobre esto me gustaría escribir algo detallado. Probablemente le dedique otro texto largo pues es un criterio de selección de procedimientos clave en teoría de decisión estadística.

Lo que quiero hacer ahora es un poco distinto y probablemente bastante sucio: con los datos de entrenamiento a mano (¡y solo con ellos!) hay métodos que permiten calcular una función de densidad falsa que se conoce como kernel density estimate (KDE). Para calcular el KDE, además de los datos $x_1, x_2,\ldots,x_m\in\mathbb{R}^n$ se necesita un $h>0$ y una función $$K:\mathbb{R}^n\to \mathbb{R}$$ simétrica con respecto al origen pero no necesariamente positiva con la condición de que $$\int_{\mathbb{R}^n} K = 1.$$ El $K$ estándar es una multinormal centrada en cero y con la matriz identidad como covarianza.

Dados estos ingredientes, a el KDE es una función $f_{K,h}$ definida como sigue: $$f_{K,h}(x)=\frac{1}{mh}\sum_{i=1}^m K\left(\frac{x-x_i}{h}\right).$$

O sea una suerte de promedio local ponderado de las distancias del punto a los datos disponibles.

En el código que sigue tomo la muestra de entrenamiento y para cada color calculo un KDE con $K$ normal bivariada y $h$ elegido con un método estándar para estos menesteres que está descrito acá. Después grafico los KDE en tres dimensiones para ver las dos montañas:


# Primero extraigo las poblaciones de cada color:
blues <- training_sample[training_sample$color == "blue", 1:2]
oranges <- training_sample[training_sample$color == "orange", 1:2]

# Después calculo los KDE en 2D:
dblues <- kde2d(blues[,1], blues[,2], n=200, lims = c(-3, 3, -3, 3))
doranges <- kde2d(oranges[,1], oranges[,2], n=200, lims = c(-3, 3, -3, 3))

# Y finalmente lo grafico:
persp3D(x=doranges$x, y=doranges$y, z= doranges$z, shade=.5, col="orange", 
    alpha= .5, phi=20, box=F, contour=T, theta=-30)
persp3D(x=dblues$x, y=dblues$y, z= dblues$z, shade=.5, col="dodgerblue2", 
    alpha= .5, add=T, phi=20, contour=T,  theta=-30)

Densidades Falsas

Para comparar, aquí las densidades reales de las poblaciones:

Densidad Real

Ahora lo que se puede hacer, por jugar, es proponer una clasificación basada en estas densidades falsas y mirar cómo se diferencia de la clasificación de Bayes. Grafiquemos ambas fronteras, Bayes (rojo) y Bayes Falsa (verde), sobre la muestra de evaluación:


# Primero organizo los KDE en un dataframe:
density.df <- expand.grid(dblues$x, dblues$y)
density.df$blues <- as.vector(dblues$z)
density.df$oranges <- as.vector(doranges$z)

# Ahora defino el Bayes falso (1 es naranja y 0 es azul):
density.df$fake.bayes <- as.numeric(density.df$blues < density.df$oranges)

# Finalmente lo grafico:
grid$color <- apply(grid[,1:2], 1, bayes_classifier)
ggplot(grid, aes(X1, X2, z=color)) +  
    geom_point(aes(X1, X2, fill=as.factor(color)), 
        size=1, col="white", shape=21, alpha=0.5) + 
    geom_point(data=test_sample, aes(x=X1, y=X2, col=color)) +
    stat_contour(bins=1, color="red", size=2) +
    stat_contour(data=density.df, aes(x=Var1, y=Var2, z=fake.bayes), 
        color="forestgreen", bins=1, size=2) +
    scale_color_manual(guide="none", values=c("dodgerblue2", "orange")) +
    scale_fill_manual(guide="none", values=c("dodgerblue2", "orange")) +
    theme_bw() + xlab("x") + ylab("y")

bayes.vs.fake.bayes

Obviamente el Bayes falso es muchísimo más débil que Bayes: se deja manipular demasiado por la muestra de entrenamiento.

Ejercicio: calcular el error de este modelo de Bayes falso en la muestra de evaluación.

Prevenir lo inevitable

La columna de hoy habla sobre lo que los sismólogos pueden hacer y lo que no pueden hacer. Resumiendo: Los sismólogos no pueden predecir terremotos pero sí pueden estimar (muy burdamente) una probabilidad de que un terremoto ocurra en un cierto lugar, así como su rango de magnitudes posibles y otros detalles así. Por eso es que aunque los sismólogos italianos no podían predecir el terremoto de L’Aquila, estaba cantado que La Aquina podía recibir un susto como ese cualquier día. La falta de prevención (y no de predicción) fue lo que produjo la tragedia. Quienes dicen que los científicos fueron condenados por entregar información errada se equivocan: fueron condenados por no decir algo que no tenían cómo decir (i.e., que el terremoto venía.) En Colombia el caso de Tumaco es emblemático tanto por el tamaño legendario del terremoto de 1906 como por la situación de abandono estatal de la ciudad. En un artículo para Universo Centro ahondaré más en la discusión sobre las capacidades de los sismólogos para predecir y algunos detalles relevantes del caso italiano. Por lo pronto dejo unos cuantos enlaces.

Sobre Tumaco

  • Aquí se puede bajar el plan de contingencia contra tsnunamis de la ciudad de Tumaco. Incluye párrafos como el siguiente:
    El sistema principal de alarma para tsunami de origen cercano será el acordado por el CLOPAD en 2002, es decir el sistema de alarma personal. Este sistema consiste en la decisión personal de cada habitante de dirigirse a la zona segura más cercana cuando sienta un sismo que le dificulte mantenerse en pie.

    Por si las moscas, también hay (o al menos fueron recomendadas) once sirenas.

  • Aquí una foto de Tumaco desde un satélite para entender mejor la gravedad de la vaina:
  • Este es un artículo ya clásico donde se estudia la relación entre el terremoto de 1906 y los otros tres terremotos posteriores. El artículo sugiere que los tres posteriores no son suficientes para descargar la energía del de 1906 y por tanto todavía hay un golpecito pendiente por ahí más bien pronto.
  • Y este artículo, más reciente, usa modelos computacionales y concluye lo contrario: los tres terremotos posteriores sí parecen haber quemado otro cartucho parecido al de 1906, pero en cuotas más “amables” (aunque “amable” sea un decir, porque en el de 1979 murieron como 500 personas en Tumaco.)
  • Convocatoria al último simulacro general de evacuación de la ciudad, en 2010.

Sobre el caso de L’Aquila

Víctimas

La Operación Doorstep, una explosión nuclear de dieciséis kilotones en el desierto de Nevada en marzo de 1953, estudiaba, entre otras cosas, las probabilidades de supervivencia de la familia americana tradicional bajo un ataque nuclear dados los medios a su disposición. Como parte de la operación, se plantaron maniquís en el área de impacto representando escenas cotidianas de la época. El informe instructivo de la defensa civil (PDF) (que recomienda la construcción de los populares refugios antiatómicos en los patios de las casas como única opción segura) incluye fotos de los maniquís antes y después de la explosión. Son tristísimas.