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

17

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

entrez la description de l'image ici

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

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

Cam.Davidson.Pilon
la source
6
Astuce : utilisez la symétrie pour doubler trivialement votre efficacité.
Cardinal
3
Oh comme: si la valeur est (0,0), cela peut être mappé à (1,1)? J'adore cette idée
Cam.Davidson.Pilon
@cardinal Ne devrait-il pas 4x multiplier l'efficacité? Vous pouvez échantillonner dans , puis le mettre en miroir sur l'axe x, l'axe y et l'origine. [0,,1]×[0,,1]
Martin Krämer
1
@Martin: Dans les quatre régions symétriques, vous avez des chevauchements, que vous devez traiter plus attentivement.
Cardinal
3
@ Martin: Si je comprends bien ce que vous décrivez, qui n'augmente pas l'efficacité du tout . (Vous avez trouvé un point, et maintenant vous en connaissez trois autres --- dans une zone quatre fois plus grande --- qui se trouvent ou non dans le disque unitaire avec une probabilité selon que fait. Comment cela aide-t-il?) Le point d'augmenter l'efficacité est d'augmenter la probabilité d'acceptation pour chaque ( x , y ) généré. Peut-être que je suis celui qui est dense? (x,y)(x,y)
Cardinal

Réponses:

10

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,θ)Θ à 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,Oui)θ(0,0),(1,0),(1,tanθ) à ( cos θ , sin θ ) . Il est ainsi proportionnel à(1,0)(cosθ,sinθ)

FΘ(θ)=Pr(Θθ)12bronzer(θ)-θ2,

d'où sa densité

FΘ(θ)=θFΘ(θ)bronzer2(θ).

Nous pouvons échantillonner à partir de cette densité en utilisant, par exemple, une méthode de rejet (qui a une efficacité de ).8/π-254,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.Rrrr=1r=secondeθ

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.)(rje,θje)(Xje,yje)

Chacune de ces réalisations de nécessite, en moyenne, une variable uniforme (pour R ) plus 1 / ( 8 π(X,Oui)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π-2)Θ4/(π4)4.66

Figure

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")
whuber
la source
1
Je ne comprends pas cette phrase: "Parce que les échantillons sont indépendants, l'échange systématique des coordonnées chaque deuxième échantillon produit un échantillon aléatoire indépendant du premier quadrant, comme souhaité." Il me semble que l'échange systématique des coordonnées tous les deux échantillons produit des échantillons très dépendants. Par exemple, il me semble que votre implémentation en code génère un demi-million d'échantillons d'affilée à partir du même octant?
A. Rex du
7
À strictement parler, cette approche ne fonctionne pas tout à fait (pour les points iid) car elle génère un nombre identique d'échantillons dans les deux octants: Les points d'échantillonnage sont donc dépendants. Maintenant, si vous lancez des pièces non biaisées pour déterminer l'octant de chaque échantillon ...
cardinal
1
@Cardinal vous avez raison; Je vais résoudre ce problème - sans (asymptotiquement) augmenter le nombre de variables aléatoires à générer!
whuber
2
n2n
2
2sin(θ)2(4π)/(π2)75%
13

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:

u1UnjeF(0,1)u2UnjeF(0,1).

min{u1,u2},max{u1,u2}

[Xy]=[11]+[22-122-10][min{u1,u2}max{u1,u2}].

Xyu1>u2

X2+y2<1.

L'intuition derrière cet algorithme est montrée dans la figure. entrez la description de l'image ici

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

X=1+22min(u1,u2)-u2y=1+22min(u1,u2)-u1

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.

Luca Citi
la source
3
+1 Très bien fait! Merci de partager une solution réfléchie, intelligente et simple.
whuber
Bonne idée! Je pensais à un mappage de l'unité sq à cette partie, mais je n'ai pas pensé à un mappage imparfait , puis à un schéma de rejet. Merci d'avoir élargi mon esprit!
Cam.Davidson.Pilon
7

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

XX

F(X)=1-1-X2.

Wolfram vous aide à intégrer cela :

0XF(y)y=-12X1-X2+X-12arcsinX.

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

Maintenant, pour générer votre X valeur, choisissez un nombre aléatoire t, uniformément répartie entre 0 et 1. Trouvez ensuiteX 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 1-X2 et 1.

Ci-dessous est le code R. Notez que je pré-évalue le CDF sur une grille deX 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="")

randoms

S. Kolassa - Réintégrer Monica
la source
Je me demande si l'utilisation des polynômes de Chebyshev pour approximer le CDF améliorerait la vitesse d'évaluation.
Sycorax dit Réintégrer Monica le
@Sycorax, non sans modifications; voir par exemple le traitement chebfun des singularités algébriques aux points limites .
JM n'est pas statisticien le