Jak wygenerować równomiernie rozmieszczone punkty na powierzchni sfery trójwymiarowej?

68

Zastanawiam się, jak wygenerować równomiernie rozmieszczone punkty na powierzchni sfery jednostki 3-d? Również po wygenerowaniu tych punktów, jaki jest najlepszy sposób na wizualizację i sprawdzenie, czy są one naprawdę jednolite na powierzchni x2+y2+z2=1 ?

Qiang Li
źródło
Jeśli przez mundur masz na myśli „regularny”, nie ma innego sposobu na zrobienie tego poza n = 2, 4, 6, 8, 12, 20.
Marcos
1
co jest nie tak z próbką z MultiVariateGaussian i tym wektorze po prostu go znormalizować: v = MultivariateNormal(torch.zeros(10000), torch.eye(10000))a następnie v = v/v.norm(10000)
Pinokio

Odpowiedzi:

72

Standardową metodą jest wygenerowanie trzech standardowych norm i zbudowanie z nich wektora jednostkowego. To znaczy, gdy i λ 2 = X 2 1 + X 2 2 + X 2 3 , wówczas ( X 1 / λ , X 2 / λ , X 3 / λ ) jest równomiernie rozłożone na kula. Ta metoda działa również dobrze w przypadku sfer wymiarowych d .XiN(0,1)λ2=X12+X22+X32(X1/λ,X2/λ,X3/λ)d

W 3D możesz użyć próbkowania odrzucenia: narysuj z równomiernego rozkładu [ - 1 , 1 ], aż długość ( X 1 , X 2 , X 3 ) będzie mniejsza lub równa 1, a następnie - podobnie jak w przypadku poprzednia metoda - znormalizuj wektor do długości jednostkowej. Oczekiwana liczba prób na punkt kulisty wynosi 2 3 / ( 4 π / 3 ) = 1,91. W wyższych wymiarach oczekiwana liczba prób staje się tak duża, że ​​szybko staje się to niewykonalne.Xi[1,1](X1,X2,X3)23/(4π/3)

Istnieje wiele sposobów sprawdzenia jednolitości . Zgrabnym, choć nieco wymagającym obliczeniowo, jest funkcja K Ripleya . Oczekiwana liczba punktów w (3D euklidesowej) odległości dowolnego miejsca na kuli jest proporcjonalna do powierzchni kuli w odległości ρ , która jest równa π ρ 2 . Obliczając wszystkie odległości między punktami, możesz porównać dane z tym ideałem.ρρπρ2

Ogólne zasady konstruowania grafiki statystycznej sugerują, że dobrym sposobem na porównanie jest wykreślenie reszt stabilizowanych wariancją względem i = 1 , 2 , , n ( n - 1 ) / 2 = m gdzie d [ i ] jest i th najmniejsza wzajemnych odległościach i e i = 2 ei(d[i]ei)i=1,2,,n(n1)/2=md[i]ith . Fabuła powinna być bliska zeru. (To podejście jest niekonwencjonalne).ei=2i/m

Oto obraz 100 niezależnych losowań z jednolitego rozkładu sferycznego uzyskanych pierwszą metodą:

100 jednolitych punktów kulistych

Oto wykres diagnostyczny odległości:

Wykres diagnostyczny

Skala y sugeruje, że wszystkie te wartości są bliskie zeru.

Oto suma 100 takich wykresów, aby zasugerować, jakie odchylenia wielkości mogą być znaczącymi wskaźnikami niejednorodności:

Symulowane wartości

(Te fabuły wyglądają okropnie jak mosty Browna ... mogą tu się kryć ciekawe odkrycia teoretyczne.)

Na koniec, oto wykres diagnostyczny dla zestawu 100 jednolitych losowych punktów plus kolejne 41 punktów równomiernie rozmieszczonych tylko na górnej półkuli:

Symulowane wartości nierównomierne

W stosunku do rozkładu równomiernego pokazuje znaczny spadek średnich odległości między punktami w zakresie jednej półkuli. To samo w sobie nie ma znaczenia, ale użyteczną informacją jest to, że coś jest nierównomierne w skali jednej półkuli. W efekcie ten wykres łatwo wykrywa, że ​​jedna półkula ma inną gęstość niż druga. (Prostszy test chi-kwadrat zrobiłby to z większą mocą , gdybyś z góry wiedział, którą półkulę przetestować spośród nieskończenie wielu możliwych).

Whuber
źródło
(X1/λ,X2/λ,X3/λ)
23
XN(0,In)Inn×nQQXN(0,In)XY=X/X2YQ=QX/QX2=QX/X2Q. Ponieważ jest niezmienny dla rotacji, podobnie jest , a ponieważ prawie na pewno, to musi być równomiernie rozmieszczony na kuli. XYY2=1
kardynał
3
@Mike Nie, ponieważ równomierny rozkład szerokości geograficznej nie daje równomiernego rozkładu na kuli. (Większość powierzchni kuli leży na niższych szerokościach geograficznych w pobliżu równika, z dala od biegunów. Zamiast tego potrzebujesz równomiernego rozkładu .)ϕcos(ϕ)
whuber
1
@Ahsan Ponieważ macierze ortogonalne tworzą przechodnią grupę transformacji zachowujących powierzchnię sfery, rozkład jest jednolity na podzbiorze sfery formy : ale to jest cała sfera. X/||X||2
whuber
1
@Cesar „Jednolity rozkład” (na kuli).
whuber
19

Oto trochę raczej prosty kod R.

n     <- 100000                  # large enough for meaningful tests
z     <- 2*runif(n) - 1          # uniform on [-1, 1]
theta <- 2*pi*runif(n) - pi      # uniform on [-pi, pi]
x     <- sin(theta)*sqrt(1-z^2)  # based on angle
y     <- cos(theta)*sqrt(1-z^2)     

Z konstrukcji bardzo łatwo jest zauważyć, że a więc ale jeśli trzeba to przetestować, tox2+y2=1z2x2+y2+z2=1

mean(x^2+y^2+z^2)  # should be 1
var(x^2+y^2+z^2)   # should be 0

i łatwe do sprawdzenia, że każdy z i są równomiernie rozmieszczone na ( pewnością jest), zxy[1,1]z

plot.ecdf(x)  # should be uniform on [-1, 1]
plot.ecdf(y)
plot.ecdf(z)

Oczywiście, biorąc pod uwagę wartość , i są równomiernie rozmieszczone wokół okręgu o promieniu i można to sprawdzić, patrząc na rozkład arcus tangensu ich stosunku. Ale ponieważ ma taki sam rozkład krańcowy jak i jak , podobne stwierdzenie jest prawdziwe dla każdej pary i to również można przetestować. zxy1z2zxy

plot.ecdf(atan2(x,y)) # should be uniform on [-pi, pi]
plot.ecdf(atan2(y,z))
plot.ecdf(atan2(z,x))

Jeśli nadal nie będzie przekonany, następnym krokiem byłoby przyjrzenie się dowolnemu obrotowi trójwymiarowemu lub liczbie punktów mieszczących się w danym kącie bryłowym, ale zaczyna się to komplikować i myślę, że nie jest to konieczne.

Henz
źródło
Zastanawiam się tylko, czy twoja metoda generowania punktów (x, y, z) jest zasadniczo taka sama jak metoda Whubera?
Qiang Li
3
Nie, nie jest: whuber używa trzech liczb losowych, a ja dwóch. Mój jest szczególnym przypadkiem „wygenerowania punktu na o odpowiedniej gęstości [proporcjonalnej do ], a następnie obniżenia wymiaru”. Tutaj dogodnie ponieważ formalnie jest to 2-kula . ( 1 - z 2 ) n / 2 - 1 n = 2[1,1](1z2)n/21n=2
Henry
3
Lub, bardziej ogólnie, wygeneruj jednolite punkty na mapie za pomocą dowolnego rzutu równego obszaru (twój jest cylindryczny rzutu równego obszaru), a następnie rzutuj z powrotem. (+1)
whuber
@whuber: Rzeczywiście. Offtopic, ale dla wszystkich zainteresowanych mam interaktywny wybór światowych odwzorowań tutaj , z których niektóre są równe obszar
Henry
2
Jest to prawie standardowe podejście stosowane w grafice komputerowej, oparte na Twierdzeniu Hat-Boxa Archimedesa: mathworld.wolfram.com/ArchimedesHat-BoxTheorem.html
Edward KMETT
10

Jeśli chcesz próbkować punkty równomiernie rozmieszczone na kuli 3D (tj. Powierzchni kuli 3D), użyj prostego odrzucenia lub metody Marsaglii (Ann. Math. Statist., 43 (1972), s. 645– 646). W przypadku niskich wymiarów współczynnik odrzucenia jest dość niski.

Jeśli chcesz generować losowe punkty ze sfer i kulek o wyższych wymiarach, zależy to od celu i skali symulacji. Jeśli nie chcesz wykonywać dużych symulacji, skorzystaj z metody Mullera (Commun. ACM, 2 (1959), s. 19–20) lub jej wersji „kulowej” (patrz artykuł Harman & Lacko cytowany powyżej). To jest:

aby uzyskać próbkę równomiernie rozłożoną na kuli n (powierzchni) 1) wygeneruj X z n-wymiarowego standardowego rozkładu normalnego 2) podziel każdy składnik X przez euklidesową normę X

aby uzyskać próbkę równomiernie rozłożoną na kuli n (wewnętrznej) 1) wygeneruj X z (n + 2) wymiarowego standardowego rozkładu normalnego 2) podziel każdy składnik X przez euklidesową normę X i weź tylko pierwsze n składników

Jeśli chcesz wykonywać duże symulacje, powinieneś zbadać bardziej wyspecjalizowane metody. Na życzenie mogę wysłać artykuł Harmana i Lacko na temat metod dystrybucji warunkowej, który zapewnia klasyfikację i uogólnienie niektórych algorytmów wspomnianych w tej dyskusji. Kontakt jest dostępny na mojej stronie internetowej (http://www.iam.fmph.uniba.sk/ospm/Lacko)

Jeśli chcesz sprawdzić, czy twoje punkty są naprawdę jednolite na powierzchni lub we wnętrzu kuli, spójrz na marginesy (wszystkie powinny być takie same, ze względu na niezmienność rotacyjną, kwadratowa norma rzutowanej próbki jest rozkładana w wersji beta).

V Lacko
źródło
co jest nie tak z próbką z MultiVariateGaussian i tym wektorze po prostu go znormalizować: v = MultivariateNormal(torch.zeros(10000), torch.eye(10000))a następnie v = v/v.norm(10000)
Pinokio
8

Miałem podobny problem (sferę n) podczas mojego doktoratu i jeden z lokalnych „ekspertów” zasugerował odrzucenie próbkowania z n-sześcianu! To oczywiście zajęłoby wiek wszechświata, gdy patrzyłem na n w kolejności hunderdów.

Algorytm, którego użyłem, jest bardzo prosty i opublikowany w:

WP Petersen i A. Bernasconic Jednolite pobieranie próbek z kuli n: Metoda izotropowa Raport techniczny, TR-97-06, Szwajcarskie Centrum Informatyki Naukowej

Mam również ten artykuł w swojej bibliografii, na który nie spojrzałem. Może ci się przydać.

Harman, R. & Lacko, V. Na decompositional algorytmów jednolitego próbkowania od -spheres i -balls Journal of analizie wieloczynnikowej, 2010nn

emakaliczny
źródło
czy można zamieścić linki, w których można znaleźć pełny tekst tych odniesień? dzięki.
Qiang Li
Nie mam na sobie papieru, ale ta strona wydaje się opisywać algorytm (i kilka innych) mlahanas.de/Math/nsphere.htm
emakalic
3
Jak rozumiem (z artykułu Petersena i Bernasconica) dla kuli dwuwymiarowej można wygenerować promień, podnosząc U (0,1) zmieniając na moc (1 / d), a ostatni kąt jako U (0,2 ) zmienia się. Kąty pośrednie można uzyskać jako , gdzie to . Dla mnie to brzmi dość prosto. Zastanawiam się nad tym: jeśli zastosuję quasi-losową sekwencję dla moich mundurów, czy dostanę także uprzejmość w piłce? πC.asin(uk)C1
πΓ(k2+0.5)Γ(k2+1)
Mohit
3

Miałem już ten problem, a oto alternatywa, którą znalazłem,

Jeśli chodzi o sam rozkład, stwierdziłem, że wzór, który działa przyzwoicie, polega na użyciu współrzędnych biegunowych (faktycznie używam opracowanej odmiany współrzędnych biegunowych), a następnie przekształceniu na współrzędne kartezjańskie.

Promień jest oczywiście promieniem kuli, na której drukujesz. Następnie masz drugą wartość kąta na płaszczyźnie płaskiej, a następnie trzecią wartość, która jest kątem powyżej lub poniżej tej płaszczyzny.

Aby uzyskać przyzwoity rozkład, załóżmy, że U jest równomiernie rozmieszczoną liczbą losową, r jest promieniem, a jest drugą współrzędną biegunową, a b jest trzecią współrzędną biegunową,

a = U * 360 b = U + U-1, a następnie konwersja na kartezjański za pomocą x = r * sin (b) sin (a) z = r sin (b) cos (a) y = r sin (b)

Niedawno znalazłem następujący, który jest lepszy matematycznie: a = 2 (pi) * U b = cos ^ -1 (2U-1)

Właściwie niewiele różni się od mojej pierwotnej formuły, chociaż moje są stopniami vs radianami.

Ta najnowsza wersja podobno może być używana do hiperfer, choć nie wspomniano, jak to osiągnąć.

Chociaż wizualnie sprawdzam jednolitość za pomocą dość taniej metody tworzenia map dla Homeworld 2, a następnie „odtwarzania” tych map. W rzeczywistości, ponieważ mapy są tworzone za pomocą skryptów lua, możesz wbudować formułę bezpośrednio w mapę, a tym samym sprawdzić wiele próbek bez wychodzenia z gry. Być może nie naukowy, ale to dobra metoda wizualnego oglądania wyników.

TheDerpyAlicorn
źródło
2

Oto pseudokod:

  1. vMultiVariateGaussian(μ,σI)
  2. v=vv

W pytorch:

v = MultivariateNormal(torch.zeros(10000), torch.eye(10000))
v = v/v.norm(2)

Nie rozumiem tego wystarczająco dobrze, ale whuber powiedział mi, że:

v = torch.normal(torch.zeros(10000), torch.eye(10000))
v = v/v.norm(2)

jest również poprawny, tj. próbkowanie z jednowymiarowej normalnej dla każdej współrzędnej.

Pinokio
źródło
0

Moim najlepszym pomysłem byłoby najpierw wygenerowanie zestawu równomiernie rozmieszczonych punktów w dwuwymiarowej przestrzeni, a następnie rzutowanie tych punktów na powierzchnię kuli za pomocą pewnego rodzaju projekcji.

Prawdopodobnie będziesz musiał mieszać i dopasowywać sposób generowania punktów do sposobu ich mapowania. Jeśli chodzi o generowanie punktów 2D, myślę, że zaszyfrowane sekwencje o niskiej rozbieżności byłyby dobrym miejscem do rozpoczęcia (tj. Zaszyfrowana sekwencja Sobola), ponieważ zwykle wytwarza punkty, które nie są „zlepione”. Nie jestem pewien, jakiego rodzaju mapowania użyć, ale Woflram wyskoczył z projekcji Gnonomic ... więc może to może zadziałać?

MATLAB ma przyzwoitą implementację sekwencji o niskiej rozbieżności, które można generować za pomocą q = sobolset(2)i szyfrować za pomocą q = scramble(q). W MATLAB-ie znajduje się także zestaw narzędzi do mapowania z szeregiem różnych funkcji projekcji, których możesz użyć, gdybyś nie chciał sam kodować mapowania i grafiki.

Berk U.
źródło
1
czy którakolwiek z tych prognoz nadal może zachować jednolitość losowości? Ponownie, jak mogę sprawdzić, czy ostateczny rozkład tych punktów jest rzeczywiście równomiernie rozmieszczony na powierzchni kuli? Dzięki.
Qiang Li,
Przepraszam, mówiłem tylko hipotetycznie ... Myślę, że funkcje mapowania w MATLABie pozwoliłyby ci to sprawdzić, ponieważ mają w sobie pewne wizualizacje. Jeśli nie, znalazłem też fajną stronę internetową, która mówi o tym, jak generować równomiernie rozmieszczone punkty na kuli w 3D za pomocą takich rzeczy, jak losowe kąty itp. Mają też trochę kodu C. Spójrz
Berk U.
3
Jednolite losowe punkty na rzucie gnomonicznym nie będą jednorodne na kuli, ponieważ gnomon nie jest równym obszarem. Rzut zaproponowany przez Henry'ego -> (od długości i szerokości geograficznej do prostokąta w ), ma równą powierzchnię. ( λ , sin ( φ ) ) R 2(λ,ϕ)(λ,sin(ϕ))R2
whuber