C++: Ereignisse mit Wahrscheinlichkeit p auswählen

[pest]

Hallo :)

Ich möchte mit Wahrscheinlichkeit p ein zufälliges Ereigniss steuern

nun gibt es dafür verschiedene Möglichkeiten (mache das auch nicht zum Ersten mal)

und ich habe im Netz Einiges (auch viel Falsches) gesehen

ich stelle erst einmal meine Basisfunktionen vor (Verbesserungsvorschläge bzgl. unnötigen und/oder falschen Casts bzw. Fehlern sind willkommen)

wähle gleichverteilte ZV im Intervall [0,1]

static double RandDouble()
{
return (double)rand()/(double(RAND_MAX));
}


wähle gleichverteilte ZV im Intervall [0,1)
(hier bin ich mir nicht sicher ob die Intervalle um jede Zahl korrekt sind)

static double RandDouble1()
{
return (double)rand()/(double(RAND_MAX)+1.0);
}


wähle gleichverteilte ZV im Intervall (0,1)
(hier bin ich mir sicher, dass es korrekt ist (http://www.forum-3dcenter.org/vbulletin/showthread.php?t=475597))

static double RandDouble2()
{
return (double(rand())+0.5)/(double(RAND_MAX)+1.0);
}




die übliche Vorgehensweise (in Code den ich so über die Jahre gesehen habe und viele Seiten im Netz) ist nun Folgende

Ist ZV aus [0,1] < p dann Ereigniss=true, also

static bool BoolEvent(double p)
{
return (RandDouble()<p);
}


Hier gibt es aber ein Problem. Das sichere Ereigniss (also p=1.0) liefert manchmal false zurück, nun könnte man das "<" durch ein "<=" ersetzen,
holt sich damit aber das unmögliche Ereigniss mit ins Boot, und die evilness des Equal-vergleichs zwischen Gleitpunktzahlen, die hier aber kein Problem darstellen sollte

diese Variante scheint daher das gewünschte zu Erledigen
Ist ZV aus [0,1) < p dann Ereigniss=true, also

static bool BoolEvent(double p)
{
return (RandDouble1()<p);
}


oder Diese
Ist ZV aus (0,1) <= p dann Ereigniss=true, also

static bool BoolEvent(double p)
{
return (RandDouble2()<=p);
}



Welche ist korrekt, was benutzt ihr?

Gruss
pest

[Aquaschaf]

Wie wäre es mit Boost.Random? ;)

[Coda]

Wollte ich gerade auch vorschlagen. Vor allem ist das auch ordentlich getestet.

[del_4901]

Wenn es schnell sein muss:


float frand( int *seed )
{
union
{
float fres;
unsigned int ires;
};

seed[0] *= 16807;

ires = ((((unsigned int)seed)>>9 ) | 0x3f800000);
return fres - 1.0f;
}
uniform floats [0,1]

[pest]

Die Verwendung von Boost löst mein (eigentliches) Problem nicht
Im Code steht auch

/** Returns a value uniformly distributed in the range [min, max). */


ich verwende eh den Mersenne-Twister, dem Code nach zu urteilen ist meine obige Generierung auch i.O.

die Frage ist halt

Ist ZV aus [0,1) < p dann Ereigniss=true
oder
Ist ZV aus (0,1) <= p dann Ereigniss=true

[Coda]

Mit Boost kannst du dir auch Würfel bauen, die mit einer gewissne Wahrscheinlichkeit p eine bestimmte Zahle ausgeben. Das Ding hat nicht nur Random-Generatoren.

[pest]

ach ist doch overkill, ne verteilungsfunktion kann ich mir selber basteln

boost macht das

WeightType val = uniform_01<WeightType>()(urng);
WeightType sum = 0;
std::size_t result = 0;
for(typename std::vector<WeightType>::const_iterator
iter = parm._probabilities.begin(),
end = parm._probabilities.end();
iter != end; ++iter, ++result)
{
sum += *iter;
if(sum > val) {
return result;
}
}


entspricht also der Variante

Ist p > ZV aus [0,1) dann Ereigniss=true


Jetzt haben wir drei Möglichkeiten?! ;) :D

[hell_bird]

Ich denke du musst auch deinen anwendungsfall genauer spezifizieren. Es gibt keine eindeutig richtige oder eindeutig falsche Methode.

Beispiel: Nehmen wir an RAND_MAX ist 9. Dann kann jeder Wert von 0-9 mit 10% Chance vorkommen. Wenn jetzt deine Wahrscheinlichkeit 0.01 ist würde bei keiner deiner Methoden jemals das Ereignis ausgelöst werden. Das ist mindestens genauso falsch wie wenn es in 10% der Fälle ausgelöst wird.

[pest]

Es gibt keine eindeutig richtige oder eindeutig falsche Methode.


für p=1.0 ist "Ist ZV aus [0,1] < p dann Ereigniss=true" nicht korrekt


Beispiel: Nehmen wir an RAND_MAX ist 9. Dann kann jeder Wert von 0-9 mit 10% Chance vorkommen. Wenn jetzt deine Wahrscheinlichkeit 0.01 ist würde bei keiner deiner Methoden jemals das Ereignis ausgelöst werden. Das ist mindestens genauso falsch wie wenn es in 10% der Fälle ausgelöst wird.

Das Ereigniss 0 liefert mit RandDouble() und RandDouble1() den Wert 0.0
und 0.0 < 0.01, also Ereigniss=true, allerdings mit starker Verzerrung

[hell_bird]

für p=1.0 ist "Ist ZV aus [0,1] < p dann Ereigniss=true" nicht korrekt
Ja schon, aber das sagtest du ja schon. Ich bin mir nicht ganz sicher worauf du hinaus willst. Mir scheinen deine Möglichkeit 2 und 3 beide richtig zu sein. Möglichkeit 3 ist die genauere. Die Probleme sind die Genauigkeit und Randfälle, aber anscheinend ist das nicht das was dich interessiert hat.

Das Ereigniss 0 liefert mit RandDouble() und RandDouble1() den Wert 0.0
und 0.0 < 0.01, also Ereigniss=true, allerdings mit starker Verzerrung
Richtig. Ich habe falsch gedacht. Korrekt wäre, dass alles über 90% wie 100% behandelt würde.

[Trap]

Alle 3 von dir in Post 1 genannten Versionen um ein double aus rand() zu bauen sind schlecht. Ein double im Bereich [0,1) kann 52 Bit Informationen speichern, die hat RAND_MAX nicht - garantiert sind nur 15 Bit.

Eine gute Lösung wäre http://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution oder auf [0,1) spezialisiert dann http://en.cppreference.com/w/cpp/numeric/random/generate_canonical

[pest]

ich habe doch schon geschrieben, dass ich überlicherweise den Mersenne-Twister verwende

die Vorgehensweise ist dann trotzdem analog

[Trap]

Warum rechnest du überhaupt mit Gleitkommazahlen?

[pest]

was meinst du? wie denn sonst?

[Trap]

Du kannst die p doch auch als Ganzzahl relativ zu RAND_MAX ausdrücken.

Also statt die Zufallszahlen umzurechnen die Grenzen umrechnen.

[pest]

klar, also irgendwie sowas


static bool BoolEvent(double p)
{
return (rand() < (p*(RAND_MAX+1))
}


wüsste nich was mir das bringt, da muss ich mit der rundung wieder aufpassen

[Trap]

Das hängt vom Rest deines Programms ab. Wenn du wirklich zwangsweise eine double->bool Funktion brauchst, bringt es fast nichts.

Eventuell macht es die Fehlersuche einfacher, da alle Fehler dann auf der nicht-zufälligen Seite des Vergleichs sind.

[pest]

naja, es geht ja Jetzt :), m.M. nach korrekt

die Basis-Funktionen brauche ich ja noch für andere Sachen
mit der [0,1]-Funktion bastle ich mir zufällige Intervalle
und mit der (0,1)-Funktion erstelle Normalverteilte ZV
Der BoolEvent richtig sich auch nach einer stetigen Verteilung p(i)

Das ist so schon rel. praktisch

[Milchkanne]

Wo Trap das jetzt sagt, stell ich mir gerade die Frage: Jetzt sind die Zahlen als Double ja alles andere als gleichverteilt. Berücksichtigen das die Zufallsgeneratoren überhaupt? Nicht dass die Ungleichung die Randfälle dann korrekt berücksichtigt, aber die Zahlen überhaupt nicht gleichverteilt sind...

[pest]

doch deswegen z.B.

return (double(rand())+0.5)/(double(RAND_MAX)+1.0);


statt

return (double(rand())+1.0)/(double(RAND_MAX)+2.0);

[pest]

nochmal hierzu

Möglichkeit 3 ist die genauere.

warum?


Die Probleme sind die Genauigkeit und Randfälle, aber anscheinend ist das nicht das was dich interessiert hat.


deswegen mach ich ja den ganzen Spass hier ;(
ob nun rand() der Mersenne-Twister mit RAND_MAX=0xffffffff ist, spielt doch erst einmal keine Rolle

[Coda]

ach ist doch overkill, ne verteilungsfunktion kann ich mir selber basteln
Ja, kannst du ist aber NIH-Syndrom.

Und bei so Stochastik-Zeug macht man sehr leicht kleine Fehler, deshalb würde ich den getesteten Code immer bevorzugen, selbst wenn ich es wohl selber machen könnte.

Deine Entscheidung :)

[pest]

Und bei so Stochastik-Zeug macht man sehr leicht kleine Fehler

das stimmt, deswegen nehme ich es genau, auch wenn es für etwaige Simulationsergebnisse im Mittel wenig Relevanz besitzt


deshalb würde ich den getesteten Code immer bevorzugen


naja ich habe jetzt in Boost keinen Code für [0,1] und (0,1)-Intervall gefunden - Ersteres bräuchte ich zwingend

[Coda]

http://www.boost.org/doc/libs/1_51_0/doc/html/boost/random/uniform_01.html
http://www.boost.org/doc/libs/1_51_0/doc/html/boost/random/uniform_real_distribution.html

Wenn du wirklich das geschlossene Interval brauchst, dann benutzt

boost::math::nextafter(1.0, std::numeric_limits<double>::max())

für das Ende.