Efektywnie generuj punkty między kołem jednostki a kwadratem jednostki

Chciałbym wygenerować próbki z zdefiniowanego tutaj niebieskiego regionu:

Naiwnym rozwiązaniem jest stosowanie próbkowania odrzucania w jednostce kwadratowej, ale zapewnia to jedynie wydajność $1-\pi/4$ (~ 21,4%).

Czy jest jakiś sposób na bardziej efektywne próbkowanie?

Czy zrobią to dwa miliony punktów na sekundę?

Rozkład jest symetryczny: musimy tylko opracować rozkład dla jednej ósmej pełnego koła, a następnie skopiować go wokół innych oktantów. We współrzędnych biegunowych skumulowany rozkład kąta Θ dla losowej lokalizacji ( X , Y ) przy wartości θ jest określony przez obszar między trójkątem ( $(r,\theta)$$\Theta$$(X,Y)$$\theta$ i łuk koła rozciągający się od$(0,0), (1,0), (1,\tan\theta)$ do ( cos θ , sin θ ) . Jest zatem proporcjonalny do$(1,0)$$(\cos\theta,\sin\theta)$

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

stąd jego gęstość jest

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

Możemy próbkować z tej gęstości, stosując, powiedzmy, metodę odrzucania (która ma wydajność ).$8/\pi-2 \approx 54.6479\%$

Gęstość warunkowa współrzędnej promieniowej jest proporcjonalna do r d r między r = 1 i $R$$rdr$$r=1$ . Można to próbkować przez łatwą inwersję CDF.$r=\sec\theta$

Jeśli wygenerujemy niezależne próbki , konwersja z powrotem do współrzędnych kartezjańskich ( x $(r_i,\theta_i)$ próbki tego oktant. Ponieważ próbki są niezależne, losowa zamiana współrzędnych tworzy niezależną losową próbkę z pierwszej ćwiartki, zgodnie z potrzebami. (Swapy losowe wymagają wygenerowania tylko jednej zmiennej dwumianowej, aby określić liczbę realizacji do zamiany).$(x_i,y_i)$

Każda taka realizacja wymaga średnio jednej zmiany jednolitej (dla R ) plus 1 / ( 8 π - 2 ) razy dwóch zmian jednolitej (dla Θ ) i niewielkiej ilości (szybkiego) obliczenia. To 4 / ( π - 4 ) ≈ 4,66 zmienna na punkt (który oczywiście ma dwie współrzędne). Pełne szczegóły znajdują się w poniższym przykładzie kodu. Ta liczba przedstawia 10 000 z ponad pół miliona wygenerowanych punktów.$(X,Y)$$R$$1/(8\pi-2)$$\Theta$$4/(\pi-4) \approx 4.66$

Oto Rkod, który wyprodukował tę symulację i dokonał jej pomiaru czasu.

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

Proponuję następujące rozwiązanie, które powinno być prostsze, wydajniejsze i / lub tańsze obliczeniowo niż inne południowe firmy @cardinal, @whuber i @ stephan-kolassa.

Obejmuje następujące proste kroki:

1) Narysuj dwie standardowe jednolite próbki:

$$u_1 \sim Unif(0,1)\\ u_2 \sim Unif(0,1).$$

2a) Zastosuj następującą transformację ścinania do punktu (punkty w prawym dolnym trójkącie są odzwierciedlone w lewym górnym trójkącie i będą „nieodbite” w 2b): [ x y ] = [ 1 1 ] + [ $\min\{u_1,u_2\}, \max\{u_1,u_2\}$

$$\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) Zamień i y, jeśli u 1 > u $x$$y$$u_1 > u_2$ .

3) Odrzuć próbkę, jeśli znajduje się w okręgu jednostkowym (akceptacja powinna wynosić około 72%), tj .:

$$x^2 + y^2 < 1.$$

Intuicja tego algorytmu pokazano na rysunku.

Kroki 2a i 2b można połączyć w jeden krok:

2) Zastosuj transformację ścinającą i zamień

$$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$$

Poniższy kod implementuje powyższy algorytm (i testuje go przy użyciu kodu @ 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")

Niektóre szybkie testy dają następujące wyniki.

Algorytm /stats//a/258349 . Najlepsze 3: 0,33 sekundy na milion punktów.

Ten algorytm Najlepsze 3: 0,18 sekundy na milion punktów.

Dobrze, bardziej efektywnie można to zrobić, ale mam nadzieję, że nie szukasz szybszego .

Pomysł polegałby na próbkowaniu $x$ wartość pierwsza, o gęstości proporcjonalnej do długości pionowego niebieskiego wycinka nad każdym $x$ wartość:

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

Wolfram pomaga zintegrować, że :

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

Więc funkcja skumulowanego rozkładu $F$ byłoby to wyrażenie skalowane w celu integracji do 1 (tzn. podzielone przez $\int_0^1 f(y)dy$).

Teraz wygeneruj swój $x$ wartość, wybierz liczbę losową $t$, równomiernie rozmieszczone między $0$ i $1$. Następnie znajdź$x$ takie, że $F(x)=t$. Oznacza to, że musimy odwrócić CDF ( próbkowanie z transformacją odwrotną ). Można to zrobić, ale nie jest to łatwe. Ani szybko.

Wreszcie podane $x$wybierz losowy $y$ który jest równomiernie rozmieszczony między $\sqrt{1-x^2}$ i $1$.

Poniżej znajduje się kod R. Zauważ, że wstępnie oceniam CDF na poziomie$x$ wartości, a nawet wtedy zajmuje to kilka minut.

Prawdopodobnie możesz znacznie przyspieszyć inwersję CDF, jeśli zainwestujesz trochę myślenia. Z drugiej strony myślenie boli. Osobiście wybrałbym próbkowanie odrzucania, które jest szybsze i znacznie mniej podatne na błędy, chyba że mam bardzo dobre powody, aby tego nie robić .

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