Ich ging so vor: Wenn nach Schritt n die Wahrscheinlichkeit, dass k Kekse von einer bestimmten Sorte vorhanden sind, p_n(k) beträgt, dann kann man daraus auch die entsprechende Verteilung für den nächsten Schritt berechnen, indem man 1001 verschiedene Binomialverteilungen mit jeweils p=k/1000 aufsummiert, nach p_n(k) gewichtet.
Also p_(n+1)(j) = sum_(k=0)^1000 p_n(k) * (1000 choose j) (k/1000)^j (1-k/1000)^(1000-j)
Da alles außer p_n(k) hier nicht von n abhängt, kann man das auch so betrachten, dass man die ursprüngliche Verteilung (eine Binomialverteilung mit p=0.5) 479-mal mit der Matrix aus den restlichen Termen (abhängig von k und j) multipliziert.
Ich habe mir aber natürlich nicht die Mühe gemacht, eine 1001x1001-Matrix zu diagonalisieren, um sie potenzieren zu können, sondern einfach alle Iterationsschritte mit dem Computer durchgeführt. Hier sieht man auch schön, wie sich die Verteilung immer mehr auf die beiden Extremwerte fokussiert. Aus der endgültigen Verteilung kann man dann (ebenfalls mit dem Computer) die Wahrscheinlichkeit von etwa 0.7 für zwei gleiche Kekse ausrechnen.
Gut zu wissen, dass es mit dem richtigen Modell aber auch ohne Computer oder riesige Matrizen geht. Ich ging auch davon aus, dass das der Fall sei, habe aber einen entsprechenden Lösungsweg nicht gefunden und freue mich, diesen hier zu sehen.
