Erzeugen Sie effizient Punkte zwischen Einheitskreis und Einheitsquadrat

Ich möchte Samples aus der hier definierten blauen Region generieren:

Die naive Lösung besteht darin, eine Zurückweisungsabtastung im Einheitenquadrat zu verwenden, dies liefert jedoch nur einen Wirkungsgrad von (~ 21,4%).$1-\pi/4$

Gibt es eine Möglichkeit, wie ich effizienter probieren kann?

Reichen zwei Millionen Punkte pro Sekunde aus?

Die Verteilung ist symmetrisch: Wir müssen die Verteilung nur für ein Achtel des vollen Kreises berechnen und dann um die anderen Oktanten kopieren. In Polarkoordinaten ist die kumulative Verteilung des Winkels Θ für den zufälligen Ort ( X , Y ) beim Wert θ durch die Fläche zwischen den Dreiecken ( 0 , 0 ) , ( 1 , 0 ) , ( 1 , tan θ ) und der Kreisbogen von $(r,\theta)$$\Theta$$(X,Y)$$\theta$$(0,0), (1,0), (1,\tan\theta)$ bis ( cos θ , sin θ ) . Es ist dabei proportional zu$(1,0)$$(\cos\theta,\sin\theta)$

woher seine Dichte ist

Wir können aus dieser Dichte beispielsweise eine Rückweisungsmethode (mit einem Wirkungsgrad von ) entnehmen .$8/\pi-2 \approx 54.6479\%$

Die bedingte Dichte der Radialkoordinate ist proportional zu r d r zwischen r = 1 und r = sec θ . Das kann mit einer einfachen Inversion der CDF abgetastet werden.$R$$rdr$$r=1$$r=\sec\theta$

Wenn wir unabhängige Abtastwerte erzeugen , tastet die Umwandlung in kartesische Koordinaten ( x i , y i ) diesen Oktanten ab. Da die Abtastwerte unabhängig sind, wird durch zufälliges Vertauschen der Koordinaten je nach Wunsch ein unabhängiger Zufallsabtastwert aus dem ersten Quadranten erzeugt. (Für die zufälligen Auslagerungen muss nur eine einzige Binomialvariable generiert werden, um zu bestimmen, wie viele Realisierungen ausgetauscht werden müssen.)$(r_i,\theta_i)$$(x_i,y_i)$

Jede solche Realisierung von erforderlich ist , im Durchschnitt einer einheitlichen Veränderlichen (für R ) sowie 1 / ( 8 π - 2 ) mal zwei einheitliches variates (für Θ ) und eine geringe Menge (schnell) Berechnung. Das sind 4 / ( π - 4 ) ≈ 4,66 Variationen pro Punkt (der natürlich zwei Koordinaten hat). Vollständige Details finden Sie im folgenden Codebeispiel. Diese Zahl zeigt 10.000 von mehr als einer halben Million generierten Punkten.$(X,Y)$$R$$1/(8\pi-2)$$\Theta$$4/(\pi-4) \approx 4.66$

Hier ist der RCode, der diese Simulation erstellt und zeitgesteuert hat.

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

Ich schlage die folgende Lösung vor, die einfacher, effizienter und / oder rechnerisch billiger sein sollte als andere Lösungen von @cardinal, @whuber und @ stephan-kolassa.

Es umfasst die folgenden einfachen Schritte:

1) Zeichnen Sie zwei einheitliche Standardproben:

2a) Wende die folgende Scherumwandlung auf den Punkt (Punkte im unteren rechten Dreieck werden zum oberen linken Dreieck reflektiert und sie werden "nicht reflektiert" in 2b): [ x y ] = [ 1 1 ] + [ $\min\{u_1,u_2\}, \max\{u_1,u_2\}$

2b) Vertausche und y, wenn u 1 > u 2 .$x$$y$$u_1 > u_2$

3) Die Probe ablehnen, wenn sie sich innerhalb des Einheitskreises befindet (die Akzeptanz sollte bei 72% liegen), dh:

Die Intuition hinter diesem Algorithmus ist in der Abbildung dargestellt.

Die Schritte 2a und 2b können zu einem einzigen Schritt zusammengefasst werden:

2) Scherumwandlung anwenden und x = 1 + √ tauschen

Der folgende Code implementiert den obigen Algorithmus (und testet ihn mit @ whubers Code).

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

Einige Schnelltests ergeben die folgenden Ergebnisse.

Algorithmus /stats//a/258349 . Best of 3: 0,33 Sekunden pro Million Punkte.

Dieser Algorithmus. Best of 3: 0,18 Sekunden pro Million Punkte.

Gut, effizienter geht es nicht, aber ich hoffe, Sie suchen nicht schneller .

$x$$x$

Wolfram hilft Ihnen dabei, das zu integrieren :

$F$$\int_0^1 f(y)dy$

$x$$t$$0$$1$$x$$F(x)=t$

$x$$y$$\sqrt{1-x^2}$$1$

$x$

Sie können die CDF-Inversion wahrscheinlich erheblich beschleunigen, wenn Sie ein wenig nachdenken. Andererseits tut das Denken weh. Ich persönlich würde für die Ablehnung Sampling gehen, die schneller ist und weniger fehleranfällig, es sei denn , ich hatte sehr gute Gründe , nicht zu tun .

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