Zufällige 3D Vektoren innerhalb einer Kugel

Werner Salomon

SeppJ schrieb:

@Werner: Bist du da sicher, dass die Verteilung deiner Längen stimmt? Dein Code ist mir zu verschlungen, um ihn jetzt noch zu entziffern. Du ziehst nur die dritte Wurzel, oder? Das würde nämlich nicht stimmen, da du dann nämlich entweder die Verteilung bei R1 als 0 hast oder die Verteilung verschoben (und somit die falsche Steigung) hast, je nachdem, was du da genau machst

Ja ich ziehe die 3.Wurzel aus einer Verteilung von x0^3 bis x1^3 (s. Cubic<>) - das ist genau das, was otze vorschlägt

otze schrieb:

Mein Ansatz wäre: wir teilen das Volumen der Kugel in unendlich viele Kugeloberflächen auf. dann wird ein Wert des Radius mit dem Verhältnis gezogen, den seine (infinitesimal dicke) Kugeloberfläche am Integral über alle Oberflächen hat. Das ist praktischerweise das Volumen der Schicht

was dann kommt habe ich zwar nicht verstanden - insbesondere die Ausgangsgleichung, aber das Ergebnis

otze schrieb:

<=> r = (q*(1-a^3)+a3)^(1/3)

ist identisch mit dem, was die Verteilung class Cubic<> macht - siehe Zeile 74 und 80. Man ersetze a in otzes Ergebnis durch das Verhältnis x0/x1 aus Cubic<> und rechne es selber nach.

otze schrieb:

vermutlich habe ich mich irgendwo verrechnet.

nö - bis hierher anscheinend nicht.

otze schrieb:

Und dann 2 Winkel gleichverteilt ziehen, zusammen mit dem Radius müsste das einen gleichverteilten Vektor auf der Kugelschicht mit Radius r erzeugen.

warum das falsch ist, ist am Anfang von http://mathworld.wolfram.com/SpherePointPicking.html erklärt und durch eine Graphik demonstriert.

SeppJ schrieb:

Die Längen kann man dann mit Punkten auf einer Einheitssphäre kombinieren, wie hier schon gezeigt wurde. Herleitung findet man hier:
http://mathworld.wolfram.com/SpherePointPicking.html
Das ist zum Glück ein ganzes Ende einfacher.

Danke für den Link SeppJ,
Dort wird genau das beschrieben, was ich in der Klasse Sphere<> codiert habe, nur das die Bezeichner andere sind und die Intervalle der Winkel phasenverschoben sind - z.B. [0,2PI) statt [-PI,+PI). Aber das sollte das Ergebnis nicht verändern.

[Rewind] schrieb:

@Werner: wie lange hast du für den Code gebraucht?

weiß ich nicht mehr genau. Weniger als eine Stunde auf jeden Fall - vielleicht 45min, die meiste Zeit habe ich dafür gebraucht, aus der Formel für die Oberfläche der Kugelkalotte die Verteilung für den Einheitsvektor zu zaubern. Das ist das, was in den Zeilen 55 bis 57 steht.

Gruß
Werner

Werner Salomon

SeppJ schrieb:

Ich habe jetzt stundenlang an der exakten Verteilung der Längen rumgetüftelt, ich kann nicht glauben, dass du so einen einfachen Code heraus bekommst, ich hingegen dieses Monstrum:

Ich habe noch mal ein kleines Testprogramm geschrieben.
Es ist wirklich so, dass Dein 'Monstrum' exakt die selbe Verteilung liefert wie meine Verteilung class Cubic<>. Bei der identischen Zufallszahl liefern unsere Algorithmen zwar unterschiedliche Radien, aber bei näherer Betrachtung muss das auch so sein, aber die Verteilung ist die selbe.
Also Du hast Dich nicht verrechnet.

SeppJ schrieb:

Die Grundzüge wurden hier im Thread schon erklärt oder Verlinkt, ...

Das kann ich in keiner Weise nachvollziehen. Kannst Du da noch ein wenig konkreter werden, wie die Ausgangsformel oder Idee ist, die dann zu Deiner Lösung führt. Das würde mich interessieren.

Gruß
Werner

SeppJ

Werner Salomon schrieb:

SeppJ schrieb:

Die Grundzüge wurden hier im Thread schon erklärt oder Verlinkt, ...

Das kann ich in keiner Weise nachvollziehen. Kannst Du da noch ein wenig konkreter werden, wie die Ausgangsformel oder Idee ist, die dann zu Deiner Lösung führt. Das würde mich interessieren.

Im Nachhinein gesehen war's eher umständlich, auf otzes Trick hätte ich kommen müssen. Mein's ist im Prinzip die Inversionsmethode direkt mit dem Kopf voran durch die Betonwand gerechnet, anstatt die Tür zu nehmen .

Die Verteilung die ich haben möchte ist

P(R) = R^2/(R2^3/3 - R1^3/3) * Theta(R-R1) * Theta(R2-R)

Ich habe eine GLeichverteilung 0-1 die ich dazu umformen möchte. Wenn ich nun für die Inversionsmethode die kumulierte Wahrscheinlichkeitsverteilung zu berechnen, bekomme ich eine Funktion die bis R1 0 ist und dann wie x^3 bis 1 bei R2 steigt. Leider nicht direkt umkehrbar, weil mehrdeutig bei 0.

Und an dieser Stelle hätte ich es wie otze machen sollen . Was ich stattdessen gemacht habe: Das Problem, wie ich es sah, war, dass die Verteilung bei R1 nicht 0 ist. Wenn's einfach nur quadratisch von R1 bis R2 ginge und P(R1)=0, wäre das Problem einfach. Aber es ist im Prinzip ein "Block" von Gleichverteilung unter der quadratischen Verteilung. Daher habe ich gesagt, ich nehme aus meinen gleichverteilten Zufallszahlen einen Block eben dieser Größe (nämlich P(R1)*(R2-R1)) heraus und strecke diesen, so dass er diesem Block unter meiner Parabel entspricht. Das ist der erste if-Zweig bei mir.

Den Rest der gleichverteilten Zufallszahlen habe ich dann wieder gestreckt, so dass ich wieder eine Gleichvertelung habe, so dass ich die Inversionsmethode nach Schulbuch anwenden kann (Zeile 27). Nun hatte ich nur noch den parabelförmigen Teil von P(R) um den ich mich kümmern musste. Auf den Nullpjunkt verschoben und die Inversionsmethode direkt angewendet. Doch dummerweise haben sich in der Rechnung bis dahin derart viele Vorfaktoren angesammelt, dass das dann doch nicht mehr so einfach war, wie der Plan vorsah. Als kumulierte Wahrscheinlichkeitsfunktion erhielt ich

F(R) = (r^2 (r + 3 R1))/((R1 - R2)^2 (2 R1 + R2))

Nun, und das habe ich dann umgekehrt . Und während draußen Gewitterstürme tobten und der strenge Geruch von Formaldehyd von den eingemachten Organen im Regal aufstieg, erwachte so in meinem Labor beim kalten Lichte der Gasüberschlagslampen das MONSTER . Muahuahuahuahua!!! Sie haben mich wahnsinnig genannt!! Ich werde ihnen zeigen, wer hier wahnsinnig ist!!!!!

Werner Salomon

SeppJ schrieb:

Die Verteilung die ich haben möchte ist
P(R) = R^2/(R2^3/3 - R1^3/3) * Theta(R-R1) * Theta(R2-R)

wo kommt diese Gleichung her. Und was ist die Funktion Theta()? In Statistik bin ich leider ziemlich schwach.

Gruß
Werner

SeppJ

Heaviside-Stufenfunktion .

Die Verteilungsfunktion ergibt sich dadurch, dass die Kugelflächen wie R^2 gehen, die Verteilung der Radien die wir ziehen, muss daher auch proportional zu R^2 sein, damit man hinterher eine gleichmäßige Verteilung im Raum bekommt. Die Werte sollen durch R1 und R2 eingeschränkt sein, daher die Stufenfunktionen. Und die Gesamtwahrscheinlichkeit muss 1 sein, das ist die Normierung unter dem Bruch.

*Rewind*

SeppJ schrieb:

Nun, und das habe ich dann umgekehrt . Und während draußen Gewitterstürme tobten und der strenge Geruch von Formaldehyd von den eingemachten Organen im Regal aufstieg, erwachte so in meinem Labor beim kalten Lichte der Gasüberschlagslampen das MONSTER . Muahuahuahuahua!!! Sie haben mich wahnsinnig genannt!! Ich werde ihnen zeigen, wer hier wahnsinnig ist!!!!!

Zu viel Sinn City geguckt?

volkard

SeppJ schrieb:

Ob sich die Rechnung nun lohnt sei dahingestellt. Bei der Wegwerfmethode kann man höchstens (4/3 r^3)/(2r)3 = 1/6 der Zufallszahlen benutzen, wenn man eine Vollkugel hat.

Da ist Π verlorengegangen.
(4/3 Π r^3)/(2r)3 = 1/6 Π = 52%

Im Fall, daß die innere Kugel den halben Radius hat, sind's noch 46%, wie pumuckl auf Seite 2 festgestellt hatte.

otze

Werner Salomon schrieb:

was dann kommt habe ich zwar nicht verstanden - insbesondere die Ausgangsgleichung

Okay, ich geb mir mal mehr Muehe

Eine vereinfachende Definition zuerst, ich bezeichne die Oberflaeche einer Kugel mit Radius r mit:
A(r) = 4* pi * r^2

was ich nun machen will ist, wie du bereits sagtest: die Auftrittswahrscheinlichkeit eines r muss proportional zu A(r) sein. Da wir zugleich die Nebenbedingung haben, dass das \int_a^1 p(r) dr = 1 sein muss, ergibt sich der Proportionalitaetsfaktor als die Normalisierungskonstante.

p(r) = A(r) / \int_a^1 A(r) dr = A(r)/ (4/3*pi - 4/3*pi*a^3)

und dann nurnoch A(r) einsetzen.

Jetzt bin ich aber auf deinen Gedankengang gespannt. Du scheinst irgendwie straight forward gewusst zu haben, wie die r Verteilt sein muessen.

brotbernd

volkard schrieb:

SeppJ schrieb:

Ob sich die Rechnung nun lohnt sei dahingestellt. Bei der Wegwerfmethode kann man höchstens (4/3 r^3)/(2r)3 = 1/6 der Zufallszahlen benutzen, wenn man eine Vollkugel hat.

Da ist Π verlorengegangen.
(4/3 Π r^3)/(2r)3 = 1/6 Π = 52%

Im Fall, daß die innere Kugel den halben Radius hat, sind's noch 46%, wie pumuckl auf Seite 2 festgestellt hatte.

Zieht man in einer Würfelschale, die die Kugelschale enthält sind es auch wieder 52%. Man kann also auch bei sehr dünnen Schalen genauso viele Zufallszahlen nutzen.

Dobi

Wenn die Kugelschale sehr dünn wird, wird man auch bei der Würfelschalenmethode ziemlich viele Zahlen wegwerfen müssen. Man könnte so eine dünne Kugelschale dann mit Würfeln approximieren. Ist dann dann aber vermutlich kaum noch einfacher als die hier diskutierte elegante/komplizierte Lösung.

brotbernd

Dobi schrieb:

Wenn die Kugelschale sehr dünn wird, wird man auch bei der Würfelschalenmethode ziemlich viele Zahlen wegwerfen müssen. Man könnte so eine dünne Kugelschale dann mit Würfeln approximieren. Ist dann dann aber vermutlich kaum noch einfacher als die hier diskutierte elegante/komplizierte Lösung.

Es ist immer 0,52.
4/3 * pi * (r1^3 - r2^3) / (8 * (r1^3 - r2^3) = 1/6 pi

Edit: Falscher Fehler. Der innere Würfel ist natürlich innerhalb der inneren Kugel. Das war also quatsch...

brotbernd

...

otze

Was dem Thread jetzt fehlt ist die Verallgemeinerung auf Hyperkugeln.

krümelkacker

Der Trick mit den Gaußverteilungen für eine Gleichverteilung auf einer Kugeloberfläche klappt immer -- unabhängig von der Dimension n. Für den Radius r benötigt man dann aber eine W'keitsdichte P(r), die von n abhängt:
$P(r) \propto r^{n-1}$

krümelkacker

Ich sehe gerade, dass Boost sogar schon eine Gleichverteilung über einer Kugeloberfläche anbietet mit demselben Trick, den ich auch jetzt schon mehrfach vorschlug: boost::random::uniform_on_sphere.

Das mit dem Radius ist auch ganz einfach: Für den n-dimensionalen Raum zwischen zwei Kugeloberflächen mit Radius a und b: Einfach eine Zahl aus dem Intervall [a^n, b^n] per Gleichverteilung auswählen und daraus die n-te Wurzel ziehen. Fertig.

Dobi

Kann man sich den ganzen Aufwand mit der Gleichverteilung bei der Polarkoordinatendarstellung nicht einfach sparen, indem mal einen Zufallspunkt im Würfel erzeugt, und dessen Abstand zum Mittelpunkt dann mit der von krümelkacker beschriebenen Methode für den Radius korrigiert?

krümelkacker

Dobi schrieb:

Kann man sich den ganzen Aufwand mit der Gleichverteilung bei der Polarkoordinatendarstellung nicht einfach sparen, indem mal einen Zufallspunkt im Würfel erzeugt, und dessen Abstand zum Mittelpunkt dann mit der von krümelkacker beschriebenen Methode für den Radius korrigiert?

Nein. Das, was Du dafür brauchst ist eine Normalverteilung für jede Achse, keine Gleichverteilung im Würfel. Sonst bekommst in den Richtungen, wo der Würfel seine Ecken hat, Häufungen. Aber scheinbar muss ich das 5mal sagen, dass das mit der Normalverteilung elegant einfach wird.... :p

Dobi

Ah, OK, stimmt.
Man könnste es dann nur so machen, dass man alle Punkte wegwirft, die ausserhalb der Kugel liegen. Damit würde man dann auch bei dünnen Kugelschalen nicht mehr Verlust haben als bei einer Vollkugel.

Vermutlich musst du es für mich auch noch ein 6stes Mal sagen, weil ich noch nicht verstanden hab, wie das funktioniert, und deshalb nach was mathematisch einfacherem suche.

Werner Salomon

otze schrieb:

Jetzt bin ich aber auf deinen Gedankengang gespannt. Du scheinst irgendwie straight forward gewusst zu haben, wie die r Verteilt sein muessen.

Hallo miteinander,

entschuldigt bitte, dass meine Antwort so lange auf sich warten lies, aber ich bin z.Zt. offline stark beschäftigt.

@otze Ich habe Deine Herleitung inzwischen verstanden - mit funktionierendem Latex wäre es auch ein wenig einfacher gewesen.

Mein Gedankengang war einfach ein zufälliges Volumen zwischen zwei gegebenen Volumina - also dem der inneren Kugel und dem der äußeren Kugel - zu wählen und dann daraus den Radius zu bestimmen. Das ist schon alles. Finde ich eher noch naheliegender, als Deine Herleitung.
Bei der Herleitung des Breitenwinkels (in meinem Code der Winkel 'phi' in class Sphere) bin ich genauso vorgegangen. Ich wähle eine zufällige Oberfläche und bestimmen daraus den Winkel (bzw.den sin des Winkels).

Was SeppJ dort gemacht hat, habe ich nicht entschlüsseln können.

krümelkacker schrieb:

Ich sehe gerade, dass Boost sogar schon eine Gleichverteilung über einer Kugeloberfläche anbietet mit demselben Trick, den ich auch jetzt schon mehrfach vorschlug: boost::random::uniform_on_sphere.

Das mit dem Radius ist auch ganz einfach: Für den n-dimensionalen Raum zwischen zwei Kugeloberflächen mit Radius a und b: Einfach eine Zahl aus dem Intervall [a^n, b^n] per Gleichverteilung auswählen und daraus die n-te Wurzel ziehen. Fertig.

@krümelkacker: Danke für diesen Hinweis. Ich habe in Deinem Beitrag zu spät das Wort normalverteilt gelesen, und hatte ihn daher daher zunächst ignoriert.

Kann denn jemand (nachvollziehbar!) beweisen, dass dieser Algorithmus tatsächlich zu gleich verteilten Punkten auf einer Kugelschale führt. Ich habe ein kleines Testprogramm geschrieben - rein numerisch ist das so.

Faszinierend ist, dass SeppJ, ich und auch boost::uniform_sphere zu identischen Ergebnissen kommen; obwohl in jedem Fall der gewählte Algorithmus völlig unterschiedlich ist.

Gruß
Werner

krümelkacker

Werner Salomon schrieb:

Kann denn jemand (nachvollziehbar!) beweisen, dass dieser Algorithmus tatsächlich zu gleich verteilten Punkten auf einer Kugelschale führt.

Ja. Sei N die Anzahl der Dimensionen und X_k mit k=1...N unabhängige Zufallsvariablen mit derselben Verteilung: normalverteilt mit Erwartungswert 0 und Streuung 1. Wenn wir X als Vektor [X_1 X_2 X_3 ...]^T betrachten, lässt sich aufgrund der Unabhängigkeit der einzelnen Komponenten die Wahrscheinlichkeitsdichte für X per Multiplikation angeben:
P(x) = (1/sqrt(2*pi))^N * exp(-1/2 x_1^2) * exp(-1/2 x_2^2) * ... * exp(-1/2 x_N^2)
Eines der Potenzgesetze sagt nun, dass wie das so umschreiben können
P(x) = (1/sqrt(2*pi))^N * exp(-1/2 (x_1^2+x_22+x_3^2+...+x_N2) )
Und x_1^2+x_22+x_3^2+...+x_N2 ist einfach das Quadrat der Läng des Vektors. Man sieht also, dass die Wahrscheinlichkeitsdichte nicht mehr von der Richtung des Vektors x abhängt, sondern nur noch von seiner 2-Norm. Die Wahrscheinlichkeitsverteilung auf der Kugeloberfläche nach einem Normalisierungsschritt ist also uniform.