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 ( 0( r , θ)Θ(X,Y)θ i łuk koła rozciągający się od (( 0 , 0 ) , ( 1 , 0 ) ,(1,tanθ) do ( cos θ , sin θ ) . Jest zatem proporcjonalny do( 1 , 0 )(cosθ ,sinθ )
faΘ(θ)=Pr(Θ≤θ)∝12tan(θ)−θ2,
stąd jego gęstość jest
faΘ( θ ) = dreθfaΘ( θ ) ∝ tan2)( θ ) .
Możemy próbkować z tej gęstości, stosując, powiedzmy, metodę odrzucania (która ma wydajność ).8 / π- 2 ≈ 54,6479 %
Gęstość warunkowa współrzędnej promieniowej jest proporcjonalna do r d r między r = 1 i rRr drr = 1 . Można to próbkować przez łatwą inwersję CDF.r = sekθ
Jeśli wygenerujemy niezależne próbki , konwersja z powrotem do współrzędnych kartezjańskich ( x i( rja, θja) 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).( xja, yja)
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)R1 / ( 8 π- 2 )Θ4 / ( π- 4 ) ≈ 4,66
Oto R
kod, 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")