R y la distribución de Rayleigh

En la reunión de usuarios de R de Madrid de ayer, Carlos Ortega estudió la distribución en el tiempo del número de bugs que aparecen en el código de R en cada versión. Indicó que es plausible que sigan una distribución de Rayleigh, relativamente frecuente en ese tipo de contextos. E indicó que esta distribución, no tan conocida, tiene que ver (he olvidado lo que dijo exactamente) con dos normales independientes.

Efectivamente, según la Wikipedia, la distribución de Rayleigh (de parámetro \sigma)admite la caracterización

\sqrt{ X^2 + Y^2 } donde X, Y \sim N(0, \sigma).

Es decir, es el módulo de un vector bidimensional aleatorio cuyas componentes son normales con la misma varianza.

En efecto,

P( \sqrt{ X^2 + Y^2 } < a ) = \frac{1}{2\pi\sigma^2}\int_{x^2 + y^2 < a} \exp \left(-\frac{x^2 + y^2}{2\sigma^2} \right)= c \int_0^a r \exp \left( -\frac{r^2}{2\sigma^2} \right)

mediante un cambio de variables (a coordenadas polares) de los afortunados. Y la integral resultante ya puede resolverse mediante el truco, desconocido en mi época de estudiante, de copiar y pegar en Wolfran Alpha. (Los trucos de mi época de estudiante, hay que reconocerlo, eran sustancialmente más tediosos; y es que las ciencias avanzan una barbaridad). Y el resultado es (ajustando lo que hay que ajustar en la constante para que en el infinito dé lo que tiene que dar)

P( \sqrt{ X^2 + Y^2 } < a ) = 1 - \exp \left( -\frac{a^2}{2\sigma^2} \right),

nuestra distribución de Rayleigh.

La integral anterior la encontré hace muchos años, cuando tenía 19, en segundo de carrera, cuando nos explicaron —y aún lo recuerdo porque me pareció, en su día, pura taumaturgia— cómo muestrear la distribución normal.

En efecto, aunque muestrear un valor normal es difícil, un razonamiento basado en lo que aparece más arriba permite muestrear dos valores normales independientes, es decir, un punto del plano (X,Y) donde X e Y son normales: primero se obtiene el módulo de un punto en el plano usando la distribución de Rayleigh y luego el ángulo usando la distribución uniforme sobre la circunferencia.

¿Y cómo se muestrea la distribución de Rayleigh? Invirtiendo la función de distribución: dado que por definición, P(X<a)=F(a)[/latex], se tiene que [latex]P(F(X)<a)= P( X < F^{-1}(a)) = F^{-1}F(a)= a[/latex]; es decir, [latex]F(X)[/latex] es uniforme. Invirtiendo [latex]F[/latex], se puede muestrear [latex]X[/latex] haciendo   <p style="text-align: center;">[latex]X = \sigma \sqrt{ -2 \log U},

donde U es uniforme.

¿Vale esto para algo? Pues mírese el código de R y, en particular, el fichero snorm.c y uno encontrará

if(BM_norm_keep != 0.0) { /* An exact test is intentional */
	    s = BM_norm_keep;
	    BM_norm_keep = 0.0;
	    return s;
} else {
	    theta = 2 * M_PI * unif_rand();
	    R = sqrt(-2 * log(unif_rand())) + 10*DBL_MIN; /* ensure non-zero */
	    BM_norm_keep = R * sin(theta);
	    return R * cos(theta);
}

¡Y mirad lo ahorrativo que es R, que guarda el valor obtenido con el seno y con el coseno!

3 comentarios sobre “R y la distribución de Rayleigh

  1. Jorge 23 marzo, 2012 21:40

    Lo he mirado un poco por encima pero tiene pinta de ser el mismo generado de normales que se conoce como Box-Muller, no?

    http://es.wikipedia.org/wiki/M%C3%A9todo_de_Box-Muller

    Saludos y gracias por el blog.

    PD: El pasado día quise ir a la reunión pero me fue imposible estaré atento para próximas.

  2. datanalytics 23 marzo, 2012 22:04

    @Jorge Sí, sí, tienes razón, es el mismo. Lo que pasa que si uno dice «es el método de Smith y Krasiliov», la gente excusa mirar debajo del capó. Y precisamente, queria hacer lo contrario.

    Creo que también lo he visto llamar método del seno-coseno. Pero bueno.

    Y respecto a la reunión de usuarios… ¡tranquilo que habrá otra(s)!

  3. Jorge 23 marzo, 2012 22:37

    Tienes razón, me ha gustado más esta explicación que la que conocía. Y además una distribución nueva al cajón.

    Saludos.

Los comentarios están desabilitados.