Symuluj równomierny rozkład na dysku

24

Próbowałem symulować wstrzykiwanie losowych punktów w kółko, tak aby dowolna część koła miała takie samo prawdopodobieństwo wystąpienia wady. Spodziewałem się, że liczenie na pole wynikowego rozkładu będzie zgodne z rozkładem Poissona, jeśli podzielę okrąg na prostokąty o równej powierzchni.

Ponieważ wymaga to jedynie umieszczenia punktów w obszarze kołowym, wstrzyknąłem dwa równomierne rozkłady losowe we współrzędnych biegunowych: R (promień) i (kąt biegunowy).θ

Ale po wykonaniu tego zastrzyku wyraźnie dostaję więcej punktów na środku koła w porównaniu do krawędzi.

wprowadź opis zdjęcia tutaj

Jaki byłby prawidłowy sposób wykonania tego wstrzyknięcia przez okrąg, tak aby punkty były losowo rozmieszczone w kręgu?

Jonjilla
źródło
To pytanie ma dokładny analog na forum geometrii: math.stackexchange.com/questions/87230/…
Aksakal

Odpowiedzi:

35

Chcesz, aby proporcja punktów była jednakowo proporcjonalna do powierzchni, a nie do odległości od początku. Ponieważ powierzchnia jest proporcjonalna do kwadratu odległości, wygeneruj jednolite losowe promienie i weź ich pierwiastki kwadratowe. Połącz to z jednolitym kątem biegunowym.

Jest to szybkie i proste do kodowania, wydajne w wykonaniu (szczególnie na platformie równoległej) i generuje dokładnie określoną liczbę punktów.

Przykład

To działa Rkod ilustrujący algorytm.

n <- 1e4
rho <- sqrt(runif(n))
theta <- runif(n, 0, 2*pi)
x <- rho * cos(theta)
y <- rho * sin(theta)
plot(x, y, pch=19, cex=0.6, col="#00000020")

wprowadź opis zdjęcia tutaj

Whuber
źródło
3

Można zastosować próbkowanie odrzucenia . Oznacza to, że możemy próbkować z jednorodnego rozkładu 2D i wybierać próbki spełniające warunki dysku.

Oto przykład.

x=runif(1e4,-1,1)
y=runif(1e4,-1,1)

d=data.frame(x=x,y=y)
disc_sample=d[d$x^2+d$y^2<1,]
plot(disc_sample)

wprowadź opis zdjęcia tutaj

Haitao Du
źródło
3
Jest to dobra alternatywa dla podejścia przyjętego przez PO. Prosty i wydajny. Jednak tak naprawdę nie odnosi się do pytania, które dotyczy sposobu modyfikacji metody współrzędnych biegunowych w celu uzyskania równomiernie rozłożonych zmiennych. Dlaczego nas to obchodzi? Ze względu na implikacje: gdy już wiesz, jak generować równomiernie rozmieszczone punkty we współrzędnych biegunowych, możesz użyć próbkowania odrzucenia (i innych znanych metod) we współrzędnych biegunowych, aby pobrać próbki z regionów, które mogą być nadmiernie skomplikowane w próbkach we współrzędnych kartezjańskich (pomyśl o hipoklikoidach , na przykład).
whuber
1
π/4
@ Whuber dziękuję za edukowanie mnie poprzez komentowanie mojej odpowiedzi!
Haitao Du
3

Dam ci ogólną n-wymiarową odpowiedź, która oczywiście działa również w przypadku dwuwymiarowej. W trzech wymiarach analog dysku jest objętością litej kuli (kuli).

Omówię dwa podejścia. Jeden z nich nazwałbym „precyzyjnym” , a dostaniesz z nim kompletne rozwiązanie w R. Drugi, który nazywam heurystyczny , i to tylko pomysł, nie ma pełnego rozwiązania.

„Precyzyjne” rozwiązanie

Moje rozwiązanie oparte jest na pracach Marsaglii i Mullera . Zasadniczo dzieje się tak, że wektor Gaussa znormalizowany do swojej normy dałby równomiernie rozmieszczone punkty na hipersferze d-wymiarowej:

wprowadź opis zdjęcia tutaj

re1/re

n <- 1e4
rho <- sqrt(runif(n))
# d - # of dimensions of hyperdisk
d = 2
r = matrix(rnorm(n*d),nrow=n,ncol=d)
x = r/rep(sqrt(rowSums(r^2))/rho,1)
plot(x[,1], x[,2], pch=19, cex=0.6, col="#00000020")

wprowadź opis zdjęcia tutaj

Oto fragment kodu dla przypadku 3D, czyli litej kuli:

library(scatterplot3d)
n <- 1e3
# d - # of dimensions of hyperdisk

d=3
rho <- (runif(n))^(1/d)
r = matrix(rnorm(n*d),nrow=n,ncol=d)
x = r/rep(sqrt(rowSums(r^2))/rho,1)

scatterplot3d(x[,1], x[,2], x[,3])

wprowadź opis zdjęcia tutaj

Podejście heurystyczne

V.n(R)=πn2)Γ(n2)+1)Rn
Rn

Dlaczego ma to znaczenie dla naszego problemu? Załóżmy, że chcesz wygenerować d losowe liczby jednolite, byłyby to losowe punkty wewnątrz hipersześcianu d-wymiarowego. Następnie stosuje się próbkowanie odrzucenia, aby wybrać punkty wewnątrz hipersfery (aka n-ball): ja=1rexja2)<R2)

1re+2)

Aksakal
źródło
@Silverfish, masz rację,
poprawiłem
@Silverfish, jest powolny ze względu na stosowanie zmiennych Gaussa, ale może być szybszy niż zwykłe próbkowanie odrzucenia w przypadku wielowymiarowym, co dla wielu nie jest oczywiste, chociaż jest to inny temat
Aksakal
1/re,re
@ whuber, kopiowałem wklejanie, poprawiłem literówkę dotyczącą mocy kostki. Jeśli użyjemy Gaussa, próbka odrzucenia nie będzie lepsza, więc musielibyśmy użyć czegoś w kształcie dzwonu, który jest szybszy niż Gaussa, masz rację
Aksakal
0

Oto alternatywne rozwiązanie w R:

n <- 1e4
## r <- seq(0, 1, by=1/1000)
r <- runif(n)
rho <- sample(r, size=n, replace=T, prob=r)
theta <- runif(n, 0, 2*pi)
x <- rho * cos(theta)
y <- rho * sin(theta)
plot(x, y, pch=19, cex=0.6, col="#00000020")

wprowadź opis zdjęcia tutaj

Q_Li
źródło
4
Czy możesz wyjaśnić tę odpowiedź zwykłym angielskim? Tak naprawdę nie jesteśmy witryną pomocy dla kodu i odpowiedzi tylko na kod powinny być odradzane.
gung - Przywróć Monikę
5
To rozwiązanie jest intrygujące. Niestety nie jest to całkiem poprawne, choć fabuła wygląda dobrze. Powodem jest to, że ogranicza promienie w próbce do dyskretnego zestawu wartości równo rozmieszczonych między nimi0 i 1. To nie to samo, co zamierzona dystrybucja. W szczególności jest wyjątkowy : nie ma nawet pliku PDF! Jeśli nie jesteś przekonany, uruchom go ponownie, r <- seq(0, 1, by=1/10)aby wyraźniej zobaczyć jego dyskretną naturę.
whuber
1
@whuber Dzięki za zwrócenie na to uwagi. To właściwie moja główna idea rozwiązania. Moje podejście polegało na wygenerowaniu wielu jednolitych kół o różnych promieniach, a dla każdego koła liczba punktów jest proporcjonalna do długości jego promienia. Dlatego na jednostkowej długości okręgów o różnych promieniach liczba punktów jest taka sama. Aby uniknąć dyskretnej natury, możemy pobrać próbki rz Uniform (0,1).
Q_Li