Me gustaría generar muestras de la región azul definida aquí:
La solución ingenua es utilizar el muestreo de rechazo en el cuadrado de la unidad, pero esto proporciona solo una eficiencia de (~ 21.4%).
¿Hay alguna manera de que pueda probar más eficientemente?
probability
sampling
monte-carlo
random-generation
Cam.Davidson.Pilon
fuente
fuente
Respuestas:
¿Serán dos millones de puntos por segundo?
La distribución es simétrica: solo necesitamos calcular la distribución para un octavo del círculo completo y luego copiarla alrededor de los otros octantes. En coordenadas polares , la distribución acumulativa del ángulo Θ para la ubicación aleatoria ( X , Y ) en el valor θ viene dada por el área entre el triángulo ( 0 , 0 ) , ( 1 , 0 ) , ( 1 , tan θ ) y el arco del círculo que se extiende desde(r,θ) Θ (X,Y) θ (0,0),(1,0),(1,tanθ) (1,0) a . Por lo tanto, es proporcional a(cosθ,sinθ)
de donde es su densidad
Podemos tomar muestras de esta densidad usando, por ejemplo, un método de rechazo (que tiene una eficiencia de ).8/π−2≈54.6479%
La densidad condicional de la coordenada radial es proporcional a r d r entre r = 1 y r = sec θ . Eso se puede muestrear con una fácil inversión del CDF.R rdr r=1 r=secθ
Si generamos muestras independientes , la conversión de nuevo a coordenadas cartesianas ( x i , y i ) muestrea este octante. Debido a que las muestras son independientes, el intercambio aleatorio de las coordenadas produce una muestra aleatoria independiente del primer cuadrante, según se desee. (Los intercambios aleatorios requieren generar solo una variable binomial única para determinar cuántas de las realizaciones intercambiar).(ri,θi) (xi,yi)
Cada realización de requiere, en promedio, una variante uniforme (para R ) más 1 / ( 8 π - 2 ) veces dos variables uniformes (para Θ ) y una pequeña cantidad de cálculo (rápido). Eso es 4 / ( π - 4 ) ≈ 4.66 variantes por punto (que, por supuesto, tiene dos coordenadas). Los detalles completos se encuentran en el siguiente ejemplo de código. Esta cifra representa 10,000 de más de medio millón de puntos generados.(X,Y) R 1/(8π−2) Θ 4/(π−4)≈4.66
Aquí está el
R
código que produjo esta simulación y la cronometró.fuente
Propongo la siguiente solución, que debería ser más simple, más eficiente y / o computacionalmente más barata que otras soluciones de @cardinal, @whuber y @ stephan-kolassa hasta ahora.
Implica los siguientes pasos simples:
La intuición detrás de este algoritmo se muestra en la figura.
Los pasos 2a y 2b se pueden combinar en un solo paso:
El siguiente código implementa el algoritmo anterior (y lo prueba usando el código de @ whuber).
Algunas pruebas rápidas arrojan los siguientes resultados.
Algoritmo /stats//a/258349 . Lo mejor de 3: 0,33 segundos por millón de puntos.
Este algoritmo Lo mejor de 3: 0.18 segundos por millón de puntos.
fuente
Bueno, se puede hacer de manera más eficiente , pero espero que no estés buscando más rápido .
Wolfram te ayuda a integrar eso :
Probablemente pueda acelerar bastante la inversión de CDF si invierte algo de pensamiento. Por otra parte, pensar duele. Yo personalmente elegiría muestras de rechazo, que es más rápido y mucho menos propenso a errores, a menos que tenga muy buenas razones para no hacerlo.
fuente