guter zufall (multiply-with-carry prng)

erik.vikinger

Hallo,

Kurzfassung: es gibt Leute die wissen was sie tun,

Richtig, die können selber prüfen was sie wollen. Wir beide sind bestimmt in der Lage einen ?RNG grob auf seine Tauglichkeit zu prüfen.

und es gibt Leute die haben keinen Tau.

Und gerade um die geht es mir. Warum gibt es die Kennzeichnungspflichtig für Eier im Einzelhandel? Damit auch jemand der keine Ahnung vom Fach hat in der Lage ist eine vernünftige Entscheidung zu treffen. Und genau das wird verhindert wenn auf jedem X-beliebigen Algorithmus "Zufall" drauf steht obwohl gar kein Zufall drin ist. Der durchschnittliche Programmierer, der vom seltenen Spezialgebiet der "Nutzbarmachung des Zufalls" keine Ahnung hat, wird sich nicht die Mühe machen und prüfen ob Library X seinen Ansprüchen wirklich genügt sondern sich auf das Wort "Zufall" verlassen. Genau deshalb bin ich dafür das diese PRNGs das R im Nahmen nicht haben dürften. Und nur weil jeder auf seinen Kram etwas draufschreibt was gar nicht drin ist wird das deshalb nicht weniger Falsch.

Wenn jemand weiss was er tut, wird er die Warnung nicht brauchen.

Möglicherweise, aber sie könnte ihm die Mühe einer eigenen Prüfung sparen.

Wenn er nicht weiss was er tut, wird er damit nichts anfangen können.

Da möchte ich nur auf die Eier verweisen, das System funktioniert IMHO recht gut (inwiefern sich auch alle Erzeuger korrekt dran halten kann ich natürlich nicht beurteilen).

Grüße
Erik

erik.vikinger schrieb:

Genau deshalb bin ich dafür das diese PRNGs das R im Nahmen nicht haben dürften. Und nur weil jeder auf seinen Kram etwas draufschreibt was gar nicht drin ist wird das deshalb nicht weniger Falsch.

Dafür ist doch das P da. Mal ne andere Frage wie bekomme ich von einem "Zufalls"generator x "Zufalls"zahlen im Bereich [10 100)?

Basher_ist_nicht_BashAr

Tim schrieb:

Mal ne andere Frage wie bekomme ich von einem "Zufalls"generator x "Zufalls"zahlen im Bereich [10 100)?

meinste das: http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html ?

hustbaer

Ich gehe mal davon aus dass der Output des Generators gleichverteilt ist, und dass man das Ergebnis auch wieder gleichverteilt haben will.

----

Im ersten Schritt bringt man normalerweise den Anfang der Range auf 0, d.h. wir suchen eine Zahl im Bereich [0, 90). (Auf die addieren wir nachher einfach 10 drauf, damit wir wieder auf [10, 100) kommen.)

Nun könnte man einfach eine Zahl ziehen, und sie mittels Modulo in den gewünschten Bereich bringen. Bloss würde das die Gleichverteilung zerstören.

Angenommen der Generator hat eine Output-Range von [0, 200).
Dann würden 0, 90 und 180 mittels Module auf 0 abgebildet,
1, 91 und 181 auf 1,
...
19, 109 und 199 auf 19,
20 und 110 auf 20,
21 und 111 auf 21,
etc.

Kurz: die Zahlen 0 bis 19 würden 1,5 mal häufiger vorkommen als die Zahlen 20 bis 89.

Die IMO einfachste (und AFAIK "üblichste") Möglichkeit dem zu begegnen, ist, gezogene Zahlen von 180 bis 199 einfach zu verwerfen.

Verallgemeinert, in Pseudo-Code, könnte das ca. so aussehen:

int GetNumber(Generator& g, int lowerBound, int upperBoundExclusive)
{
    int newUpperBoundExclusive = upperBoundExclusive - lowerBound;
    int generatorUpperBoundExclusive = g.GetUpperBoundExclusive(); // wir gehen davon aus dass die untere Schwelle des Generators immer 0 ist
    int n = generatorUpperBoundExclusive / newUpperBoundExclusive; // integer Division, Ergebnis immer abgerundet
    assert(n > 0); // sonst haben wir ein Problem
    int acceptableUpperBoundExclusive = n * newUpperBoundExclusive;

    int x = 0;
    do
    {
        x = g.GetNumber();
    }
    while (x >= acceptableUpperBoundExclusive);

    x = x % newUpperBoundExclusive;
    x = x + lowerBound;
    return x;
}

Das aufwendigste an der ganzen Sache ist dabei, sicherzustellen, dass nirgends ein Überlauf passiert. Und, falls die geforderte Range grösser werden kann, als die die der Generator ausspuckt, muss man sich natürlich statt dem einfachen assert() etwas besseres einfallen lassen. Beispielsweise in jedem Durchlauf zwei (oder mehr) Zahlen zu ziehen, und diese zu einer grösseren zusammenzusetzen.

Blue-Tiger

Basher schrieb:

Tim schrieb:

Mal ne andere Frage wie bekomme ich von einem "Zufalls"generator x "Zufalls"zahlen im Bereich [10 100)?

meinste das: http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html ?

Oder sowas: http://www.mathematik.uni-ulm.de/stochastik/lehre/ss03/markov/skript/node28.html ?

SeppJ

@hustbaer: In diesem Fall eine ungümstige Methode, da du unnötigerweise Zahlen wegwirfst. Blue-Tigers Link zeigt auf ein paar bessere, sehr allgemeine Methoden, die jedoch für den Fall einer Gleichverteilung [10,100) wie Kanonen auf Spatzen sind. Wenn man (wie im Computer üblich) von einer gleichmäßigen Verteilung von Zufallszahlen X im Bereich [0,RMAX) ausgeht, bekommt man diese direkt mittels (X/RMAX)*90+10 auf den fraglichen Bereich transformiert.

hustbaer

@SeppJ:
Nö, das kannst du knicken. Üblicherweise werfen Zufallszahlengeneratoren Ganzzahlen aus. Und wenn du deine Transformation auf Ganzzahlen loslässt (und das Ergebnis wieder auf Ganzzahlen "beschneidest"), dann zerstörst du damit die Gleichverteilung.
(Mal davon abgesehen dass bei Integer-Arithmetik die Division nach der Multiplikation kommen muss, damit überhaupt irgendwas annähernd sinnvolles rauskommen kann.)

Das Ergebnis ist gleich schlecht, wie wenn du Modulo (ohne Verwerfen) als Transformation verwendest. Der einzige Unterschied ist, dass die "bevorzugten" Zahlen sich nicht wie bei Modulo alle am Anfang der Range sammeln, sondern schön über die ganze Range verteilt werden. Bringt mir aber auch nicht viel, wenn ich will, dass alle Zahlen in der Output-Range genau die gleiche Wahrscheinlichkeit haben.

----

Ich persönlich glaube dass das Verwerfen von Zahlen die einfachste und daher beste Methode ist.

Natürlich gäbe es auch andere Möglichkeiten die die Gleichverteilung nicht zerstören. z.B. könnte man einen "Range-Coder" (den Decoder-Teil genauergenommen) verwenden, und damit einen Stream aus zufälligen Bits dekodieren. (Die Wahrscheinlichkeiten/Häufigkeiten würde man dazu alle gleich annehmen, um die gewünschte Gleichverteilung zu erreichen.)

http://de.wikipedia.org/wiki/Arithmetische_Kodierung

Das würde ich dann aber wirklich mit Kanonen auf Spatzen schiessen nennen.

----

BTW: die Implementierung in der Boost macht genau das was ich vorhin skizziert habe: Zahlen verwerfen.

Basher_ist_nicht_BashAr

hustbaer schrieb:

Nö, das kannst du knicken. Üblicherweise werfen Zufallszahlengeneratoren Ganzzahlen aus.

hihi, vielleicht war [a b**)** nur ein tippfehler. damit war weder das mathematische 'halboffene intervall' noch reelle zahlen gemeint (das komma fehlt ja auch).

SeppJ

hustbaer schrieb:

@SeppJ:
Nö, das kannst du knicken. Üblicherweise werfen Zufallszahlengeneratoren Ganzzahlen aus. Und wenn du deine Transformation auf Ganzzahlen loslässt (und das Ergebnis wieder auf Ganzzahlen "beschneidest"), dann zerstörst du damit die Gleichverteilung.

Das sehe ich überhaupt nicht ein, kannst du das mal erläutern?

X aus [0,RMAX) gleichverteilt, meinetwegen nur Ganzzahlen

=> X/RMAX in [0,1) gleichverteilt (und man MUSS vorher dividieren, da normalerweise RMAX auch gleichzeitig die größte darstellbare Zahl ist). Ich nehme mal an, jeder im Forum weiß, wie man mit einem Computer irchtig dividiert.

=> (X/RMAX) * 90 in [0,90) gleichverteilt, jedoch als Fließkommazahlen

=> (X/RMAX) * 90 + 10 in [10,100) gleichverteilt, jedoch als Fließkommazahlen

Jetzt kann man natürlich nur den ganzzahligen Teil der Zahlen nehmen. Da die Werte aber gleichverteilt waren, wird dabei kein bestimmter Wert bevorzugt.

volkard

SeppJ schrieb:

Jetzt kann man natürlich nur den ganzzahligen Teil der Zahlen nehmen. Da die Werte aber gleichverteilt waren, wird dabei kein bestimmter Wert bevorzugt.

Falsch.
Wenn ich 100 Gäste habe und nur 90 Zimmer, MÜSSEN gelegentliche Doppelbelegungen vorkommen und die Gleichverteilung ist kaputt.
Aber auch wenn ich 1000000 Betten habe und 90 Zimmer, kann es keine Gleichverteilung geben, weil 90 einfach kein Teiler von 1000000 ist.
Wenn ich aber nur 10 Gäste ablehne (mit einem Einnahmeverlust von (10/1000000) einem Tausendstel Prozent, kann ich die verbleibenden 999990 Gäste prima auf 90 Zimmer gleichverteilen.

volkard

SeppJ schrieb:

Das sehe ich überhaupt nicht ein, kannst du das mal erläutern?

Das Skriptum war ein wenig Weltfremd. Das ist normal.

X aus [0,RMAX) gleichverteilt, meinetwegen nur Ganzzahlen

OK. Also wie rand().

=> X/RMAX in [0,1) gleichverteilt

Nö. Nur RMAX ausgesuchte Werte kommen vor und nicht alle reellen Zahlen zwischen 0 und 1. Das ist ein riesiger Unterschied.

=> (X/RMAX) * 90 in [0,90) gleichverteilt, jedoch als Fließkommazahlen

Nicht haltbar, da die Vorbedingung schon falsch war.

=> (X/RMAX) * 90 + 10 in [10,100) gleichverteilt, jedoch als Fließkommazahlen

Nicht haltbar, da die Vorbedingung schon falsch war.

Jetzt kann man natürlich nur den ganzzahligen Teil der Zahlen nehmen. Da die Werte aber gleichverteilt waren, wird dabei kein bestimmter Wert bevorzugt.

Nicht haltbar, da die Vorbedingung schon falsch war.

SeppJ

Ah, jetzt verstehe ich, was ihr gemeint habt. Ja, da habt ihr wohl Recht. Ich habe nicht bedacht, dass ich mit der Transformation nur (quasi-)kontinuierliche Verteilungen erzeugen kann, was für Ganzzahlen zwischen 10 und 100 nicht gegeben ist.

Aber es muss doch was besseres geben als das was hustbaer vorgeschlagen hat. Wenn ich im ungünstigsten Fall Ganzzahlen zwischen 0 und RMAX/2+1 suche, müsste ich ja fast die Hälfte der Zahlen wegwerfen. Mir fällt spontan jedoch auch nichts besseres ein.

life

Wenn Java das auch so macht, wirds schon richtig sein ;):

public int nextInt(int n)

    Returns a pseudorandom, uniformly distributed int value between 0 (inclusive) and the specified value (exclusive), drawn from this random number generator's sequence. The general contract of nextInt is that one int value in the specified range is pseudorandomly generated and returned. All n possible int values are produced with (approximately) equal probability. The method nextInt(int n) is implemented by class Random as follows:

         public int nextInt(int n) {
             if (n<=0)
        		throw new IllegalArgumentException("n must be positive");

             if ((n & -n) == n)  // i.e., n is a power of 2
                 return (int)((n * (long)next(31)) >> 31);

             int bits, val;
             do {
                 bits = next(31);
                 val = bits % n;
             } while(bits - val + (n-1) < 0);
             return val;
         }

    The hedge "approximately" is used in the foregoing description only because the next method is only approximately an unbiased source of independently chosen bits. If it were a perfect source of randomly chosen bits, then the algorithm shown would choose int values from the stated range with perfect uniformity.

    The algorithm is slightly tricky. It rejects values that would result in an uneven distribution (due to the fact that 2^31 is not divisible by n). The probability of a value being rejected depends on n. The worst case is n=2^30+1, for which the probability of a reject is 1/2, and the expected number of iterations before the loop terminates is 2.

    The algorithm treats the case where n is a power of two specially: it returns the correct number of high-order bits from the underlying pseudo-random number generator. In the absence of special treatment, the correct number of low-order bits would be returned. Linear congruential pseudo-random number generators such as the one implemented by this class are known to have short periods in the sequence of values of their low-order bits. Thus, this special case greatly increases the length of the sequence of values returned by successive calls to this method if n is a small power of two.

    Parameters:
        n - the bound on the random number to be returned. Must be positive. 
    Returns:
        a pseudorandom, uniformly distributed int value between 0 (inclusive) and n (exclusive). 
    Throws:
        IllegalArgumentException - n is not positive.
    Since:
        1.2

volkard

SeppJ schrieb:

Aber es muss doch was besseres geben als das was hustbaer vorgeschlagen hat. Wenn ich im ungünstigsten Fall Ganzzahlen zwischen 0 und RMAX/2+1 suche, müsste ich ja fast die Hälfte der Zahlen wegwerfen. Mir fällt spontan jedoch auch nichts besseres ein.

Du sagst "Nur die Hälfte? :("
Sag mal "Nur die Hälfte! :)"

Außerdem braucht man das selten. Ungefähr immer geht rand()%49, außer ganz manchmal, dann nimmt man rand()/double(RAND_MAX+1)*49, und nur wenn das auch nicht geht, nimmt man die 1337-Version. Normalerweise stört an der Modulo-Version nämlich gar nicht die Ungleichverteilung, sondern die Schieflage, weil sich die volleren Zimmer links sammeln und dort Slums entstehen, da bringt die double-Version schon Abhilfe.

volkard

//Multiply-With-Carry Pseudo Random Number Generator
class Random {
private:
	static UInt64 const a=4294967118;
	UInt64 x;
public:
	Random(UInt64 seed=0) {
		x=seed+!seed;
	}
	UInt32 operator()() {
		x=a*(x&0xffffffff)+(x>>32);
		return x;
	}
};

//Numerical Receipes sagt:
// Wenn der Seed in 1 <= x <= 2^32-1 liegt 
// und sowohl a*b-1 als auch (a*b-2)/2 Primzahlen sind, 
// dann ist die Periodenlänge gleich (a*b-2)/2. 
//Volkard sagt: 
// Die größte Zahl, die diese Bedingungen erfüllt, ist 
// obiges a=2^32-178 mit beinahe 2^63 als Periodenlänge.

Basher_ist_nicht_BashAr

nimm doch den: http://burtleburtle.net/bob/rand/isaacafa.html