Muestreando la distribución uniforme sobre la esfera unidad en n dimensiones

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.

Para muestrear la bola n-dimensional, hay que muestrear primero la esfera (como en el párrafo anterior) y luego generar m variables aletorias $U_i$ con la distribución uniforme. La muestra en la esfera unidad será entonces $U_i^{1/n} X_i/\| X_i \|$.

Efectivamente, $X_i/\| X_i \|$ proporciona la dirección. Y en cuanto a la distancia con respecto al centro hay que tener encuenta que la bola de radio r < 1 contiene sólo un [latex]r^n[/latex] del volumen de la bola unidad. Como [latex]P( U_i < r ) = r[/latex], entonces [latex]P( U^{1/n} < r ) = P( U < r^n ) = r^n[/latex].

En R,

n <- 100
d <- 2

x <- matrix( rnorm(d * n), n, d )
x.norma <- sqrt( rowSums( x^2 ) )
u <- runif( n )

plot( x / x.norma, col = "red", xlab = "", ylab = "" )
points( x * u^{1/d} / x.norma, col = "blue" )

Y si alguien quiere ver las rodajas de una distribución uniforme sobre la esfera 5-dimensional, por ejemplo, puede ejecutar

n <- 100
d <- 5
 
x <- matrix( rnorm(d * n), n, d )
x.norma <- sqrt( rowSums( x^2 ) )
u <- runif( n )
 
plot( data.frame(x * u^{1/d}/ x.norma ))

¿Le sugiere el gráfico (p.e., los rangos de las proyecciones sobre los distintos ejes) algún comentario?