Générez efficacement des points entre le cercle unitaire et le carré unitaire

Je voudrais générer des échantillons à partir de la région bleue définie ici:

La solution naïve consiste à utiliser l'échantillonnage de rejet dans le carré unitaire, mais cela n'offre qu'une efficacité de $1-\pi/4$ (~ 21,4%).

Existe-t-il un moyen de échantillonner plus efficacement?

Est-ce que deux millions de points par seconde suffiront?

La distribution est symétrique: il suffit de calculer la distribution pour un huitième du cercle complet, puis de la copier autour des autres octants. En coordonnées polaires , la distribution cumulée de l'angle Θ pour l'emplacement aléatoire$(r,\theta)$$\Theta$ à la valeur θ est donnée par l'aire entre le triangle ( 0 , 0 ) , ( 1 , 0 ) , ( 1 , tan θ ) et l'arc de cercle s'étendant de$(X,Y)$$\theta$$(0,0), (1,0), (1,\tan\theta)$ à ( cos θ , sin θ ) . Il est ainsi proportionnel à$(1,0)$$(\cos\theta,\sin\theta)$

d'où sa densité

Nous pouvons échantillonner à partir de cette densité en utilisant, par exemple, une méthode de rejet (qui a une efficacité de ).$8/\pi-2 \approx 54.6479\%$

La densité conditionnelle de la coordonnée radiale est proportionnelle à r d r entre r = 1 et r = sec θ . Cela peut être échantillonné avec une inversion facile du CDF.$R$$rdr$$r=1$$r=\sec\theta$

Si nous générons des échantillons indépendants , la conversion en coordonnées cartésiennes ( x i , y i ) échantillonne cet octant. Du fait que les échantillons sont indépendants, un échange aléatoire des coordonnées produit un échantillon aléatoire indépendant du premier quadrant, comme souhaité. (Les échanges aléatoires nécessitent de générer une seule variable binomiale pour déterminer le nombre de réalisations à échanger.)$(r_i,\theta_i)$$(x_i,y_i)$

Chacune de ces réalisations de nécessite, en moyenne, une variable uniforme (pour R ) plus 1 / ( 8 $(X,Y)$$R$ fois deux variables uniformes (pour Θ ) et une petite quantité de calcul (rapide). C'est 4 / ( π - 4 ) ≈ 4,66 varie par point (qui, bien sûr, a deux coordonnées). Tous les détails sont dans l'exemple de code ci-dessous. Ce chiffre représente 10 000 des plus d'un demi-million de points générés.$1/(8\pi-2)$$\Theta$$4/(\pi-4) \approx 4.66$

Voici le Rcode qui a produit cette simulation et l'a chronométré.

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")

Je propose la solution suivante, qui devrait être plus simple, plus efficace et / ou moins coûteuse que les autres soutiens de @cardinal, @whuber et @ stephan-kolassa jusqu'à présent.

Cela implique les étapes simples suivantes:

$\min\{u_1,u_2\}, \max\{u_1,u_2\}$

$x$$y$$u_1 > u_2$

L'intuition derrière cet algorithme est montrée dans la figure.

Les étapes 2a et 2b peuvent être fusionnées en une seule étape:

Le code suivant implémente l'algorithme ci-dessus (et le teste en utilisant le code 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")

Quelques tests rapides donnent les résultats suivants.

Algorithme /stats//a/258349 . Meilleur de 3: 0,33 seconde par million de points.

Cet algorithme. Meilleur de 3: 0,18 seconde par million de points.

Eh bien, plus efficacement peut être fait, mais j'espère que vous ne cherchez pas plus rapidement .

$x$$x$

Wolfram vous aide à intégrer cela :

Ainsi, la fonction de distribution cumulative $F$ serait cette expression, mise à l'échelle pour s'intégrer à 1 (c.-à-d. divisée par $\int_0^1 f(y)dy$).

Maintenant, pour générer votre $x$ valeur, choisissez un nombre aléatoire $t$, uniformément répartie entre $0$ et $1$. Trouvez ensuite$x$ tel que $F(x)=t$. Autrement dit, nous devons inverser le CDF ( échantillonnage à transformée inverse ). Cela peut être fait, mais ce n'est pas facile. Ni rapide.

Enfin, étant donné $x$, choisissez un hasard $y$ qui est uniformément réparti entre $\sqrt{1-x^2}$ et $1$.

Ci-dessous est le code R. Notez que je pré-évalue le CDF sur une grille de$x$ valeurs, et même alors cela prend quelques minutes.

Vous pouvez probablement accélérer l'inversion du CDF un peu si vous réfléchissez. Là encore, la pensée fait mal. Personnellement, je choisirais l'échantillonnage de rejet, qui est plus rapide et beaucoup moins sujet aux erreurs, sauf si j'avais de très bonnes raisons de ne pas le faire.

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="")