6.2 El método de aceptación/rechazo

El algoritmo de aceptación-rechazo es el mismo que el del caso univariante descrito en la Sección 4.2, la única diferencia es que las densidades son multidimensionales. Supongamos que la densidad objetivo \(f\) y la densidad auxiliar \(g\) verifican: \[f\left( x_1,x_2,\ldots,x_d\right) \leq c\cdot g\left( x_1,x_2,\ldots,x_d\right) \text{, }\forall \mathbf{x} = \left( x_1,x_2,\ldots,x_d\right)\in \mathbb{R}^d\text{.}\] para una constante \(c>0\). El algoritmo sería:

  1. Generar \(U\sim \mathcal{U}\left( 0,1\right)\).

  2. Generar \(\mathbf{T} = \left( T_1,T_2,\ldots,T_d\right) \sim g\).

  3. Si \(c\cdot U\cdot g\left( T_1,T_2,\ldots,T_d\right) \leq f\left( T_1,T_2,\ldots,T_d\right)\) devolver \(\mathbf{X}=\mathbf{T}\).

    En caso contrario volver al paso 1.

Por ejemplo, de forma análoga al caso unidimensional, en el caso de una densidad acotada en un hipercubo (intervalo cerrado multidimensional) siempre podríamos considerar una uniforme como densidad auxiliar.

Ejemplo 6.2 (distribución bidimensional acotada)

Supongamos que estamos interesados en generar valores de una variable aleatoria bidimensional \(\left( X,Y\right)\) con función de densidad: \[f(x,y)=\left\{ \begin{array}{cl} \frac{3}{16}\left( 2-\left( x^2+y^2\right) \right) & \text{si }x\in \lbrack -1,1]\text{ e }y\in \lbrack -1,1] \\ 0 & \text{en otro caso} \end{array} \right.\]

Podríamos considerar como densidad auxiliar la uniforme en \(\left[ -1,1\right] \times\left[ -1,1\right]\):

\[g\left( x, y \right) =\left\{ \begin{array}{ll} \frac{1}{4} & \text{si }x\in \lbrack -1,1]\text{ e }y\in \lbrack -1,1] \\ 0 & \text{en otro caso} \end{array}\right.\]

Como \(f(x, y) \leq M = f(0,0) = \frac38\), tomando \(c=\frac{M}{g(x,y)} = \frac32\) tendríamos que \(f(x,y) \leq cg(x,y) = M\) y el algoritmo sería:

  1. Generar \(U \sim \mathcal{U}\left( 0, 1\right)\).

  2. Generar \(T_1, T_2 \sim \mathcal{U}\left( -1, 1 \right)\).

  3. Si \(M \cdot U\leq f\left( T_1, T_2 \right)\) devolver \(\mathbf{X} = \left( T_1, T_2 \right)\).

  4. En caso contrario volver al paso 1.

En este caso, la condición de aceptación del paso 3 simplificada sería: \(U \leq 1 - \left( T_1^2 + T_2^2 \right) / 2\).

Ejercicio 6.1 (distribución bidimensional acotada)

Escribir una función que implemente el algoritmo del ejemplo anterior, emplearla para generar 1000 observaciones y representar gráficamente su distribución (puede resultar de interés la función sm::sm.density() o las herramientas del paquete ks).

Ejemplo 6.3 (distribución uniforme en la esfera)

Supongamos que el objetivo es simular puntos uniformemente distribuidos sobre la “esfera” unitaria \(d\)-dimensional (ver Figura D.1): \[C_d=\left\{ \left( x_1, x_2, \ldots, x_d \right) \in \mathbb{R}^d : x_1^2 + x_2^2 + \cdots + x_d^2 \leq1 \right\}.\]

Denotando por \(V_d\left( 1\right)\), el “volumen” (la medida) de la esfera \(d\)-dimensional de radio \(1\) (en general, la de radio \(r\) verifica \(V_d\left( r\right) =r^{d}V_d\left( 1\right)\)), se tiene: \[f\left( x_1,x_2,\ldots,x_d\right) =\left\{ \begin{array}{ll} \frac{1}{V_d\left( 1\right) } & \text{si } \left( x_1, x_2, \ldots ,x_d\right) \in C_d\\ 0 & \text{si } \left( x_1,x_2,\ldots,x_d\right) \notin C_d \end{array} \right.\]

Para simular valores en \(\mathbb{R}^{d}\), con densidad \(f\), podemos utilizar como distribución auxiliar una \(\mathcal{U}\left( \left[ -1,1\right] \times\left[ -1,1\right] \times\overset{\text{d}}{\cdots}\times\left[ -1,1\right] \right) = \mathcal{U}\left( \left[ -1,1\right]^{d}\right)\), dada por: \[g\left( x_1,x_2,\ldots,x_d\right) =\left\{ \begin{array}{ll} \frac{1}{2^{d}} & \text{si } x_i\in\left[ -1,1\right], \text{ para todo } i=1,2,\ldots,d\\ 0 & \text{en otro caso} \end{array}\right.\]

La constante \(c\) óptima para la utilización del método de aceptación/rechazo es: \[c=\max_{\{\mathbf{x}:g\left( \mathbf{x}\right) > 0\}} \frac{f\left( \mathbf{x}\right) }{g\left( \mathbf{x}\right) } =\frac{\frac{1}{V_d\left( 1\right) }}{\frac{1}{2^{d}}} =\frac{2^{d}}{V_d\left( 1\right)}\] y la condición de aceptación \(cUg\left( \mathbf{T}\right) \leq f\left( \mathbf{T}\right)\) se convierte en: \[\frac{2^{d}}{V_d\left( 1\right) }U\frac{1}{2^{d}}1_{\left[ -1,1\right] ^{d}}\left( \mathbf{T}\right) \leq\frac{1}{V_d\left( 1\right) }1_{C_d}\left( \mathbf{T}\right),\] o, lo que es lo mismo, \(U1_{\left[ -1,1\right]^{d}}\left( \mathbf {T}\right) \leq1_{C_d}\left( \mathbf{T}\right)\). Dado que el número aleatorio \(U\) está en el intervalo \((0,1)\) y que las funciones indicadoras valen \(0\) ó \(1\), esta condición equivale a que \(1_{\left[ -1,1\right] ^{d}}\left( \mathbf{T}\right) =1_{C_d}\left( \mathbf{T}\right)\), es decir, a que \(\mathbf{T}\in C_d\), por tanto, a que se verifique: \[T_1^2+T_2^2+\cdots+T_d^2\leq1.\]

Por otra parte, la simulación de \(T \sim \mathcal{U}\left( \left[ -1,1\right] ^{d}\right)\) puede hacerse trivialmente mediante \(T_i \sim \mathcal{U}\left( -1, 1 \right)\) para cada \(i=1,2,\ldots,d\), ya que las componentes son independientes. Como el valor de \(U\) es superfluo en este caso, el algoritmo queda:

  1. Simular \(V_1,V_2,\ldots,V_d \sim \mathcal{U}\left( 0,1\right)\) independientes.

  2. Para \(i = 1, 2, \ldots, d\) hacer \(T_i = 2V_i - 1\).

  3. Si \(T_1^2 + T_2^2 + \cdots + T_d^2 > 1\) entonces volver al paso 1.

  4. Devolver \(\mathbf{X} = \left( T_1, T_2, \ldots, T_d \right)^t\).

Ver el Ejercicio 1.1 para el caso de \(d=2\).

Usando las fórmulas del “volumen” de una “esfera” \(d\)-dimensional: \[V_d\left( r\right) =\left\{ \begin{array}{ll} \dfrac{\pi^{d/2}r^{d}}{\left( d/2\right) !} & \text{si } d \text{ es par}\\ \dfrac{2^{\left\lfloor \frac{d}{2}\right\rfloor +1}\pi^{\left\lfloor \frac{d}{2}\right\rfloor }r^{d}}{1\cdot3\cdot5\cdots d} & \text{si } d \text{ es impar} \end{array}\right.\] puede verse que el número medio de iteraciones del algoritmo, dado por la constante \(c=\frac{2^{d}}{V_d\left(1 \right)}\), puede llegar a ser enormemente grande. Así, si \(d=2\) se tiene \(c=1.27\), si \(d=3\) se tiene \(c=1.91\), si \(d=4\) entonces \(c=3.24\) y para \(d=10\) resulta \(c=401.5\) que es un valor que hace que el algoritmo sea tremendamente lento en dimensión \(10\). Esto está relacionado con la maldición de la dimensionalidad (curse of dimensionality), a medida que aumenta el número de dimensiones el volumen de la “frontera” crece exponencialmente (ver p.e. Fernández-Casal et al, 2021, Sección 1.4).