Generar puntos de manera eficiente entre el círculo unitario y el cuadrado unitario

Question

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 sólo proporciona una $1-\pi/4$ (~21,4%) de eficiencia.

¿Hay alguna manera de que pueda tomar muestras de manera más eficiente?

Answer 1

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.

Se trata de los siguientes pasos sencillos:

1) Dibuja dos muestras uniformes estándar: $$ u_1 \sim Unif(0,1)\\ u_2 \sim Unif(0,1). $$

2a) Aplica la siguiente transformación de cizalladura al punto $\min\{u_1,u_2\}, \max\{u_1,u_2\}$ (los puntos del triángulo inferior derecho se reflejan en el triángulo superior izquierdo y serán "no reflejados" en 2b): $$ \begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} 1\\1 \end{bmatrix} + \begin{bmatrix} \frac{\sqrt{2}}{2} & -1\\ \frac{\sqrt{2}}{2} - 1 & 0\\ \end{bmatrix} \, \begin{bmatrix} \min\{u_1,u_2\}\\ \max\{u_1,u_2\}\\ \end{bmatrix}. $$

2b) Intercambio $x$ y $y$ si $u_1 > u_2$ .

3) Rechazar la muestra si está dentro del círculo de la unidad (la aceptación debe estar en torno al 72%), es decir: $$ x^2 + y^2 < 1. $$

La intuición detrás de este algoritmo se muestra en la figura.

Los pasos 2a y 2b pueden fusionarse en un solo paso:

2) Aplicar la transformación de cizalladura e intercambiar $$ x = 1 + \frac{\sqrt{2}}{2} \min(u_1, u_2) - u_2\\ y = 1 + \frac{\sqrt{2}}{2} \min(u_1, u_2) - u_1 $$

El siguiente código implementa el algoritmo anterior (y lo prueba utilizando el código de @whuber).

n.sim <- 1e6
x.time <- system.time({
    # Draw two standard uniform samples
    u_1 <- runif(n.sim)
    u_2 <- runif(n.sim)
    # Apply shear transformation and swap
    tmp <- 1 + sqrt(2)/2 * pmin(u_1, u_2)
    x <- tmp - u_2
    y <- tmp - u_1
    # Reject if inside circle
    accept <- x^2 + y^2 > 1
    x <- x[accept]
    y <- y[accept]
    n <- length(x)
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

Algunas pruebas rápidas arrojan los siguientes resultados.

Algoritmo https://stats.stackexchange.com/a/258349 . Al mejor de 3: 0,33 segundos por millón de puntos.

Este algoritmo. El mejor de 3: 0,18 segundos por millón de puntos.

Answer 2

¿Servirán dos millones de puntos por segundo?

La distribución es simétrica: sólo tenemos que calcular la distribución para un octavo del círculo completo y luego copiarla alrededor de los otros octantes. En coordenadas polares $(r,\theta)$ la distribución acumulativa del ángulo $\Theta$ para la ubicación aleatoria $(X,Y)$ en el valor $\theta$ viene dada por el área entre el triángulo $(0,0), (1,0), (1,\tan\theta)$ y el arco del círculo que se extiende desde $(1,0)$ a $(\cos\theta,\sin\theta)$ . Por lo tanto, es proporcional a

$$F_\Theta(\theta) = \Pr(\Theta \le \theta) \propto \frac{1}{2}\tan(\theta) - \frac{\theta}{2},$$

por lo que su densidad es

$$f_\Theta(\theta) = \frac{d}{d\theta} F_\Theta(\theta) \propto \tan^2(\theta).$$

Podemos muestrear a partir de esta densidad utilizando, por ejemplo, un método de rechazo (que tiene una eficiencia $8/\pi-2 \approx 54.6479\%$ ).

La densidad condicional de la coordenada radial $R$ es proporcional a $rdr$ entre $r=1$ y $r=\sec\theta$ . Eso se puede muestrear con una fácil inversión de la FCD.

Si generamos muestras independientes $(r_i,\theta_i)$ , conversión de nuevo a coordenadas cartesianas $(x_i,y_i)$ muestra este octante. Dado que las muestras son independientes, el intercambio aleatorio de las coordenadas produce una muestra aleatoria independiente del primer cuadrante, tal como se desea. (Los intercambios aleatorios requieren la generación de una sola variable binomial para determinar cuántas de las realizaciones hay que intercambiar).

Cada una de estas realizaciones de $(X,Y)$ requiere, por término medio, una variante uniforme (para $R$ ) y además $1/(8\pi-2)$ por dos variantes uniformes (para $\Theta$ ) y una pequeña cantidad de cálculos (rápidos). Es decir $4/(\pi-4) \approx 4.66$ variantes por punto (que, por supuesto, tiene dos coordenadas). Los detalles completos se encuentran en el ejemplo de código más abajo. Esta figura muestra 10.000 de los más de medio millón de puntos generados.

Aquí está el R código que produjo esta simulación y la cronometró.

n.sim <- 1e6
x.time <- system.time({
  # Generate trial angles `theta`
  theta <- sqrt(runif(n.sim)) * pi/4
  # Rejection step.
  theta <- theta[runif(n.sim) * 4 * theta <= pi * tan(theta)^2]
  # Generate radial coordinates `r`.
  n <- length(theta)
  r <- sqrt(1 + runif(n) * tan(theta)^2)
  # Convert to Cartesian coordinates.
  # (The products will generate a full circle)
  x <- r * cos(theta) #* c(1,1,-1,-1)
  y <- r * sin(theta) #* c(1,-1,1,-1)
  # Swap approximately half the coordinates.
  k <- rbinom(1, n, 1/2)
  if (k > 0) {
    z <- y[1:k]
    y[1:k] <- x[1:k]
    x[1:k] <- z
  }
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

Answer 3

Bueno, más eficientemente se puede hacer, pero espero que no estés buscando más rápido .

La idea sería tomar una muestra de un $x$ con una densidad proporcional a la longitud del corte azul vertical sobre cada $x$ valor:

$$ f(x) = 1-\sqrt{1-x^2}. $$

Wolfram le ayuda a integrar esa :

$$ \int_0^x f(y)dy = -\frac{1}{2}x\sqrt{1-x^2}+x-\frac{1}{2}\arcsin x.$$

Así que la función de distribución acumulativa $F$ sería esta expresión, escalada para integrarse en 1 (es decir, dividida por $ \int_0^1 f(y)dy$ ).

Ahora, para generar su $x$ valor, escoge un número aleatorio $t$ distribuido uniformemente entre $0$ y $1$ . Entonces encuentra $x$ tal que $F(x)=t$ . Es decir, tenemos que invertir la FCD ( muestreo por transformación inversa ). Esto se puede hacer, pero no es fácil. Ni rápido.

Por último, teniendo en cuenta $x$ , elija una al azar $y$ que se distribuye uniformemente entre $\sqrt{1-x^2}$ y $1$ .

A continuación se muestra el código R. Tenga en cuenta que estoy pre-evaluando el CDF en una cuadrícula de $x$ valores, e incluso entonces esto lleva bastantes minutos.

Probablemente puedas acelerar la inversión de la FCD bastante si inviertes un poco de pensamiento. Por otra parte, pensar duele. Yo personalmente optaría por el muestreo de rechazo, que es más rápido y mucho menos propenso a errores, a menos que tuviera muy buenas razones para no hacerlo.

epsilon <- 1e-6
xx <- seq(0,1,by=epsilon)
x.cdf <- function(x) x-(x*sqrt(1-x^2)+asin(x))/2
xx.cdf <- x.cdf(xx)/x.cdf(1)

nn <- 1e4
rr <- matrix(nrow=nn,ncol=2)
set.seed(1)
pb <- winProgressBar(max=nn)
for ( ii in 1:nn ) {
    setWinProgressBar(pb,ii,paste(ii,"of",nn))
    x <- max(xx[xx.cdf<runif(1)])
    y <- runif(1,sqrt(1-x^2),1)
    rr[ii,] <- c(x,y)
}
close(pb)

plot(rr,pch=19,cex=.3,xlab="",ylab="")