c++ - tipos - ejercicios resueltos de distribucion uniforme continua pdf

Puntos aleatorios de distribución uniforme y rápida en la superficie de un hemisferio unitario

Question

Estoy intentando generar puntos aleatorios uniformes en la superficie de una esfera unitaria para un programa de trazado de rayos de Monte Carlo. Cuando digo uniforme me refiero a que los puntos se distribuyen uniformemente con respecto al área de superficie. Mi metodología actual es calcular puntos aleatorios uniformes en un hemisferio que apunta en el eje z positivo y la base en el plano xy.

El punto aleatorio en el hemisferio representa la dirección de emisión de la radiación térmica para un emisor gris difuso.

Logro el resultado correcto cuando uso el siguiente cálculo:

Nota: dsfmt * is devolverá un número aleatorio entre 0 y 1.

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
zenith = asin(sqrt(dsfmt_genrand_close_open(&dsfmtt)));

// Calculate the cartesian point
osRay.c._x = sin(zenith)*cos(azimuthal); 
osRay.c._y = sin(zenith)*sin(azimuthal);
osRay.c._z = cos(zenith);

Sin embargo, esto es bastante lento y la creación de perfiles sugiere que ocupa una gran parte del tiempo de ejecución. Por eso busqué algunos métodos alternativos:

El método de rechazo de Marsaglia 1972

do {
   x1 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;
   x2 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;
   S = x1*x1 + x2*x2;
} while(S > 1.0f);


osRay.c._x = 2.0*x1*sqrt(1.0-S);
osRay.c._y = 2.0*x2*sqrt(1.0-S);
osRay.c._z = abs(1.0-2.0*S);

Cálculo analítico de coordenadas cartesianas.

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
u = 2*dsfmt_genrand_close_open(&dsfmtt) -1;
w = sqrt(1-u*u);

osRay.c._x = w*cos(azimuthal);
osRay.c._y = w*sin(azimuthal);
osRay.c._z = abs(u);

Mientras estos dos últimos métodos se ejecutan a veces más rápido que el primero, cuando los uso obtengo resultados que indican que no están generando puntos aleatorios uniformes en la superficie de una esfera, sino que dan una distribución que favorece al ecuador.

Además, los dos últimos métodos dan resultados finales idénticos; sin embargo, estoy seguro de que son incorrectos, ya que estoy comparando con una solución analítica.

Todas las referencias que he encontrado indican que estos métodos producen distribuciones uniformes, sin embargo no logro el resultado correcto.

¿Hay algún error en mi implementación o me he perdido una idea fundamental en el segundo y tercer método?

Answer 1

¿Has intentado deshacerte de asin ?

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
sin2_zenith = dsfmt_genrand_close_open(&dsfmtt);
sin_zenith = sqrt(sin2_zenith);

// Calculate the cartesian point
osRay.c._x = sin_zenith*cos(azimuthal); 
osRay.c._y = sin_zenith*sin(azimuthal);
osRay.c._z = sqrt(1 - sin2_zenith);

Answer 2

1er intento (mal)

point=[rand(-1,1),rand(-1,1),rand(-1,1)];
len = length_of_vector(point);

EDITADO:

¿Qué pasa?

while(1)
 point=[rand(-1,1),rand(-1,1),rand(-1,1)];
 len = length_of_vector(point);
 if( len > 1 )
     continue;
 point = point / len
     break

La aceptación es aquí aprox. 0.4. Eso significa que rechazarás el 60% de las soluciones.

Answer 3

Creo que el problema que tiene con los resultados no uniformes es que en las coordenadas polares, un punto aleatorio en el círculo no está distribuido uniformemente en el eje radial. Si observa el área en [theta, theta+dtheta]x[r,r+dr] , para theta y dtheta , el área será diferente de los diferentes valores de r . Intuitivamente, hay "más área" más alejada del centro. Por lo tanto, necesitas escalar tu radio aleatorio para tener en cuenta esto. No tengo la prueba por ahí, pero la escala es r=R*sqrt(rand) , siendo R el radio del círculo y rand comienzan el número aleatorio.

Answer 4

De hecho, el segundo y tercer método producen puntos aleatorios distribuidos uniformemente en la superficie de una esfera, y el segundo método ( Marsaglia 1972 ) produce los tiempos de ejecución más rápidos con aproximadamente el doble de la velocidad en un Intel Xeon 2.8 GHz Quad-Core.

Como señaló hay un método adicional que utiliza la distribución normal que se expande a n-esferas mejor que los métodos que he presentado.

Este enlace le dará más información sobre cómo seleccionar puntos aleatorios distribuidos uniformemente en la superficie de una esfera.

Mi método inicial, según lo señalado por TonyK, no produce puntos distribuidos uniformemente y más bien sesga los polos al generar los puntos aleatorios. Esto es requerido por el problema que estoy tratando de resolver, pero simplemente asumí que generaría puntos aleatoriamente uniformes. Según lo sugerido por este método se puede optimizar eliminando la llamada asin () para reducir el tiempo de ejecución en alrededor del 20%.

Answer 5

Esto debería ser rápido si tienes un RNG rápido:

// RNG::draw() returns a uniformly distributed number between -1 and 1.

void drawSphereSurface(RNG& rng, double& x1, double& x2, double& x3)
{
    while (true) {
        x1 = rng.draw();
        x2 = rng.draw();
        x3 = rng.draw();
        const double radius = sqrt(x1*x1 + x2*x2 + x3*x3);
        if (radius > 0 && radius < 1) {
            x1 /= radius;
            x2 /= radius;
            x3 /= radius;
            return;
        }
    }   
}

Para acelerarlo, puede mover la llamada sqrt dentro del bloque if .

Answer 6

La forma más sencilla de generar una distribución uniforme en la esfera unitaria (cualquiera que sea su dimensión) es dibujar distribuciones normales independientes y normalizar el vector resultante.

De hecho, por ejemplo en dimensión 3, e ^ (- x ^ 2/2) e ^ (- y ^ 2/2) e ^ (- z ^ 2/2) = e ^ (- (x ^ 2 + y ^ 2 + z ^ 2) / 2) por lo que la distribución conjunta es invariante por las rotaciones.

Esto es rápido si utiliza un generador de distribución normal rápido (Ziggurat o Ratio-Of-Uniforms) y una rutina de normalización rápida (google para "raíz cuadrada inversa rápida). No se requiere una llamada de función trascendental.

Además, la Marsaglia no es uniforme en la media esfera. Tendrá más puntos cerca del ecuador ya que el punto de correspondencia en el punto 2D del disco 2D en la media esfera no es isométrico. Sin embargo, el último parece correcto (sin embargo, no hice el cálculo para asegurar esto).

Answer 7

Si toma un corte horizontal de la esfera unitaria, de altura h , su área de superficie es de solo 2 pi h . (Así es como Arquímedes calculó el área de superficie de una esfera). Por lo tanto, la coordenada z se distribuye uniformemente en [0,1] :

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
osRay.c._z = dsfmt_genrand_close_open(&dsfmtt);

xyproj = sqrt(1 - osRay.c._z*osRay.c._z);
osRay.c._x = xyproj*cos(azimuthal); 
osRay.c._y = xyproj*sin(azimuthal);

También es posible que pueda ahorrar algo de tiempo calculando cos(azimuthal) y sin(azimuthal) juntos: vea esta pregunta de para discusión.

Editado para agregar: OK, ahora veo que esto es solo una pequeña modificación de su tercer método. Pero se corta un paso.