Simulación

Me entretuve el otro día en cómo muestrear uniformemente dentro de triángulos motivado por Randomly selecting points inside a triangle de John D. Cook.

Hay uno que se le ocurriría a cualquiera: el del rechazo. Se inserta el triángulo en un cuadrado y se seleccionan solo aquellos valores que caigan dentro del triángulo.

Hay otro, que no está en esa entrada, y que consiste en transformar el triángulo en un triángulo rectángulo mediante una transformación lineal que preserve el área (shear o cizallamiento), del tipo

Un viejo amigo me escribe y me propone (simplificándolo) el siguiente problema:

• Tengo una normal de parámetros desconocidos.
• De ella solo dispongo de un histograma.
• ¿Cómo puedo reconstruir la normal original? Es decir, ¿cómo puedo estimar los $\mu$ y $\sigma$ originales?

En el caso concreto, la normal tiene una media próxima a 255 y los valores del histograma proceden de una muestra suya redondeada al entero más próximo.

Aquí va mi solución.

El otro día publiqué un pequeño fragmento de código,

a <- 2.89
b <- 36.81
sample_dist <- function() rbeta(1, a, b)

sample_p <- function(y){
  candidate <- sample_dist()
  my_sample <- runif(1)

  if (y == 1) if (my_sample < candidate) return(candidate)
  if (y == 0) if (my_sample < 1 - candidate) return(candidate)

  sample_p(y)
}

p1 <- replicate(100000, sample_p(1))
p0 <- replicate(100000, sample_p(0))

auc <- mean(p1 > p0)
auc

que había usado antes aquí, para muestrear unas distribuciones relacionadas con el cálculo del AUC en modelos perfectamente calibrados. Lo había escrito meses atrás y supongo que me pasó como a la mayoría de mis lectores: darlo por bueno primero y usarlo después suponía todo un acto de fe (en mí, además).

[Menos mal que se me ha ocurrido buscar en mi propio blog sobre el asunto y descubrir —no lo recordaba— que ya había tratado el asunto previamente en entradas como esta, esta o esta.]

El problema de la propagación de errores lo cuentan muy bien Iñaki Úcar y sus coautores aquí. Por resumirlo: tienes una cantidad, $X$ conocida solo aproximadamente —en concreto, con cierto error— e interesa conocer y acotar el error de una expresión $f(X)$.

Está extraído de aquí y dice los siguiente:

Las Proyecciones de Población constituyen una simulación estadística de la población que residiría en España, sus comunidades autónomas y provincias en los próximos años, así como de la evolución de cada uno de los fenómenos demográficos básicos asociados, en caso de mantenerse las tendencias y comportamientos demográficos actualmente observados.

Para interpretar correctamente los resultados de las Proyecciones de Población es importante distinguir entre previsiones y proyecciones demográficas. Si bien pueden emplear el mismo método de cálculo, difieren en la filosofía.

Esta entrada es casi una referencia para mí. Cada vez tiro más de lme4 en mis modelos y en uno en concreto que tengo entre manos toca simular escenarios. Para lo cual, simulate.merMod.

Véamoslo en funcionamiento. Primero, datos (ANOVA-style) y el modelo que piden a gritos:

library(plyr)
library(lme4)

a <- c(0,0,0, -1, -1, 1, 1, -2, 2)
factors <- letters[1:length(a)]

datos <- ldply(1:100, function(i){
    data.frame(x = factors, y = a + rnorm(length(a)))
})
modelo <- lmer(y ~ (1 | x), data = datos)

El resumen del modelo está niquelado:

summary(modelo)

# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ (1 | x)
# Data: datos
#
# REML criterion at convergence: 2560.3
#
# Scaled residuals:
#     Min      1Q  Median      3Q     Max
# -3.6798 -0.6442 -0.0288  0.6446  3.3582
#
# Random effects:
#     Groups   Name        Variance Std.Dev.
# x        (Intercept) 1.5197   1.2328
# Residual             0.9582   0.9789
# Number of obs: 900, groups:  x, 9
#
# Fixed effects:
#     Estimate Std. Error t value
# (Intercept) -0.009334   0.412212  -0.023

En particular,

Por referencia y afán de completar dos entradas que hice hace un tiempo sobre el método delta, esta y esta, dejo constar mención al paquete propagate, que contiene métodos para la propagación de la incertidumbre.

Para desavisados: si $x \sim N(5,1)$ e $y \sim N(10,1)$, ¿cómo sería la distribución de $x/y$? Etc.

NIMBLE ha sido uno de mis más recientes y provechosos descubrimientos. Mejor que hablar de él, que otros lo harán mejor y con más criterio que yo, lo usaré para replantear el problema asociado el método delta que me ocupó el otro día.

Casi autoexplicativo:

library(nimble)

src <- nimbleCode({
    T_half <- log(.5) / k
    k ~ dnorm(-0.035, sd = 0.00195)
})

mcmc.out <- nimbleMCMC(code = src,
    constants = list(),
    data = list(), inits = list(k = -0.035),
    niter = 10000,
    monitors = c("k", "T_half"))

out <- as.data.frame(mcmc.out)

# hist(out$T_half), sd(out$T_half), etc.

Cosas:

Tenemos en Circiter un proyecto sobre el que no puedo dar muchos detalles, pero que vamos a plantear (en versión muy resumida) como un modelo que alimenta una simulación.

El modelo no va a ser un modelo sino un modelo por sujeto (rebaños, los llamamos aquí). Los modelos serán, casi seguro, modelos mixtos (lmer/glmer).

Pero claro, si usas un modelo, por muy mixto que sea, con intención de simular, predict se queda muy corto (¡siempre da la el mismo resultado!).

Debo esta entrada a la diligencia de Juanjo Gibaja, que se tomó la molestia de ubicar los teoremas relevantes en el libro Simulation and the Monte Carlo Method de Rubinstein y Kroese.

Esencialmente, como la distribución normal multivariante (con matriz de covarianzas I) es simétrica, entonces, dadas $X_1,\dots, X_m \sim N( 0, I_n )$ independientes, los m puntos del espacion n-dimensional $X_i/| X_i |$ siguen una distribución uniforme sobre su esfera (su superficie, vale la pena reiterar) unidad.