Algorithm 678: BTPEC: sampling from the binomial distribution

[Sephiroth]

Es geht um diesen Algorithmus: http://doi.acm.org/10.1145/76909.76916 und hier der Fortran Code (http://www.netlib.org/toms/678).

Der ist in R als rbinom implementiert (http://www.koders.com/c/fid9977D34D440E61A6D92398F26ED3370F01375BAC.aspx?s=Chebyshev) und ich würde ihn gerne verstehen, da ich ihn in Java bräuchte.

Kann den mir einer erklären? Wie erzeuge ich eine binomialverteilte Zufallszahl mit gegebener Wahrscheinlichkeit p?

[Pinoccio]

Gibt es doch auch in JAVA (http://www.koders.com/java/fid198D7715E59AA6DD5486517A0DAFE2450FB38C44.aspx?s=Chebyshev)!? (runterscrollen) (Man mekr aber leider an dem schlechten Stil, daß es einen Maschinenübersetzung ist. Grausam!)

Das ganze sieht recht verworren (aka optimiert) aus, für meinen Geschmack nach einer speziell angepasten Approximation (also anhand der exakten Lösung mit definierten Fehlergrenzen Interpolations-Funktionen gebaut. Wie eben alle komplexeren Sachen ja auch intern in den Prozessoren so gerechnet werden, z.B. exp und sin. Und Binomialverteilung ist halt durch den Binomialkoeffizenten auf direktem Weg ziemlich langsam und numerisch evtl. auch nicht stabil.).

Wie erzeuge ich eine binomialverteilte Zufallszahl mit gegebener Wahrscheinlichkeit p?Wie meinst du das?
C N : NUMBER OF BERNOULLI TRIALS (INPUT)
C PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
C ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
C JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)

mfg

/Disclaimer: Ich bin kein Numeriker.

[pest]

Inversionsmethode, Verwerfungsmethode (http://de.wikipedia.org/wiki/Verwerfungsmethode)<-macht dein Programm rbinom.c


/Disclaimer: Ich bin kein Numeriker.
ich schon :uhippie:

[Pinoccio]

Inversionsmethode, Verwerfungsmethode (http://de.wikipedia.org/wiki/Verwerfungsmethode)<-macht dein Programm rbinom.cAha. Wieder was gelernt.
Welche Verteilung G (nach Wikipeida-Namen) wird in rbinom verwendet?ich schon :uhippie:<Dein genauer Studiengang> - das ist doch eher Hilfsschul-Mathe! :ulove3:

mfg

[pest]

Aha. Wieder was gelernt.
Welche Verteilung G (nach Wikipeida-Namen) wird in rbinom verwendet?


Also so weit ich das überblicke, wird für einen Erwartungswert < 30 die Inversionsmethode angewendet, und ansonsten ist die Hilfsverteilung ne Mischung aus Rechteck-und Dreieck-Verteilung.

edit: das machen sie höchstwahrscheinlich nur aus geschwindigkeitsgründen, die inversionsmethode ist eigentlich in der lage, die verteilung genau zu approximieren, wäre aber elendig langsam
weil man bei jeder auswertung (realisierung einer zg) das integral der dichtefunktion bis zu einem punkt bestimmen muss.


<Dein genauer Studiengang> - das ist doch eher Hilfsschul-Mathe!


so unrecht hast du nicht-wenn du die abschlussarbeiten der anderen siehst :(
aber wenn wir gerade dabei sind, ich habe die regenwaldmodellierung revolutioniert :D - deshalb kenne ich das zeug hier auch...brauchte die polar-methode

[Sephiroth]

Gibt es doch auch in JAVA (http://www.koders.com/java/fid198D7715E59AA6DD5486517A0DAFE2450FB38C44.aspx?s=Chebyshev)!? (runterscrollen) (Man mekr aber leider an dem schlechten Stil, daß es einen Maschinenübersetzung ist. Grausam!)

hehe, danke, wobei ich das umschreiben schon hinbekommen hätte

Inversionsmethode, Verwerfungsmethode (http://de.wikipedia.org/wiki/Verwerfungsmethode)<-macht dein Programm rbinom.c

aha :uroll:

Also so weit ich das überblicke, wird für einen Erwartungswert < 30 die Inversionsmethode angewendet
wollt ich auch grad sagen und das tät mir auch erstmal reichen :D

Aber was benutzen die für G bei np<30? Im Code müsste das qn sein (http://www.koders.com/c/fid9977D34D440E61A6D92398F26ED3370F01375BAC.aspx#L181), was (1-p)^n ist. Ist das die erzeugende Funktion (http://de.wikipedia.org/wiki/Binomialverteilung#Erzeugende_Funktion)?

function F_verteilte_Zufallszahl()
var x, u
repeat
x := G_verteilte_Zufallszahl()
u := gleichförmig_verteilte_Zufallszahl()
until u * k * g(x) < f(x)
return x
end

[pest]

Ich glaube du verwechselst was. Die Hilfsverteilung G benutzen die nur für np>30.

Die Inversionsmethode (hier für np<30) geht prinzipiell so.

- bestimme eine Zufallszahl p zwischen [0,1]
- bestimme Quantil Q_p=F^-1(p), wobei F die Verteilungsfunktion der Binomialverteilung ist.
- Q_p ist deine neue Zufallszahl

F^-1(p) wird im Abschnitt " L_np_small:" berechnet, halt en bissl schlauer als die direkte Berechnung des Integrals, was aber zum Testen ausreicht

[Sephiroth]

Alles klar, danke nochmal.

Jetzt weiß ich auch wieso die alles dran setzen, um nicht die inverse Verteilungsfunktion der Binomialverteilung (http://en.wikipedia.org/wiki/Incomplete_beta_function#Incomplete_beta_function) ausrechnen zu müssen. :o

[pest]

Jetzt weiß ich auch wieso die alles dran setzen, um nicht die inverse Verteilungsfunktion der Binomialverteilung (http://en.wikipedia.org/wiki/Incomplete_beta_function#Incomplete_beta_function) ausrechnen zu müssen. :o

das machen sie schon,
man braucht die inverse verteilungsfunktion nicht direkt, denn die verteilungsfunktion ist ja das integral über die dichtefunktion

numerisch geht man dann wie folgt vor (nur beispielhaft)


summe=0
Q_p=0

solange summe < p
summe = summe + (wert der dichtefunktion an der stelle Q_p)*schrittweite
Q_p = Q_p + Schrittweite

rückgabe: Q_p

[Pinoccio]

Evtl. ist auch noch folgendes PDF interessant: Luc Devroye, Non-Uniform Random Variate Generation, New York: Springer-Verlag, 1986. Chapter X, Discrete Univariate Distributions. (http://cg.scs.carleton.ca/~luc/chapter_ten.pdf) (ab Seite "521") Ist auch in en.WP verlinkt, da hab ich es gesehen.

@pest#9: GLaub nciht, daß du Sephiroth erklären must, wie man numerisch integriert. ;-)

mfg