Probleme beim Zufalls-Punkte erzeugen [gelöst]

[heinz]

Hallo,

es hat mir einfach keine Ruhe gelassen...

Mein derzeitiger Stand ist:

Code: Alles auswählen

int radX=2000;
int radY=1000;
int radZ=500;
double rad_faktorY=(double)radX/(double)radY;
double rad_faktorZ=(double)radX/(double)radZ;

while(count < 10000)
{
	double winkel=(double)(rand()%36000000)/100000.0;
	double winkelZ=(double)(rand()%36000000)/100000.0;
	double radius=rand()%radX;

	double XX=sin(winkelZ*(Pi/180.0)) * cos(winkel*(Pi/180.0)) * radius;
	double YY=( ( sin(winkelZ * (Pi/180.0)) * sin(winkel * (Pi/180.0)) ) / rad_faktorY ) * radius;
	double ZZ=( cos(winkelZ * (Pi/180.0)) / rad_faktorZ ) * radius;
}

Das Ergebnis sieht in der x/y-Achse jetzt gut aus

aber nun habe ich das Problem auf der z-Achse...


Werde das Projekt wohl vorerst aufgeben.
Meine Mathekenntnisse sind einfach nicht ausreichend...

Vielen Dank Euch allen fuer Eure Unterstuetzung.

Gruss heinz

[reox]

Kann es sein, dass das nur ein Artefakt ist, weil du die Achsen unterschiedlich skalierst?

Ansonsten schaut die Variante 3 Zufallszahlen: 2 Winkel + Radius gut aus.

[heinz]

Hallo reox,

Kann es sein, dass das nur ein Artefakt ist, weil du die Achsen unterschiedlich skalierst?

leider nicht.

Ich habe fuer die letzten Bilder

Code: Alles auswählen

int radX=2000;
int radY=2000;
int radZ=2000;

verwendet, weil man es dann besser erkennt. (Habe vergessen es dazu zu schreiben...)
Es ist egal welche Werte ich fuer x,y und z verwende, in der z-Achse ist eine starke haeufung bei nahe 0.

Gruss heinz

[reox]

Mhh interessant, also ich sehe dieses Phänomen auch:


Generiert mit:

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

N = 10000

r = 1

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

theta = np.pi * np.random.rand(N)
phi = 2 * np.pi * np.random.rand(N)
radius = r * np.random.rand(N)

x = np.sin(theta) * np.cos(phi) * radius;
y = np.sin(theta) * np.sin(phi) * radius;
z = np.cos(theta) * radius

ax.scatter(x, y, z)
plt.show()

Eine Frage wäre jetzt: Wie kann man zB ein Histogramm erzeugen aus den Koordinaten um zu schauen ob die Verteilung passt?

Also es schaut so aus, als ob die Verteilung an den Polen größer ist. Nur warum? Offenbar erzeugt die Verteilung von phi und radius das erwartete Bild. Nur theta scheint nicht das zu tun was wir erwarten.

Edit: Theta darf nur werte zwischen [0, pi] haben. Ändert aber auch nichts am Bild.

[heinz]

Eine Frage wäre jetzt: Wie kann man zB ein Histogramm erzeugen aus den Koordinaten um zu schauen ob die Verteilung passt?

Also wenn ich die Z-Koordinaten aus der erzeugten POV-Datei mal nach der Haeufigkeit aufliste, erkennt man sehr gut
wie sie in der naehe von 0 zunimmt.
Habe die Nachkommastellen mal weggelassen und dann die Haeufigkeit der Werte ermittelt.

Code: Alles auswählen

grep '^translate' ./zufall.pov |cut -f3 -d','|tr -d '>'|cut -f1 -d'.'|sort -g|uniq -c >./z-h.txt

40584
Erste Spalte=Haeufigkeit
Zweite Spalte=Z-Werte

Meintest Du sowas?

Gruss heinz

[Edit:]
Sorry, das war Quatsch.
Das zeigt nur die Anhaeufung um die Mitte.

[reox]

also was man schon sieht, ist ja dass x und y jeweils von 3 zufallszahlen abhängt aber z nur von zweien.
Und wenn man aus der Z Richtung drauf schaut, schauts ja gut aus.
Aber mir leuchtet das nicht ein warum das einen unterschied machen sollte...

[smutbert]

hab das einmal in r nachgestellt und im Histogramm ist der Unterschied deutlich zu sehen. Zum Nachdenken wieso bin ich noch nicht gekommen...

[reox]

hab das einmal in r nachgestellt und im Histogramm ist der Unterschied deutlich zu sehen. Zum Nachdenken wieso bin ich noch nicht gekommen...

Ok, 3 histogramme (dichteschätzer) sind besser als einer
Gibt es eigentlich eine Möglichkeit Kreis oder Kugelkoordinaten (oder sagen wir n-dimensionale koordinaten) in einer einzigen Zahl zusammenzufassen? Gibt es da eine Ordnung die zB Koordinaten in Rationale Zahlen umbaut?

[Lohengrin]

Code: Alles auswählen

theta = np.pi * np.random.rand(N)
phi = 2 * np.pi * np.random.rand(N)
radius = r * np.random.rand(N)

x = np.sin(theta) * np.cos(phi) * radius;
y = np.sin(theta) * np.sin(phi) * radius;
z = np.cos(theta) * radius

Also es schaut so aus, als ob die Verteilung an den Polen größer ist. Nur warum? Offenbar erzeugt die Verteilung von phi und radius das erwartete Bild. Nur theta scheint nicht das zu tun was wir erwarten.

Beachte die Fläche auf der Kugel zwischen dem theta-ten und dem (theta + epsilon)-ten Breitengrad!
Das hängt von theta ab. Je näher am Pol, desto kleiner diese Fläche.

[smutbert]

@heinz

Darf man zwischendurch einmal fragen wie die Verteilung aussehen soll?
Die größere Dichte in der Mitte ist gewünscht?
(Entschuldige, falls es schon im Thread steht, aber ich hab nichts gefunden.)

[Lohengrin]

Darf man zwischendurch einmal fragen wie die Verteilung aussehen soll?
Die größere Dichte in der Mitte ist gewünscht?

Ich denke, dass die Verteilung über eine Kugelsphäre gleich sein soll. Das mit dem Radius soll unabhängig davon sein.

Mit theta=pi*Zufall (Zufall ist gleichverteilt auf [0,1]) haut das nicht hin, weil ein Breitengrad die Länge radius*2*pi*sin(theta) hat. Die Verteilung der theta muss so sein, dass das bei gleichem Radius eine gleiche Verteilung auf der Kugelsphäre ergibt.
Nun ist in jedem Punkt der Kugelsphäre der Winkel zwischen dem Meridian und der Ebene parallel zum Äquator ebenfalls theta, also bei Veränderung des theta das Verhältnis von der Änderung in z-Richttung zur Änderung des theta sin(theta). Wenn man auf der z-Achse gleichverteilt, und das dann auf den Meridian abbildet, passt es. Also ist theta=acos(2*Zufall - 1) passend.

[heinz]

Danke erstmal an diejenigen die sich noch gedanken ueber diese "Sache" machen.

Darf man zwischendurch einmal fragen wie die Verteilung aussehen soll?
Die größere Dichte in der Mitte ist gewünscht?

Ich wäre schon zufrieden wenn die Betrachtung von oben

und von der Seite

gleich aussehen wuerde, also keine "seltsamen" Anhaufungen ausser die im Zentrum.

MMn. ist die groessere Dichte in der Mitte eine andere Baustelle um die ich mich noch nicht gekuemmert habe.
Es wurde ja schon von hikaru geschrieben:
Für eine gleichmäßige Dichte in der Fläche musst du dafür sorgen, dass deine Punkte gleichmäßig auf dem Quadrat des Radius' verteilt sind.
Nach meinem Verstaendnis, je groesser der Radius um so groesser die Wahrscheinlichkeit fuer einen gesetzten Punkt. (Quadratisch zunehmend)
(Vlt. eine rand-Funktion die (Quadratisch) mehr grosse Zahlen (Radien) liefert als kleine.)
Ich denke/hoffe, das bekomme ich hin. Wollte nur erst mal die eine "Baustelle" beenden, bevor ich mit der naechsten beginne.

Ich denke, dass die Verteilung über eine Kugelsphäre gleich sein soll. Das mit dem Radius soll unabhängig davon sein.

Am Ende waere ich froh, beides zu haben (an-/abschaltbar), einmal mit haeufung bei kleinen Radien und einmal ohne.
Aber Du hast recht, die 3 Radien (x,y,z) sollen einzeln bestimmbar sein.

Viele Gruesse,
heinz

[smutbert]

So habe ich es jetzt probiert, allerdings wieder in r, weil es zu lange dauert bis ich in C etwas lauffähiges hinbekomme (falls ich es überhaupt schaffe):

Code: Alles auswählen

library("scatterplot3d")
library("car")

N=50
rx=1
ry=2
rz=1.5

phi=2*pi*runif(N)
theta=acos(2*runif(N)-1)

x=rx*cos(phi)*sin(theta)
y=ry*sin(phi)*sin(theta)
z=rz*cos(theta)

scatter3d(x, y, z, surface=FALSE)

Damit berechne ich, wie Logengrinn es geschrieben hat, Punkte, die auf der Oberfläche eines Ellipsoids gleich(mäßig) verteilt sind.
runif erzeugt gleich verteilte Zufallszahlen zwischen 0 und 1, das heißt da kommen Vektoren mit der Länge N (50) heraus. Mit dem letzten Befehl wird das Ergebnis geplottet, allerdings erscheint es wieder als Kugel, weil die Achsen einzeln automatisch skaliert werden (und ich habe es auf die Schnelle nicht geschafft die Skalierung fix vorzugeben).



Um nun das ganz Volumen gleichmäßig mit Punkten zu füllen, führe ich das ganze in einer Schleife aus und zwar mit einem Skalierungsfaktor, den ich mit den Werten im Beispiel in Hundertstelschritten von 1/100 bis 1 durchlaufen lasse und plotte jedes Mal die erzeugten Punkte.
Mit N=10 und den Nr=100 Skalierungsfaktoren sollten am Ende 1000 Punkte herauskommen.

Code: Alles auswählen

library("rgl")

N=10
Nr=100
rx=1
ry=2
rz=1.5

for(i in 1:Nr) {
	scale=i/Nr
	phi=2*pi*runif(N)
	theta=acos(2*runif(N)-1)

	x=scale*rx*cos(phi)*sin(theta)
	y=scale*ry*sin(phi)*sin(theta)
	z=scale*rz*cos(theta)
	
	plot3d(x, y, z, add=TRUE)
}

Es gibt keine bevorzugte Richtung mehr. Innen häufen sich die Punkte, weil bei den früheren Schleifendurchläufen trotz kleinerem Skalierungsfaktor/Radius gleich viele Punkte gibt wie bei den späteren Schleifendurchgängen. Ohne es jetzt genau durchdacht zu haben, wird man die Häufung los, indem man das Quadrat des Skalierungsfaktors in die Zahl der Zufallszahlen pro Schleifendurchgang einfließen lässt

Code: Alles auswählen

library("rgl")

N=20
rx=1
ry=2
rz=1.5

Nr=100

for(i in 1:Nr) {
	scale=i/Nr
	phi=2*pi*runif(N*scale^2)
	theta=acos(2*runif(N*scale^2)-1)

	plot3d(scale*rx*cos(phi)*sin(theta), scale*ry*sin(phi)*sin(theta), scale*rz*cos(theta), add=TRUE)
}

(das mit N*scale^2 als Anzahl der Zufallszahlen funktioniert in r, anderswo müsste man diese Zahl wohl auf einen Integerwert runden)


In weiterer Folge könnte man die Abhängigkeit der Dichte vom Radius beliebig ändern, zum Beispiel normalverteilte Zufallszahlen erzeugen und die dann als Skalierungsfaktoren verwenden oder ähnliches.

[Lohengrin]

Nach meinem Verstaendnis, je groesser der Radius um so groesser die Wahrscheinlichkeit fuer einen gesetzten Punkt. (Quadratisch zunehmend)
(Vlt. eine rand-Funktion die (Quadratisch) mehr grosse Zahlen (Radien) liefert als kleine.)

Ja, du brauchst eine Zufallsvariable, die eine Verteilung hat, wo die Wahrscheinlichkeitsdichte quadratisch wächst.
So etwas kannst du aus einer Gleichverteilung auf [0,1] erzeugen, indem zu das mit etwas Passendem auf [0,R] abbildest.

Ein Fachbegriff dazu ist Bildmaß. Eine Verteilung (zB Gleichverteilung auf [0,1]) wird durch die Zufallsvariable auf einen anderen Bereich (zB [0,R]) übertragen.
Auf diese Tour kann man auch Verteilungen erzeugen, die so sind, wie etwas, das man gerade gemessen hat. Wenn man zB tausend Messwerte hat, kann man sich 500 davon zufällig auswählen, damit die passende Zufallsvariable bauen, und dann mit Zufallswerten, die so sind wie diese 500, an den anderen 500 herumspielen. Aber man darf da nicht zu viel hineininterpretieren. Wenn man Pech hat, waren die einen 500 und die anderen 500 doch nicht voneinander unabhängig, man publiziert eine Sensation, und dann korreliert das doch nur mit dem Fahrplan der Straßenbahn (pünktlich alle 10 Minuten hat da was vibriert).

Wenn du etwas haben willst, wo "größeres aus [0,1] auf größeres in [0,R]" gilt, gibt es nur eine Lösung.
Beim Kreis (2-dimensional) geht das mit Wurzel.
Beweis: 0 wird auf 0 abgebildet, 1 auf R. Irgend ein x wird auf r abgebildet, also das Stück [0,x] auf das Stück [0,r]. Das Verhältnis der Kreisflächen mit r und R ist (r/R)^2. Im Kreis mit Radius r soll (r/R)^2 der Punkte liegen. Weil bei Gleichverteilung auf [0,1] das Stück von [0,(r/R)^2] das Gewicht (r/R)^2 hat, muss x=(r/R)^2 auf r abgebildet werden. Nach r aufgelöst ergibt das r=R*Wurzel(x).

Bei der Kugel (3-dimensional) geht es genauso. Das Verhältnis der Volumen mit Radius r und R ist (r/R)^3. Also muss x=(r/R)^3 auf r abgebildet werden. Nach r aufgelöst ergibt das r=R*Kubikwurzel(x).

Ich denke, dass die Verteilung über eine Kugelsphäre gleich sein soll. Das mit dem Radius soll unabhängig davon sein.

Am Ende waere ich froh, beides zu haben (an-/abschaltbar), einmal mit haeufung bei kleinen Radien und einmal ohne.
Aber Du hast recht, die 3 Radien (x,y,z) sollen einzeln bestimmbar sein.

Was für drei Radien? Meinst du drei Radien vom Ellipsoid?
Das mit acos für theta geht nur, wenn das eine Kugel ist. Wenn du die Kugel verzerrst, dann hast du in jedem Ellipsoidsektor so viele Punkte, wie es der Raumwinkel der Elipsoidsphäre verlangt, und das ist dann je nach Richtung verschieden großes Volumen. Du müsstest da je nach Raumwinkel unterschiedlich viel reinwerfen.

Ich empfehle immernoch Montecarlo. Gleichverteilung im Raum geht mit kartesischen Koordinaten ganz leicht. Und was draußen liegt, ist leicht zu erkennen.

Damit berechne ich, wie Logengrinn es geschrieben hat, Punkte, die auf der Oberfläche eines Ellipsoids gleich(mäßig) verteilt sind.

Die sind auf den Raumwinkeln gleichmäßig verteilt. Auf der Fläche sind sie es nicht, weil du die Kugel verzerrt hast.
Wenn du es geschafft hast, Gleichverteilung auf einer Ellipsoidsphäre zu bauen, passt es im Volumen wieder nicht, weil je nach Richtung die Schalen unterschiedlich dicht liegen.

Dass du theta und phi vertauscht hast, könnte verwirren.

[smutbert]

Danke für die Erklärungen. Ich habe mir beim Übergang von der Kugel zu dem Ellipsoid zu wenig (eigentlich gar nichts) gedacht. (zumindest habe ich jetzt das phi/theta-Chaos aufgeräumt)

[heinz]

@Lohengrin
Danke nochmals fuer die ausfuehrliche Erklaerung!

Was für drei Radien? Meinst du drei Radien vom Ellipsoid?

rX=2000; rY=1000; rZ=500

Ich hoffe so wird es "klarer" was ich meine...

Ich empfehle immernoch Montecarlo.

Und ein paar Posts vorher:

Kreis ...Du nimmst solange zwei Zufallszahlen z0 und z1 zwischen 0 und 1, bis z0^2+z1^2 < 1
Ellipsoid geht genauso. Da wiederholst du, bis du (z0/a)^2+(z1/b)^2+(z2/c)^2 < 1 hast.

das wollte ich versuchen umzusetzen.
Mein Ziel ist ein Ellipsoid.
Der Kreis scheint mir verstaendlich.
Allerdings scheitere ich daran, zu verstehen was im Zusammenhang mit dem Ellipsoid a,b und c ist.
Sind das die max. Radius-Werte? Also z.B. a=2000; b=1000; c=500;?
Falls ja, muessten die Radien dann nicht auch in der Kreis-Formel auftauchen?
z.B. (z0/a)^2+(z1/b)^2 < 1

Gruss, heinz

[Lohengrin]

Kreis ...Du nimmst solange zwei Zufallszahlen z0 und z1 zwischen 0 und 1, bis z0^2+z1^2 < 1
Ellipsoid geht genauso. Da wiederholst du, bis du (z0/a)^2+(z1/b)^2+(z2/c)^2 < 1 hast.

das wollte ich versuchen umzusetzen.
Mein Ziel ist ein Ellipsoid.
Der Kreis scheint mir verstaendlich.
Allerdings scheitere ich daran, zu verstehen was im Zusammenhang mit dem Ellipsoid a,b und c ist.
Sind das die max. Radius-Werte? Also z.B. a=2000; b=1000; c=500;?

Das habe ich falsch, mindestens verwirrend, erklärt.
Beim Kreis habe ich z0 und z1 aus dem Bereich von 0 bis 1 genommen.
Beim Ellipsoid habe ich z0 aus dem Bereich von -a bis a, z1 aus dem Bereich von -b bis b und z2 aus dem Bereich von -c bis c gemeint.

Du kannst einen Ellipsoid erzeugen, wenn du die Einheits-Kugel streckst. Streckst du in alle drei Richtungen gleich, also a=b und b=c, dann wird es eine Kugel.
Wenn du etwas in den Raum streust, und wissen willst, ob es im Ellipsoid ist, dann machst du die Streckung wieder rückgängig, und schaust, ob das dann in der Einheitskugel ist.

Falls ja, muessten die Radien dann nicht auch in der Kreis-Formel auftauchen?
z.B. (z0/a)^2+(z1/b)^2 < 1

Beim Kreis ist a=b, wird da üblicherweise R genannt. (Das ist jetzt genau der Mathematikerwitz "Beweis durch Verwirrung": Sei a gleich b, wir nennen es c.) Aus (z0/a)^2+(z1/b)^2 < 1 wird dann (z0/R)^2+(z1/R)^2 < 1, mit "mal R^2" äquivalent umgeformt ergibt das z0^2+z1^2 < R^2, also "(z0,z1) liegt in Kugel mit Radius R" gemäß Pythagoras (Beachte das Dreieck mit den Eckpunkten (0,0), (z0,0), (z0,z1). Der rechte Winkel ist bei (z0,0), die lange Seite ist von (0,0) nach (z0,z1), Länge R.).

Um einen gleichmäßig gefüllten Ellipsoid zu bauen, kannst du z0,z1,z2 aus dem Bereich von -1 bis 1 wählen, und mitnehmen, wenn z0^2+z1^2+z2^2 < 1. Dann hast du eine gefüllte Einheitskugel. Das streckst du dann zum Ellipsoid, also (z0,z1,z2) wird zu (a*z0, b*z1, c*z2).

Du kannst aber auch z0 aus dem Bereich von -a bis a, z1 aus dem Bereich von -b bis -b und z2 aus dem Bereich von -c bis c nehmen. Dann hast du erst einmal einen gleichmäßig gefüllten Quader. Du nimmst den Punkt mit, wenn er im Ellipsoid liegt, also wenn (z0/a)^2 + (z1/b)^2 + (z2/c)^2 < 1.

[heinz]

Ich hoffe, ich habe es jetzt kapiert.

Der Code sieht jetzt so aus:
40591
Ist das so wie Du es meintest?
(Ups: > mit < verwechselt...)

Die Werte, die dabei entstehen:
40590

Und aussehen tut es jetzt so:

(Ist wieder ein animiertes .gif Rotiert um die x- und y-Achse
Die x-, y- und z-Radius-Werte sind alle gleich damit eine Kugel entsteht.
Die Kamera blickt von oben drauf.)

Ersetze ich:

Code: Alles auswählen

while(rX*rX + rY*rY + rZ*rZ > 1.0);

durch

Code: Alles auswählen

while(rX*rX + rY*rY + rZ*rZ < 1.0);

Sieht es so aus:



Es ist zum Maeusemelken...

Gruss, heinz

[Lohengrin]

Warum die Winkel? Lass das weg!

Code: Alles auswählen

do
        {
                rX=(double)(rand()%20000000 - 10000000)/10000000.0;
                rY=(double)(rand()%20000000 - 10000000)/10000000.0;
                rZ=(double)(rand()%20000000 - 10000000)/10000000.0;
        }
        while(rX*rX + rY*rY + rZ*rZ > 1.0);

macht solange einen Punkt (rX,rY,rZ) in den Würfel [-1,1]x[-1,1]x[-1,1], bis (rX,rY,rZ) in der Einheitskugel liegt.
Das steckst du dann.

Code: Alles auswählen

XX=radX*rX;
YY=radY*rY;
ZZ=radZ*rZ;

Das Ergebnis ist ein Punkt (XX,YY,ZZ) im Ellipsoid

[heinz]

Lieber Lohengrin ,

ES FUNKTIONIERT!

Warum die Winkel? Lass das weg!

Sorry, das ich so lange gebraucht habe um es zu kapieren...

Vielen herzlichen Dank fuer Deine Geduld und Deine ausdauernde Hilfe!

Auch nochmals vielen Dank an alle anderen, die sich damit beschaeftigt haben.

Gruss, heinz

Falls es sonst noch jemanden interessiert, hier ist die endgueltige Loesung:
40592

[reox]

Warum die Winkel? Lass das weg!

Ja die Lösung hatten wir ja auch schon für 2D. Interessant wäre es jetzt so eine Lösung zu machen ohne das man Werte wegwirft

[heinz]

Interessant wäre es jetzt so eine Lösung zu machen ohne das man Werte wegwirft

Interessieren wuerde mich das natuerlich auch, allerdings befuerchte ich, dass ich keine
Ahnung davon haette, warum es am Ende funktioniert...
Ich muss gestehen, etwas neidisch auf Euer mathematisches Verstaendnis zu sein...
Auf jeden Fall bin ich sehr Dankbar, hier (wie schon oft) eine Loesung fuer mein Problem gefunden zu haben.
Falls Ihr Euch hier weiterhin mathematische Formeln um die Ohren hauen moechtet, werde ich natuerlich
mit Begeisterung mitlesen und versuchen etwas daraus zu lernen...


Mit der derzeitigen Version habe ich mal etwas herumgespielt:
Das Programm ist in C geschrieben.
Startbedingung ist:
20000 Stellen in einer ca. 2 MB. grossen .txt Datei, durch die berechneten Werte zu ersetzen.

Code: Alles auswählen

translate <[>

wird zu

Code: Alles auswählen

translate <-91.9506 , -9.98181 , -20.8027>

Bei der Berechnung in 2D werden im Schnitt 5470 erzeugte Werte (rX, rY) verworfen.
Bei der Berechnung in 3D werden im Schnitt 18200 erzeugte Werte (rX, rY, rZ) verworfen.
Die Laufzeit des Programms liegt im Schnitt zwischen 0.08 bis 0.1 sek.
Das Ergebnis sieht immer perfekt gleichverteilt aus.

Gruss, heinz

[Lohengrin]

20000 Stellen
...
Bei der Berechnung in 2D werden im Schnitt 5470 erzeugte Werte (rX, rY) verworfen.
Bei der Berechnung in 3D werden im Schnitt 18200 erzeugte Werte (rX, rY, rZ) verworfen.

Für 2D:
Der Kreis mit Radius 1 hat die Fläche pi*1^2=pi. Das Quadrat [-1,1]x[-1,1] hat die Fläche 2*2=4. Es werden die weggeworfen, die nicht im Kreis sind. Für 20000 im Kreis braucht man 4/pi mal so viele Schüsse ins Quadrat, also etwa 25500, macht etwa 5500 zu erwartende Verworfene. 5470 passt.
Für 3D:
Die Kugel mit Radius 1 hat Volumen 4/3*pi*1^3=4/3*pi (Achtung! Gerechnet wird von links nach rechts. (4/3)*pi und 4/(3*pi) ist nicht dasselbe.). Der Würfel [-1,1]x[-1,1]x[-1,1] hat Volumen 2^3=8. Es werden die weggeworfen, die nicht in der Kugel sind. Für 20000 in der Kugel braucht man 8/(4/3*pi) Schüsse in den Würfel, also etwa 38200, macht etwa 18200 zu erwartende Verworfene. 18200 passt.

Wozu will hier jemand eine Lösung haben, wo man keine Werte verliert? Es geht bei 3D um "die Hälfte wirft man weg". Dafür würde ich keinen Tamtam machen.
Früher, als man noch nicht mal so einfach ein paar Millionen mal schießen konnte, hat man solche Probleme umgeformt, und kam dann oft nur bis zu etwas, über das es Tabellen gab. Dabei hatte man dann oft immernoch einen Approximationsfehler. Das ist heute gar nicht mehr nötig, und bei Approximationsfehler ist Montecarlo sogar besser.
Ich verstehe ja, dass das die Leute ärgert, die noch gelernt haben, wie man das umformen muss, um die Tabelle benutzen zu können. Die Generation, die noch keine Computer hatte, ist heute schon lange im Ruhestand. Es ist jetzt die nächste Generation, die noch etwas gegen Montecarlo hat, weil deren Lehrer ihnen erzählt haben, dass das Banausentum sei.
Mir hat mal ein uralter Mathematiker erzählt, dass die Banausen bevor es Computer gab mit 80g-Papier, Schere und Feinwaage integriert hätten.

[reox]

Wozu will hier jemand eine Lösung haben, wo man keine Werte verliert?

Weil es vielleicht interessant ist genau dieses Problem zu lösen

[hikaru]

Wozu will hier jemand eine Lösung haben, wo man keine Werte verliert? Es geht bei 3D um "die Hälfte wirft man weg". Dafür würde ich keinen Tamtam machen.

Weil Montecarlo "Brute Force", also unelegant ist, wie mein (natürlich hinkender) Vergleich mit dem Paintballgewehr und der Zielscheibe zeigen sollte.
Den Algorithmus nutzt du nur, weil er mit moderner Rechentechnik vermeintlich billig ist. Das mag isoliert betrachtet richtig sein, aber weil alle so denken haben wir heute haufenweise Software-Bloat für den wir gigabyteweise RAM brauchen, nur um Briefe zu schreiben.
Ich hatte Spaß am Suchen nach einer Lösung, die nicht auf Quadrate schießt, wenn man Kreise treffen will. U.a. dank deiner Erklärungen habe ich auch verstanden, wo die Probleme der vorgestellten Implementierungen liegen. Ich glaube qualitativ verstanden zu haben, was man tun müsste, um den ursprünglichen Winkelansatz zum Ziel zu führen. Leider kann ich keine Lösung präsentieren, da ich mich seit jeher mit der mathematischen Formulierung von Lösungen schwer tue.

@heinz:
Mich würde abschließend noch interessieren, wozu du das brauchst. Als es noch um Kreise und Ellipsen ging, dachte ich, du willst Galaxien für dein Weltraumspiel erzeugen. Aber wozu sind die Kugeln? Kugelsternhaufen?