Muestreo de una distribución (a partir de un modelo de galaxia)

El artículo se ve realmente mal. De hecho, hay dos errores.

Primero, tiene razón en que el método de aceptación-rechazo debe aplicarse a$\rho(r)$, y no a$m(r)$. Para entender cómo funciona esta idea, supongamos que queremos generar una función de distribución normalizada unidimensional$p(y)$. Ahora, supongamos que podemos reescribir esta función de distribución en términos de una variable$x$, tal que toma la forma de una distribución uniforme. Eso es,

$$p(x) = \begin{cases} 1& \text{for $0\leqslant x \leqslant 1,$}\\ 0& \text{elsewhere}. \end{cases}$$ Dado $p(y)$, qué es $x$? Tenemos la transformación jacobiana $$p(y)dy = p(x(y))\left|\frac{dx}{dy}\right|dy = \left|\frac{dx}{dy}\right|dy,$$ lo que implica $$p(y) = \frac{dx}{dy},$$ asumiendo que $x(y)$es una función creciente. De este modo $$x = \int_0^y p(y')dy' = F(y).$$ En otras palabras, la integral de $p(y)$(o de manera equivalente, el área bajo la curva) sigue una distribución uniforme. Con esto en mente, existen esencialmente dos formas de realizar una simulación Monte-Carlo.

La primera forma es el método de aceptación-rechazo: trace la curva$p(y)$y generar uniformemente un par de números$(a,b)$en el intervalo$([0,y_\max],[0,p_\max])$, dónde$y_\max$y$p_\max$son los límites superiores de$y$y$p(y)$. Si la coordenada$(a,b)$se encuentra debajo de la curva$p(y)$, acéptalo; de lo contrario, rechácelo. Si se acepta la coordenada,$y=a$es el punto generado.

Hay grandes inconvenientes en este método:$y_\max$y$p_\max$puede ser infinito, por lo que uno necesitaría un corte. Y si$p(y)$tiene un pico agudo, uno termina rechazando muchos puntos.

Un método mucho más eficiente es generar uniformemente$x$, y calcular el correspondiente$y$invirtiendo$x=F(y)$:

$$y = F^{-1}(x).$$ Esto llena automáticamente el área bajo la curva, sin rechazar puntos.

Si el cálculo de$F^{-1}(x)$está demasiado involucrado numéricamente, se puede usar una combinación de ambos métodos: introducir otra función (más simple)$f(y)$que yace en todas partes arriba$p(y)$. Aplicar el método de inversión a$f(y)$, generando un punto$y$. Luego genera uniformemente un valor$b$en el intervalo$[0,f(y)]$. Si$b\leqslant p(y)$, aceptar$y$; de lo contrario, rechácelo.

Ahora, considere la distribución de Hernquist. Como tiene una cúspide en el origen y la masa acumulada$m(r)$es una funcion sencilla

$$m(r) = M\frac{r^2}{(a+r)^2},$$ Definitivamente recomendaría el método de inversión. Pero hay una advertencia importante aquí, y ese es el segundo error del artículo:$\rho(r)$no es realmente una distribución unidimensional . En cambio, es una distribución en un espacio tridimensional, y es solo una función de una variable debido a la simetría esférica. Para aplicar el método Monte-Carlo, tenemos que expresar $\rho$como una función de distribución verdaderamente unidimensional, lo que podemos hacer expresándola en términos del volumen $$y = \frac{4\pi}{3}r^3.$$ Ahora tenemos $$p(y) = \rho(y) = \frac{M}{2\pi}\frac{a\,(3y/4\pi)^{-1/3}}{\left[a + (3y/4\pi)^{1/3}\right]^3},\\ F(y) = m(y) = \int_0^y\rho(y')dy' = M\frac{(3y/4\pi)^{2/3}}{\left[a + (3y/4\pi)^{1/3}\right]^2}.$$ Una vez que generamos un punto $y$, el radio correspondiente es $$r=\left(\frac{3y}{4\pi}\right)^{1/3}.$$

Hay una consecuencia importante: es probable que haya más partículas en radios grandes que alrededor del centro, aunque$\rho(r)$es mucho mayor en radios pequeños. La razón es que las partículas entre dos radios$r$y$(r+\Delta r)$ocupan una concha con volumen

$$V = \frac{4\pi}{3}\left[(r+\Delta r)^3-r^3\right].$$ Cuanto mayor sea el radio $r$, cuanto mayor sea el volumen del caparazón, lo que significa que necesita más partículas para llenarlo y obtener $\rho(r)$. Esto es obvio en el caso de una densidad constante, pero también es cierto para las densidades generales.

la FCD,$F(x)$, está relacionado con el PDF,$f(x)$, a través de la relación:

$$F(x) = \int_{-\infty}^xdx'\,f(x')$$ En el caso de distribuciones radiales, su límite inferior es obviamente 0 y no $-\infty$. Por lo tanto, su CDF es $m(r)$y el pdf es $4\pi(r')^2\rho(r')$(técnicamente debería ser $\rho(r')$con el $4\pi(r')^2$procedente de $dx\to dr$e isotropía espacial, pero lo que sea).

En el caso de funciones invertibles (por ejemplo,$f(x)=Ax$con constante de normalización$A$), puedes resolver esto como

$$F(x)=\frac{A}{2}x^2\to x=\sqrt{\frac{2F}{A}}$$ Si genera un número aleatorio, configúrelo igual a $F$y sale tu $x$que satisface el PDF.

En el caso de que sus funciones no sean invertibles (a menudo el caso que involucra distribuciones radiales), sugeriría usar un método iterativo de Newton-Raphson para encontrar raíces para este caso. Esto se puede hacer fácilmente ya que sabes$F(x)$(su FCD)$F'(x)=f(x)$(tu PDF):

$$x_{new}=x_{old}-\frac{F(x_{old})}{f(x_{old})}$$ Cuando $|x_{new}|<\epsilon$, entonces $x_{old}$es la raíz (aproximada). A menudo, solo se necesitan unas pocas iteraciones para la convergencia.

---

Aparte, no puedo enfatizar lo suficiente lo terrible que es la sugerencia del método de aceptación-rechazo. NO USE ESTE MÉTODO . Literalmente desperdicia tiempo de computación, algo que debería considerarse valioso, independientemente de la prevalencia de las computadoras. No escuches a nadie que diga que tienes que usar este método, están total y absolutamente equivocados.

El método de Newton que sugiero anteriormente no rechaza un solo número aleatorio, se ajustará al PDF con precisión y es altamente eficiente.