Python, Numpy and Wahrscheinlichkeiten

Einführung

dice

"Every American should have above average income, and my Administration is going to see they get it."

Wir beginnen mit diesem Zitat von Bill Clinton unser Kapitel über Statistik, weil es ein besonders hervorstechendes Beispiel eines weit vorherrschenden Problems ist. Ein Problem, was der Mathematiker John Allen Paulos in einem seiner Bücher als "Innumeracy" bezeichnet hatte. (Paulos, John Allen. 1988. "Innumeracy, Mathematical Illiteracy and its consequences". New York. Hill and Wang) Ein Ausdruck für den es noch keine deutsche Übersetzung gibt, aber soviel bedeutet wie "numerischer Analphabetismus". Darunter versteht er die weit-verbreitete Unfähigkeit einfachste statistische und wahrscheinlichkeitstheoretische Zusammenhänge zu begreifen. Während Analphabetismus in allen Gesellschaften und Gesellschaftskreisen geächtet wird, gilt "innumeracy" als tolerierbar. Nehmen wir an, es gäbe ein Land, in dem fast alle Leute eine grüne Hautfarbe hätten, außer einer blaufarbigen Minderheit. Nun können wir davon ausgehen, dass es unter beiden Gruppen Verbrecher geben wird. Stellen wir uns aber vor, dass innerhalb kurzer Zeit Zeitungen mit Artikeln der Form "Blauer mordet Mutter von zwei Kindern" oder "Polizei ermittelt erneut gegen Blauen im Fall XYZZ!" veröffentlicht werden. Über die gleichzeitig stattfindenden Straftaten von Menschen grüner Hautfarbe wird jedoch nicht berichtet. Numerischer Analphabetismus bedeutet nun, dass Leute auf die Strasse gehen und gegen die durch "Blaue" ausgelöste Gewalt demonstrieren, Politiker ein mögliches Problem bei Blauen untersuchen wollen oder gar zu erkenen glauben und so weiter.


Statistik und Wahrscheinlichkeites-Berechnung begegnet uns überall in unserem Leben. Wir müssen damit umgehen, immer wenn wir zwischen verschiedenen Alternativen zu wählen haben. Können wir morgen wandern gehen oder wird es regnen? Die Wettervorhersage sagt uns, dass die Niederschlags-Wahrscheinlichkeit bei 30% liegt. Was jetzt? Können wir es riskieren?
Andere Senario: Sie sind in den Ferien auf einer paradiesischen Insel, weit weg von zu Hause. Plötzlich treffen Sie ihren Nachbarn! Wie hoch sind die Chancen, dass so etwas passiert?
Sie spielen jede Woche Lotto und träumen von der weit entfernten Insel. Wie hoch ist die Wahrscheinlichkeit den Jackpot zu gewinnen, so dass Sie niemals mehr arbeiten müssen und im "Paradies" leben können?

Die Ungewissheit umgibt uns und nur einige Menschen verstehen die Grundlagen der Wahrscheinlichkeits-Theorie.

Die Programmiersprache Python und die Module Numpy und Scipy helfen uns nicht, die oben genannten alltäglichen Probleme zu lösen, jedoch bieten Python und Numpy starke Funktionalitäten, um Probleme der Statistik und Wahrscheinlichkeits-Theorie zu berechnen.



Zufallszahlen mit Python

Die Module random und secrets

Das Modul secrets wurde erst mit Python 3.6 neu eingeführt. Mit diesem Modul kann man kryptografisch starke Pseudozufallszahlen erzeugen, die sich als Passwörter, Tokens oder ähnliches eignen. Dieses Modul wurde mit diesem Ziel entwickelt. Als Generator enthält es eine CSPRNG (Cryptographically Strong Pseudo Random Number Generator). Das random-Modul von Python wurde nicht in Hinblick auf kryptografische Anwendungen entwickel. Der Fokus dieses Moduls lag auf Modellbildungen und Simulationen.

Es gibt sogar eine ausdrückliche Warnung in der Dokumentation des random-Modules:



Warnung: Beachten Sie bitte, dass der Pseudo-Zufalls-Generator im random-Modul NICHT für sicherheitsrelevante Zwecke benutzt werden sollte. Benutzen Sie secrets ab Python 3.6+ und os.urandom() bei Python 3.6+ und früher.
Im Original: "Note that the pseudo-random generators in the random module should NOT be used for security purposes. Use secrets on Python 3.6+ and os.urandom() on Python 3.5 and earlier."

Man kann allerdings die SystemRandom-Klasse des random-Moduls verwenden, die die sichere os.urandom-Systemfunktion benutzt. random.SystemRandom und secrets.SystemRandom sind identisch.

Mit der Funktion random.random können wir ein Zufallszahl im halb-offenen Interval [0, 1) erzeugen:

import random
random_number = random.random()
print(random_number)
0.4197271123707619

Nun wollen wir eine kryptografisch starke Zufallszahl erzeugen:

from secrets import SystemRandom
# from random import SystemRandom   # äquivalent
crypto = SystemRandom()
print(crypto.random())
0.8687421606095797

Erzeugen einer Liste von Zufallszahlen

Häufig benötigen wir mehr als eine Zufallszahl. Die folgende Funktion random_list kann eine Liste von Zufallszahlen mit vorgegebener Anzahl erzeugen. Mit dem dem Parameter secure können wir wir steuern, ob wir sichere, also mit SystemRandom erzeugte Zahlen wollen:

import random
def random_list(n, secure=True):
    random_floats = []
    if secure:
        crypto = random.SystemRandom()
        random_float = crypto.random
    else:
        random_float = random.random
    for _ in range(n):
        random_floats.append(random_float())
    return random_floats
print(random_list(10, secure=False))
[0.9205923800470898, 0.9656150728749113, 0.7372850291299475, 0.7318383271049008, 0.1963972136060309, 0.3210564104198963, 0.8562302503940156, 0.9395264233003696, 0.10181184275402566, 0.12427810857739963]

Im Folgenden können wir sehen, dass der Preis für die Sicherheit eine deutliche Erhöhung der Rechenzeit bedeutet:

%%timeit
random_list(100, secure=True)
147 µs ± 306 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%%timeit
random_list(100, secure=False)
6.44 µs ± 72.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Am einfachsten, was den Programmieraufwand betrifft, und am schnellsten, was die Laufzeit betrifft, lassen sich Zufallszahlen mit dem random-Untermodul von numpy erzeugen:

import numpy as np
np.random.random(10)
Führt man obigen Code aus, erhält man Folgendes:
array([0.08839793, 0.4801768 , 0.16206772, 0.2220249 , 0.7012845 ,
       0.53328173, 0.95577474, 0.15833697, 0.06686501, 0.18246386])
%%timeit
np.random.random(100)
1.47 µs ± 51.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

Allerdings werden auch hier keine sicheren Zufallszahlen erzeugt!



Zufällige Integerzahlen mit Python

Jeder ist mit der Generierung von Zufallszahlen ohne Computer vertraut. Wenn Sie einen Würfel werfen, erhalten Sie eine Zufallszahl zwischen 1 und 6. Im Zusammenhang mit der Wahrscheinlichkeits-Theorie nennen wir "Werfen eines Würfels" ein Experiment, dessen Menge der möglichen Ergebnisse {1, 2, 3, 4, 5, 6} ist. Diese Menge wird auch als "Stichprobenraum" bezeichnet.

Integer-Zahlen könnten wir relativ einfach aus dem erzeugen, was wir schon kennengelernt haben. Unsere Funktion erzeugt Zahlen im halboffenen Intervall zwischen start und stop:

def randint(start, stop, secure=True):
    if secure:
        crypto = random.SystemRandom()
        rnd = crypto.random
    else:
        rnd = random.random
    result = int(rnd() * (stop-start)) + start
    return result
[ randint(3, 7) for i in range(20)]
Wir können die folgende Ausgabe erwarten, wenn wir den obigen Python-Code ausführen:
[3, 4, 6, 6, 5, 5, 4, 3, 3, 6, 3, 6, 4, 6, 3, 3, 6, 6, 6, 6]

Es empfiehlt sich jedoch die Funktionen randint zu nutzen, die die Module random, secrets und numpy.random bieten:

import random
# randint
outcome = random.randint(1, 6)
print(outcome)
2
import secrets
c = secrets.SystemRandom()
c.randint(1, 6)
Wir erhalten die folgende Ausgabe:
4

Würfeln wir sehr oft, so sollte jede Augenzahl ungefähr in einem sechstel der Fällte auftreten. Wir testen dies im Folgenden. Das Zählen erledigen wir mit der Counter-Klasse aus dem collections-Modul:

from collections import Counter
import random
c = Counter()
anzahl_wuerfe = 100000
for wurf in range(anzahl_wuerfe):
    augenzahl = random.randint(1, 6)
    c[augenzahl] += 1
    
print(c)
for augenzahl in c:
    print(f"Augenzahl: {augenzahl}, relative Häufigkeit: {c[augenzahl]/anzahl_wuerfe}")
Counter({1: 16850, 2: 16760, 3: 16647, 5: 16638, 6: 16627, 4: 16478})
Augenzahl: 2, relative Häufigkeit: 0.1676
Augenzahl: 6, relative Häufigkeit: 0.16627
Augenzahl: 5, relative Häufigkeit: 0.16638
Augenzahl: 4, relative Häufigkeit: 0.16478
Augenzahl: 1, relative Häufigkeit: 0.1685
Augenzahl: 3, relative Häufigkeit: 0.16647

Obiges Beispiel lässt sich mit der random-Funktion aus dem Numpy-Modul einfacher realisieren.

import numpy as np
from collections import Counter
anzahl_wuerfe = 100000
outcome = np.random.randint(1, 7, size=anzahl_wuerfe)
c = Counter(outcome)
for augenzahl in c:
    print(f"Augenzahl: {augenzahl}, relative Häufigkeit: {c[augenzahl]/anzahl_wuerfe}")
Augenzahl: 4, relative Häufigkeit: 0.16552
Augenzahl: 5, relative Häufigkeit: 0.16602
Augenzahl: 1, relative Häufigkeit: 0.167
Augenzahl: 6, relative Häufigkeit: 0.1676
Augenzahl: 2, relative Häufigkeit: 0.16735
Augenzahl: 3, relative Häufigkeit: 0.16651

Sicherlich haben Sie bemerkt, dass wir in random.randint beim zweiten Parameter ein 7 statt einer 6 verwendet haben. Diese Funktion nutzt ein "halb-offenes" Intervall im Gegensatz zu Pythons random-Modul, dass ein geschlossenes Intervall erwartet.

Die formale Definition:

numpy.random.randint(low, high=None, size=None)

Diese Funktion liefert zufällige Integers zwischen "low" (inklusiv) und "high" (exklusiv). Mit anderen Worten: randint liefert zufällige Integers aus der "Diskreten Uniform" Distribution im "halb-offenen" Intervall ["low", "high").

Wenn "high" None oder nicht übergeben wird, liegen die Ergebnisse in der Range von [0, "low").

Der Parameter "size" definiert die Shape des Ergebnisses. Wenn "size" None ist, wird nur ein Integer-Wert generiert. Andernfalls ist das Ergebnis ein Array. "size" sollte ein Tupel sein. Wenn als "size" ein Integer n übergeben wird, entspricht dies dem Tupel (n,).

Folgende Beispiele verdeutlichen das Verhalten der Parameter:

import numpy as np
print(np.random.randint(1, 7))
print(np.random.randint(1, 7, size=1))
print(np.random.randint(1, 7, size=10))
print(np.random.randint(1, 7, size=(10,))) # the same as the previous one
print(np.random.randint(1, 7, size=(5, 4)))
3
[3]
[5 4 1 4 5 6 1 6 6 1]
[4 5 3 5 4 4 1 1 6 6]
[[2 3 5 1]
 [5 2 5 1]
 [2 4 6 5]
 [3 2 3 4]
 [6 2 4 4]]

Wir können das Würfel-Werfen mit dem Code np.random.randint(1, 7, size = 1) aus dem Numpy-Modul simulieren oder anhand des Codes random.randint(1, 6) aus dem Standard random-Modul. Wir nehmen an, dass unser Würfel fair ist, d.h. die Wahrscheinlichkeit für jede Seite liegt bei 1/6. Wie können wir einen krummen oder belasteten Würfel simulieren? Die randint-Methoden beider Module sind dafür nicht geeignet. Wir werden im folgenden ein paar Funktionen schreiben um das Problem zu lösen.

Zunächst schauen wir uns weitere nützliche Funktionen des random-Moduls an.



Zufällige Auswahlen (Choices) und Samples in Python


"choice" ist eine weitere extrem nützliche Funktion des random-Moduls. Sie kann genutzt werden um aus einer "nicht leeren" Sequenz ein zufälliges Element zu wählen.

Das bedeutet, wir sind dazu in der Lage aus einem String ein zufälliges Zeichen zu holen oder ein zufälliges Element aus einer Liste oder eins Tupels. Wie Sie im folgenden Beispiel sehen können:

In [ ]:
from random import choice
professions = ["scientist", "philosopher", "engineer", "priest"]
print(choice("abcdefghij"))
print(choice(professions))
print(choice(("apples", "bananas", "cherries")))

"choice" liefert ein Objekt aus einer Sequenz. Die Chancen für die Elemente ausgewählt zu werden, sind gleichmäßig verteilt. Die Chance, dass "scientist" zurück geliefert wird beim Aufruf von choice(profession) liegt bei 1/4. Das hat allerdings wenig mit der Realität zu tun. Es gibt sicherlich mehr Wissenschaftler und Ingenieure auf der Welt als Priester und Philosophen. Wie bei dem belasteten Würfel brauchen wir auch hier eine gewichtete Auswahl.

Wir definieren eine Funktion weighted_choice die ein zufälliges Element einer Sequenz zurückliefert wie random.choice, jedoch sind die Element der Sequenz gewichtet.


Zufalls-Intervalle

Bevor wir mit dem Design der gewichteten Auswahl beginnen, definieren wir eine Funktion find_interval(x, partition), die wir für unsere weighted_choice-Funktion brauchen. find_interval erwartet zwei Argumente:

Die Funktion liefert i zurück, wenn pi < x < pi+1. -1 wird zurückgeliefert, wenn x kleiner als p0 ist oder x größer oder gleich pn ist.

partitions In [ ]:
def find_interval(x, partition):
    """ find_interval -> i
        partition is a sequence of numerical values
        x is a numerical value
        The return value "i" will be the index for which applies
        partition[i] < x < partition[i+1], if such an index exists.
        -1 otherwise
    """
    
    for i in range(0, len(partition)):
        if x < partition[i]:
            return i-1
    return -1
I = [0, 3, 5, 7.8, 9, 12, 13.8, 16]
for x in [-1.3, 0, 0.1, 3.2, 5, 6.2, 7.9, 10.8, 13.9, 15, 16, 16.5]:
    print(find_interval(x, I), end=", ")



Gewichtete Zufalls-Auswahl

Jetzt können wir die weighted_choice-Funktion definieren. Nehmen wir an, dass wir 3 Gewichtungen haben, d.h. 1/5, 1/2 und 3/10. Wir bilden die kumulative Summe der Gewichtungen mit np.cumsum(weights).

In [ ]:
import numpy as np
weights = [0.2, 0.5, 0.3]
cum_weights = [0] + list(np.cumsum(weights))
print(cum_weights)

Wenn wir eine Zufallszahl x zwischen 0 und 1 mit random.random() generieren, so ist die Wahrscheinlichkeit, dass x im Intervall [0, cum_weights[0]) liegt, 1/5. Die Wahrscheinlichkeit, dass x im Intervall [cum_weights[0], cum_weights[1]) liegt, ist 1/2. Letztlich ist die Wahrscheinlichkeit, dass x im Intervall [cum_weights[1], cum_weights[2]) , bei 3/10.

Jetzt sind Sie in der Lage die grundlegende Idee zu verstehen, wie weighted_choice funktioniert:

In [ ]:
import numpy as np
def weighted_choice(sequence, weights):
    """ 
    weighted_choice selects a random element of 
    the sequence according to the list of weights
    """
    x = np.random.random()
    cum_weights = [0] + list(np.cumsum(weights))
    index = find_interval(x, cum_weights)
    return sequence[index]

Beispiel:

Wir können die Funktion weighted_choice nun für die folgende Aufgabe verwenden: Vermuten wir, dass wir einen "belasteten" Würfel haben mit P(6)=1/12 und P(1)=1/12. Die Wahrscheinlichkeit für alle anderen möglichen Ergebnisse sind gleich, d.h. P(2) = P(3) = P(4) = P(5) = p. Wir können p wie folgt ausrechnen mit

1 - P(1) - P(6) = 4 x p.

Das entspricht

p = 1/6.

Wie können wir diesen Würfel mit unserer weighted_choice-Funktion simulieren.

Wir rufen weighted_choice auf mit "Seiten des Würfels"- und der "Gewichtungs"-Liste auf. Jeder Aufruf entspricht einem Wurf des belasteten Würfels.

Wir sehen, dass nach 10.000 Würfen die geschätzte Wahrscheinlichkeit der Gewichtung entspricht:

In [ ]:
from collections import Counter
faces_of_die = [1, 2, 3, 4, 5, 6]
weights = [1/12, 1/6, 1/6, 1/6, 1/6, 3/12]
outcomes = []
n = 10000
for _ in range(n):
    outcomes.append(weighted_choice(faces_of_die, weights))
c = Counter(outcomes)
for key in c:
    c[key] = c[key] / n
    
print(sorted(c.values()))

Die Werte der Aufteilungs-Liste definieren die Unter-Aufteilungen, in denen wir den Wert x erwarten. Wenn der Wert x kleiner als p0 oder grösser/gleich pn, liefern wir -1 zurück.

Wir könnten unsere erste Unter-Aufteilung als ein Intervall definieren von -unendlich bis p0 und 0 zurückliefern. Die letzte Unter-Aufteilung wäre dann ein Intervall von pn bis unendlich.

Um zwischen beiden Fällen unterscheiden zu können, erweitern wir die Funktion find_interval um den Parameter "endpoints". "True" entspricht unserem Eingangs erwähnten Vorgehen. "False" entspricht unserem soeben beschriebenen Fall.

In anderen Worten - Wenn "endpoints" False ist, passiert folgendes:

Wir demonstrieren dies im folgenden Diagramm:

partitions without endpoints

Die neue Funktion sieht folgendermaßen aus:

In [ ]:
def find_interval(x, 
                  partition, 
                  endpoints=True):
    """ 
    find_interval -> i
    If endpoints is True, "i" will be the index for which 
    applies partition[i] < x < partition[i+1], if such an
    index exists. -1 otherwise
        
    If endpoints is False, "i" will be the smallest index 
    for which applies x < partition[i]. If no such index 
    exists 
    "i" will be set to len(partition)
    """
    for i in range(0, len(partition)):
        if x < partition[i]:
            return i-1 if endpoints else i
    return -1 if endpoints else len(partition)
I = [0, 3, 5, 7.8, 9, 12, 13.8, 16]
print("Endpoints are included:")
for x in [-1.3, 0, 0.1, 3.2, 5, 6.2, 7.9, 13.9, 15, 16, 16.5]:
    print(find_interval(x, I), end=", ")
print("\nEndpoints are not included:")
for x in [-1.3, 0, 0.1, 3.2, 5, 6.2, 7.9, 13.9, 15, 16, 16.5]:
    print(find_interval(x, I, endpoints=False), end=", ")



Zufalls-Samples mit Python


Eine Sample kann als ein repräsentativer Teil einer größeren Gruppe angesehen werden, die wir "Population" nennen.

Das Modul numpy.random beinhaltet die Funktion random_sample, die zufällige Float-Werte liefert im halb-offenen Intervall [0.0, 1.0). Die Ergebnisse sind gleichmäßig verteilt über das angegebene Intervall. Die Funktion erwartet lediglich einen Parameter "size", der die Shape der Ausgabe definiert. Wenn wir "size" bspw. mit (3, 4) angeben, erhalten wir ein Array mit der Shape (3, 4), die mit den Zufalls-Werten gefüllt sind:

In [ ]:
import numpy as np
x = np.random.random_sample((3, 4))
print(x)

Wenn random_sample mit einem Integer-Wert aufgerufen wird, erhalten wir ein eindimensionales Array. Ein Integer-Wert bewirkt den gleichen Effekt, wie ein einfaches Tupel als Argument:

In [ ]:
x = np.random.random_sample(7)
print(x)
y = np.random.random_sample((7,))
print(y)

Es können ebenfalls Arrays aus einem beliebigen Intervall generiert werden, bei dem a kleiner als b sein muss. Das kann wie folgt aussehen:

(b - a) * random_sample() + a

Beispiel:

In [ ]:
a = -3.4
b = 5.9
A = (b - a) * np.random.random_sample((3, 4)) + a
print(A)

Das Standard-Modul "random" von Python hat eine allgemeinere Funktion "sample", die Samples einer Population produziert. Die Population kann eine Sequenz oder eine Menge sein.

Die Syntax von "sample":

sample(population, k)

Die Funktion erstellt eine Liste, die k Elemente beinhaltet aus "population". Die Ergebnis-Liste beinhaltet keine Duplikate, wenn in der Population keine Duplikate vorkommen.

Wenn eine Sample aus einer Range von Integer-Werten ausgewählt werden soll, dann können - oder besser sollten Sie - range als Argument für die Population verwenden.

Im folgenden Beispiel ziehen wir sechs Zahlen aus der Range zwischen 1 und 49 (inklusive). Das entspricht der einer Lotto-Ziehung in Deutschland:

In [ ]:
import random
print(random.sample(range(1, 50), 6))
In [ ]:
def weighted_sample(population, weights, k):
    """ 
    This function draws a random sample of length k 
    from the sequence 'population' according to the 
    list of weights
    """
    sample = set()
    population = list(population)
    weights = list(weights)
    while len(sample) < k:
        choice = weighted_choice(population, weights)
        sample.add(choice)
        index = population.index(choice)
        weights.pop(index)
        population.remove(choice)
        weights = [ x / sum(weights) for x in weights]
    return list(sample)
In [ ]:
def weighted_sample_alternative(population, weights, k):
    """ 
    Alternative way to previous implementation.
        
    This function draws a random sample of length k 
    from the sequence 'population' according to the 
    list of weights
    """
    sample = set()
    population = list(population)
    weights = list(weights)
    while len(sample) < k:
        choice = weighted_choice(population, weights)
        if choice not in sample:
            sample.add(choice)
    return list(sample)

Beispiel:

Nehmen wir an, wir haben acht Zuckerstücke in den Farben rot, grün, blau, gelb, schwarz, weiß, pink, und orange. Unser Freund Peter hat für die Farben die "gewichteten" Vorlieben 1/24, 1/6, 1/6, 1/12, 1/12, 1/24, 1/8, 7/24. Peter darf sich 3 Zuckerstücke aussuchen:

In [ ]:
balls = ["red", "green", "blue", "yellow", "black", 
         "white", "pink", "orange"]
weights = [ 1/24, 1/6, 1/6, 1/12, 1/12, 1/24, 1/8, 7/24]
for i in range(10):
    print(weighted_sample(balls, weights, 3))

Erhöhen wir die Wahrscheinlichkeit, dass ein orangenes Zuckerstück im Sample enthalten ist:

In [ ]:
n = 100000
orange_counter = 0
orange_counter_alternative = 0
for i in range(n):
    if "orange" in weighted_sample(balls, weights, 3):
        orange_counter += 1
    if "orange" in weighted_sample_alternative(balls, weights, 3):
        orange_counter_alternative += 1 
        
print(orange_counter / n)
print(orange_counter_alternative / n)



Kartesische Auswahl

Die Funktion cartesian_choice ist nach dem Kartesischen Produkt aus der Mengen-Lehre benannt.

Kartesisches Produkt

Das kartesische Produkt ist eine Berechnung die eine Menge zurückliefert aus mehreren Mengen. Die Ergebnis-Menge des kartesischen Produkts wird auch "Produkt-Menge" oder einfach "Produkt" genannt.

Für die beiden Mengen A und B ist das kartesische Produkt A × B die Menge aller geordneten Paare (a, b) für die a ∈ A and b ∈ B gilt:

A x B = { (a, b) | a ∈ A and b ∈ B }

Wenn wir n Mengen haben A1, A2, ... An, können wir das kartesische Produkt entsprechend wie folgt bilden:

A1 x A2 x ... x An = { (a1, a2, ... an) | a1 ∈ A1, a2 ∈ A2, ... an ∈ An]

Das kartesische Produkt aus n Mengen wir manchmal auch n-faches kartesisches Produkt genannt.

Kartesische Auswahl: cartesian_choice

Wir schreiben nun eine Funktion cartesian_choice, die eine beliebige Anzahl von iterierbaren Argumenten erwartet und eine Liste zurückliefert, die eine zufällige Auswahl von jedem Iterator beinhaltet in der entsprechenden Reihenfolge.

Aus mathematischer Sicht können wir das Ergebnis der Funktion cartesian_choice als Element des kartesischen Produkts der übergebenen iterierbaren Argumente ansehen.

In [ ]:
import random
def cartesian_choice(*iterables):
    """
    A list with random choices from each iterable of iterables 
    is being created in respective order.
    
    The result list can be seen as an element of the 
    Cartesian product of the iterables 
    """
    res = []
    for population in iterables:
        res.append(random.choice(population))
    return res
cartesian_choice(["The", "A"],
                 ["red", "green", "blue", "yellow", "grey"], 
                 ["car", "house", "fish", "light"],
                 ["smells", "dreams", "blinks"])

Wir definieren nun eine gewichtete Version der vorher definierten Funktion:

In [ ]:
import random
def weighted_cartesian_choice(*iterables):
    """
    A list with weighted random choices from each iterable of 
    iterables is being created in respective order
    """
    res = []
    for population, weight in iterables:
        lst = weighted_choice(population, weight)
        res.append(lst)
    return res
determiners = (["The", "A", "Each", "Every", "No"], 
               [0.3, 0.3, 0.1, 0.1, 0.2])
colours = (["red", "green", "blue", "yellow", "grey"], 
           [0.1, 0.3, 0.3, 0.2, 0.2])
nouns = (["water", "elephant", "fish", "light", "programming language"], 
         [0.3, 0.2, 0.1, 0.1, 0.3])
nouns2 = (["of happiness", "of chocolate", "of wisdom", 
           "of challenges", "of air"], 
         [0.5, 0.2, 0.1, 0.1, 0.1])
verb_phrases = (["smells", "dreams", "thinks", "is made of"], 
         [0.4, 0.3, 0.3])
print("It may or may not be true:")
for i in range(10):
    res = weighted_cartesian_choice(determiners,
                                    colours,
                                    nouns,
                                    verb_phrases,
                                    nouns2)
    print(" ".join(res) + ".")

In der folgenden Version prüfen wir, ob alle "Wahrscheinlichkeiten" korrekt sind:

In [ ]:
import random
def weighted_cartesian_choice(*iterables):
    """
    A list with weighted random choices from each iterable of iterables 
    is being created in respective order
    """
    res = []
    for population, weight in iterables:
        lst = weighted_choice(population, weight)
        res.append(lst)
    return res
determiners = (["The", "A", "Each", "Every", "No"], 
               [0.3, 0.3, 0.1, 0.1, 0.2])
colours = (["red", "green", "blue", "yellow", "grey"], 
           [0.1, 0.3, 0.3, 0.2, 0.2])
nouns = (["water", "elephant", "fish", "light", "programming language"], 
         [0.3, 0.2, 0.1, 0.1, 0.3])
nouns2 = (["of happiness", "of chocolate", "of wisdom", 
           "of challenges", "of air"], 
         [0.5, 0.2, 0.1, 0.1, 0.1])
verb_phrases = (["smells", "dreams", "thinks", "is made of"], 
         [0.4, 0.3, 0.2, 0.1])
print("It may or may not be true:")
sentences = []
for i in range(10000):
    res = weighted_cartesian_choice(determiners,
                                    colours,
                                    nouns,
                                    verb_phrases,
                                    nouns2)
    sentences.append(" ".join(res) + ".")
words = ["smells", "dreams", "thinks", "is made of"]
from collections import Counter
c = Counter()
for sentence in sentences:
    for word in words:
        if word in sentence:
            c[word] += 1
  
wsum = sum(c.values())
for key in c:
    print(key, c[key] / wsum)

Wahre Zufallszahlen

Haben Sie jemals ein Würfel-Spiel gespielt und sich gefragt, ob etwas mit den Würfeln nicht stimmt? Sie würfeln so oft, jedoch erscheint bspw. niemals die 6.

Vielleicht haben Sie sich auch gefragt, ob das random-Modul Python "wahre" Zufallszahlen generieren kann, z.B. wie einen richtigen Würfel. In Wahrheit ist, dass die meisten Zufallszahlen in Computer-Programmen pseudozufällig sind. Die Zahlen werden in einer vorhersehbaren Art generiert, weil der Algorithmus deterministisch ist. Pseudozufällige Zahlen sind für viele Situation gut genug, aber ist nicht "wirklich" zufällig für Würfel oder Lotterie-Ziehungen.

Die Webseite RANDOM.ORG erhebt den Anspruch "wahre" Zufallszahlen zu generieren. Sie nutzen die Zufälligkeit aus atmosphärischen Geräuschen. Diese Zahlen sind für viele Anwendungen besser als die pseudozufälligen Zahlen aus einem Algorithmus in Computer-Programmen.



Zufalls-Samen

Random Seed

Ein Zufalls-Samen, - auch "Samen-Status" oder einfach "Samen" genannt - ist ein Zahl um einen pseudozufälligen Zahlen-Generator zu initialisieren.

Wenn wir random.random() aufrufen, so erwarten wir einen Zufalls-Wert zwischen 0 und 1. random.random() berechnet eine neuen Zufalls-Wert anhand des vorangegangenen Zufalls-Wertes. Und wenn wir random zum ersten Mal in unserem Programm verwenden? Genau, es gibt keinen vorangegangenen Zufalls-Wert. Wenn ein Zufalls-Generator zum ersten Mal aufgerufen wird, so muss der erste "Zufalls"-Wert zuerst erstellt werden.

Wenn wir einen Pseudo-Zufallszahlen-Generator säen, so liefern wir bereits einen ersten "vorangegangen" Wert. Der Samen-Wert korrespondiert mit einer Sequenz von generierten Werten eines gegebenen Zufalls-Generators. Wenn Sie den gleichen Samen-Wert wieder verwenden, so erhalten Sie wieder die gleiche Sequenz an Werten.

Die Samen-Zahl selbst muss nicht zufällig sein, so dass der Algorithmus Werte erzeugt, die einer Wahrscheinlichkeitsverteilung in einer pseudozufälligen Weise folgen. Momentan ist der Samen für die Sicherheit von Bedeutung. Wenn Sie den Samen kennen, können Sie bspw. den Sicherheits-Schlüssel generieren, der auf dem Samen basiert.

Zufällige Samen werden in vielen Programmiersprachen anhand des Systemstatus generiert, der meistens die System-Zeit ist.

Für Python trifft dies ebenfalls zu. help(random.seed) sagt, dass wenn die Funktion mit None oder keinen Argumenten aufgerufen wird, so wird anhand "der aktuellen System-Zeit oder einer anderen System-Spezifischen Zufalls-Quelle" gesät.

In [ ]:
import random
help(random.seed)

Die seed-Funktion liefert eine deterministische Sequenz von Zufalls-Zahlen. Die Sequenz kann beliebig wiederholt werden, um bspw. bestimmte Situationen zu debuggen.

In [ ]:
import random
random.seed(42)
for _ in range(10):
    print(random.randint(1, 10), end=", ")
    
print("\nLet's create the same random numbers again:")
random.seed(42)
for _ in range(10):
    print(random.randint(1, 10), end=", ")



Zufalls-Zahlen in Python mit der Gaussschen und Normalvariaten Verteilung

Wir möchten nun 1000 Zufalls-Zahlen zwischen 130 und 230 generieren die eine Gausssche Verteilung aufweisen mit dem Hauptwert mu = 550 und eine Standard-Abweichung von sigma = 30.

In [ ]:
from random import gauss
n = 1000
values = []
frequencies = {}
while len(values) < n:
    value = gauss(180, 30)
    if 130 < value < 230:
        frequencies[int(value)] = frequencies.get(int(value), 0) + 1
        values.append(value)
        
print(values[:10])

Das folgende Programm zeichnet die Zufalls-Werte, die wir soeben generiert haben. Wir haben matplotlib bis jetzt noch nicht behandelt. Das ist aber für das Verständnis des folgenden Codes nicht notwendig:

In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
freq = list(frequencies.items())
freq.sort()
plt.plot(*list(zip(*freq)))

Jetzt führen wir dies mit der Normalvariaten Verteilung durch statt mit der Gaussschen:

In [ ]:
from random import normalvariate
n = 1000
values = []
frequencies = {}
while len(values) < n:
    value = normalvariate(180, 30)
    if 130 < value < 230:
        frequencies[int(value)] = frequencies.get(int(value), 0) + 1
        values.append(value)
freq = list(frequencies.items())
freq.sort()
plt.plot(*list(zip(*freq)))



Übung mit 0en und 1en

Es ist vielleicht keine schlechte Idee folgende Funktion als Übung selbst zu schreiben. Die Funktion soll mit einem Parameter p aufgerufen werden, der einen Wahrscheinlichkeitswert zwischen 0 und 1 beinhaltet. Die Funktion liefert eine 1 mit der Wahrscheinlichkeit von p, d.h. in p Prozent der Fälle werden 1en zurückgegeben und 0en in (1-p) Prozent der Fälle:

In [ ]:
import random
def random_ones_and_zeros(p):
    """ p: probability 0 <= p <= 1
        returns a 1 with the probability p
    """
    x = random.random()
    if x < p:
        return 1
    else:
        return 0
    

Wir testen unsere kleine Funktion:

In [ ]:
n = 1000000
sum(random_ones_and_zeros(0.8) for i in range(n)) / n

Eine weitere gute Idee ist es, die Aufgabe mit einem Generator zu implementieren. Sollten Sie nicht mit der Arbeitsweise eines Python-Generators vertraut sein, empfehlen wir unser Kapitel zu Generatoren und Iteratoren aus dem Python-Tutorial.

In [ ]:
import random
def random_ones_and_zeros(p):
    while True:
        x = random.random()
        yield 1 if x < p else 0
        
def firstn(generator, n):
	for i in range(n):
		yield next(generator)
In [ ]:
n = 1000000
sum(x for x in firstn(random_ones_and_zeros(0.8), n)) / n

Unser 0en-und-1en-Generator kann wie ein Sender betrachtet werden, der 0en und 1en mit einer Wahrscheinlichkeit von p, respektive (1-p), abgibt.

Wir schreiben nun einen anderen Generator, der diesen Bitstrom empfängt. Die Aufgabe dieses neuen Generators ist es, den Bitstrom zu lesen und einen anderen Bitstrom aus 0en und 1en zu erzeugen mit einer Wahrscheinlichkeit von 0.5 ohne die Wahrscheinlichkeit p zu kennen. Es sollte für jeden beliebigen Wert p funktionieren.1

Bitstreams Lighthouses In [ ]:
def ebitter(bitstream):
    while True:
        bit1 = next(bitstream)
        bit2 = next(bitstream)
        if bit1 + bit2 == 1:
            bit3 = next(bitstream)
            if bit2 + bit3 == 1:
                yield 1
            else:
                yield 0
        
In [ ]:
def ebitter2(bitstream):
    bit1 = next(bitstream)
    bit2 = next(bitstream)
    bit3 = next(bitstream)
    while True:
        if bit1 + bit2 == 1:
            if bit2 + bit3 == 1:
                yield 1
            else:
                yield 0
        bit1, bit2, bit3 = bit2, bit3, next(bitstream)
In [ ]:
n = 1000000
sum(x for x in firstn(ebitter(random_ones_and_zeros(0.8)), n)) / n
In [ ]:
n = 1000000
sum(x for x in firstn(ebitter2(random_ones_and_zeros(0.8)), n)) / n

Grundlagen der Theorie:

Unser erster Generator erzeugt einen Bitstrom B0, B1, B2,...

Wir prüfen nun ein beliebiges Paar aufeinanderfolgender Bits Bi, Bi+1,...

Ein solches Paar kann die Werte 01, 10, 00 oder 11 haben. Die Wahrscheinlichkeit p(01) = (p-1) x p und die Wahrscheinlichkeit p(10) = p x (p-1) ergibt eine kombinierte Wahrscheinlichkeit, dass die aufeinanderfolgenden Bits entweder 01 oder 10 sind, von 2 x (p-1) x p.

Betrachten wir ein anderes Bit Bi+2. Wie is die Wahrscheinlichkeit dieser beiden

Bi + Bi+1 = 1

und

Bi+1 + Bi+2 = 1 ?

Die möglichen Ausgaben passen auf die Bedingungen und die zugehörigen Wahrscheinlichkeiten sind in der folgenden Tabelle aufgeführt:


Wahrscheinlichkeit Bi Bi+1 Bi+2
p2 x (1-p) 0 1 0
p x (1 - p)2 1 0 1


Wir bezeichnen das Ergebnis sum(Bi, Bi+1)=1 als X1 und entsprechend dazu sum(Bi+1, Bi+2)=1 als X2.

Die verknüpfte Wahrscheinlichkeit P(X1, X2) = p2 x (1-p) + p x (1 - p)2 kann zu p x (1-p) umgestellt werden.

Die bedingte Wahrscheinlichkeit von X2 gibt X1:

P(X2 | X1) = P(X1, X2) / P(X2)

P(X2 | X1) = p x (1-p) / 2 x p x (1-p) = 1 / 2



Synthetische Verkaufszahlen

In diesem Unterkapitel möchten wir eine Datei mit Verkaufszahlen erzeugen. Stellen Sie sich vor wir hätten eine Kette von Geschäften in diversen europäischen und kanadischen Städten: Frankfurt, München, Berlin, Zürich, Hamburg, London, Toronto, Strassburg, Luxemburg, Amsterdam, Rotterdam,

Wir beginnen mit einem Array "sales" mit Verkaufszahlen für das Jahr 1997:

In [ ]:
import numpy as np
sales = np.array([1245.89, 2220.00, 1635.77, 1936.25, 1002.03, 2099.13,  723.99, 990.37, 541.44, 1765.00, 1802.84, 1999.00])

Das Ziel ist eine kommaseparierte Liste zu erzeugen, wie wir sie aus Excel kennen. Die Datei soll die Verkaufszahlen enthalten, die wir nicht kennen, für alle Geschäfte, die wir nicht haben, über die Jahre von 1997 bis 2016.

Wir fügen für jedes Jahr den Verkaufszahlen noch zufällige Werte hinzu. Dafür konstruieren wir ein Array mit Wachstums-Raten. Die Wachstums-Raten können variieren zwischen einem minimalen Prozent-Wert (min_percent) und einem maximalen Prozent-Wrt (max_percent):

In [ ]:
min_percent = 0.98  # corresponds to -1.5 %
max_percent = 1.06   # 6 %
growthrates = (max_percent - min_percent) * np.random.random_sample(12) + min_percent
print(growthrates)

Für die neuen Verkaufszahlen nach einem Jahr multiplizieren wir das "sales"-Array mit dem "growthrates"-Array:

In [ ]:
sales * growthrates

Um eine nachhaltige Verkaufs-Entwicklung zu erhalten, verändern wir die Wachstums-Raten alle 5 Jahre.

Das ist unser komplettes Porgramm, welches die Daten in der Datei sales_figures.csv ablegt:

In [ ]:
import numpy as np
fh = open("sales_figures.csv", "w")
fh.write("Year, Frankfurt, Munich, Berlin, Zurich, Hamburg, London, Toronto, Strasbourg, Luxembourg, Amsterdam, Rotterdam, The Hague\n")
sales = np.array([1245.89, 2220.00, 1635.77, 1936.25, 1002.03, 2099.13,  723.99, 990.37, 541.44, 1765.00, 1802.84, 1999.00])
for year in range(1997, 2016):
    line = str(year) + ", " + ", ".join(map(str, sales))
    fh.write(line + "\n")
    if year % 4 == 0:
         min_percent = 0.98  # corresponds to -1.5 %
         max_percent = 1.06   # 6 %
         growthrates = (max_percent - min_percent) * np.random.random_sample(12) + min_percent
         #growthrates = 1 + (np.random.rand(12) * max_percent - negative_max) / 100
    sales = np.around(sales * growthrates, 2)
fh.close()

Das Ergebnis ist in der Datei sales_figures.csv zu finden.

Wir werden diese Daten in einem folgenden Kapitel (Lesen und Schreiben in Numpy) noch benutzen.



Übungen

  1. Machen wir etwas mehr mit Würfeln. Beweisen Sie empirische - indem Sie ein Simulations-Programm schreiben - dass die Wahrscheinlichkeit für das kombinierte Ereignis "egal welche Zahl wurde gewürfelt" (E) und "Eine Zahl größer als 2 wurde gewürfelt" 1/3 ist.
  2. Die Datei ["universities_uk.txt"](universities_uk.txt) beinhaltet eine Liste der Universitäten im Vereinigten Königreich zur Einschreibung zwischen 2013-2014. (Quelle: [Wikipedia](https://en.wikipedia.org/wiki/List_of_universities_in_the_United_Kingdom_by_enrollment#cite_note-1)). Schreiben Sie eine Funktion die ein Tupel (universities, enrollments, total_number_of_students) zurückliefert mit:
    • universities: Liste der Namen der Universitäten
    • enrollments: zugehöirge Liste mit Einschreibungen
    • total_number_of_students: Über alle Universitäten

    Jetzt können Sie 100,000 fiktive Studenten immatrikulieren mit einer Wahrscheinlichkeit die den "echten" Einschreibungen entspricht.
  3. Ancient Greek scene as Pythonia

    Lasst uns eine kleine Reise in die Zeit und den Raum der nächsten Übung machen. Wir reisen zurück in das antike Pythonia (Πηθωνια). Es war in der Zeit, in der der König Pysseus als gütiger Diktator regierte. Es war die Zeit, als Pysseus seine Botschafter aussandte um durch die Welt zu reisen und zu verkünden, dass es für seine Prinzen Anacondos (Ανακονδος), Cobrion (Κομπριον), Boatos (Μποατος) und Addokles (Ανδοκλης) an der Zeit ist zu heiraten. Sie veranstalteten ein Programmier-Wettbewerb zwischen den fairen und tapferen Amazonenen, besser bekannt als Pythanier aus Pythonia. Zum Schluss wurden nur 11 Amazonen ausgewählt:

    1. Die ätherische Airla (Αιρλα)
    2. Barbara (Βαρβαρα), die Eine aus dem fremden Land
    3. Eos (Ηως), Sieht in der Dämmerung göttlich aus
    4. Die süße Glykeria (Γλυκερια)
    5. Die anmutige Hanna (Αννα)
    6. Helen (Ελενη), das Licht in der Dunkelheit
    7. Der gute Engel Agathangelos (Αγαθαγγελος)
    8. Die violett getönte Wolke Iokaste (Ιοκαστη)
    9. Medousa (Μεδουσα), Die Wächterin
    10. Die selbstbestimmende Sofronia (Σωφρονια)
    11. Andromeda (Ανδρομεδα), Die eine die wie eine Mann oder ein Krieger denkt

    Der Tag, an dem Sie die Chance bekamen in der Lotterie gezogen zu werden, war für alle Amazonen der gleiche. Aber Pysseus wollte die Lotterie ein paar Tage in die Zukunft verschieben. Die Wahrscheinlichkeit änderte sich mit jedem Tag: Sie sank um 1/13 für die ersten 7 Amazonen und stieg um 1/12 für die letzten 5 Amazonen.

    Wie lange muss der König warten, bis er zu 90% sicher sein kann, dass seine Prinzen Anacondas, Cobrion, Boatos und Addokles die Amazonen Iokaste, Medouse, Sofronia und Andromeda heiraten?

    </li>

    </ol>



    Lösungen zu den Übungen

    1. In [ ]:
      from random import randint
      outcomes = [ randint(1, 6) for _ in range(10000)]
      even_pips = [ x for x in outcomes if x % 2 == 0]
      greater_two = [ x for x in outcomes if x > 2]
      combined = [ x for x in outcomes if x % 2 == 0 and x > 2]
      print(len(even_pips) / len(outcomes))
      print(len(greater_two) / len(outcomes))
      print(len(combined) / len(outcomes))
      

    2. Wir schreiben zuerst die Funktion "process_datafile" um die Daten aus der Datei zu verarbeiten:

      In [ ]:
      def process_datafile(filename):
          """ process_datafile -> (universities, 
                                   enrollments, 
                                   total_number_of_students) 
              universities: list of University names
              enrollments: corresponding list with enrollments
              total_number_of_students: over all universities
          """
          universities = []
          enrollments = []
          with open(filename) as fh:
              total_number_of_students = 0
              fh.readline() # get rid of descriptive first line
              for line in fh:
                   line = line.strip()
                   *praefix, undergraduates, postgraduates, total = line.rsplit()
                   university = praefix[1:]
                   total = int(total.replace(",", ""))
                   enrollments.append(total)
                   universities.append(" ".join(university))
                   total_number_of_students += total
          return (universities, enrollments, total_number_of_students)
      

      Lassen wir die Funktion laufen und prüfen das Ergebnis:

      In [ ]:
      universities, enrollments, total_students = process_datafile("universities_uk.txt")
      for i in range(14):
          print(universities[i], end=": ")
          print(enrollments[i])
      print("Total number of students onrolled in the UK: ", total_students)
      

      Wir wollen einen virtuellen Studenten in eine zufällige Universität einschreiben. Um eine gewichtete Liste zu erhalten, die für die weighted_choice geeignet ist, müssen wir die Werte der Liste "enrollments" normalisieren:

      In [ ]:
      normalized_enrollments = [ students / total_students for students in enrollments]
      # enrolling a virtual student:
      print(weighted_choice(universities, normalized_enrollments))
      

      Die Aufgabe war 100,000 fiktive Studenten zu "immatrikulieren". Dies kann mit einer Schleife einfach durchgeführt werden:

      In [ ]:
      from collections import Counter
      outcomes = []
      n = 100000
      for i in range(n):
          outcomes.append(weighted_choice(universities, normalized_enrollments))
      c = Counter(outcomes)
          
      print(c.most_common(20))
      

      </li>

    3. Übung 3

      Die Sammlung der Amazonen ist als Liste implementiert, während für die Menge aus Pysseusses Favoritinnen auswählen. Die Gewichtung liegt zu Beginn bei 1/11 für alle, dh. 1/len(amazons).

      Jeder Schleifen-Durchlauf entspricht einem neuen Tag. Jedes Mal wenn wir einen neuen Durchlauf starten, ziehen wir "n" Samples aus den Pythoniern um das Verhältnis zu berechnen wue oft die Sample gleich den Favoritinnen des Königs geteilgt durch die Häufigkeit, wie oft die Sample nicht der Idee einer Schwiegertochter entspricht. Dies entspricht der Wahrscheinlichkeit "prob". Wir stoppen das erste Mal, wenn die Wahrscheinlichkeit bei 0.9 oder größer liegt.

      Für die Ziehung können wir beide Funktionen "weighted_same" und "weighted_same_alternative" verwenden.

      In [ ]:
      import time
      amazons = ["Airla", "Barbara", "Eos",
                 "Glykeria", "Hanna", "Helen",
                 "Agathangelos", "Iokaste", 
                 "Medousa", "Sofronia", 
                 "Andromeda"]
      weights = [ 1/len(amazons) for _ in range(len(amazons)) ]
      Pytheusses_favorites = {"Iokaste", "Medousa", 
                              "Sofronia", "Andromeda"}
      n = 1000
      counter = 0
      prob = 1 / 330
      days = 0
      factor1 = 1 / 13
      factor2 = 1 / 12
      start = time.clock()
      while prob < 0.9:
          for i in range(n):
              the_chosen_ones = weighted_sample_alternative(amazons, weights, 4)
              if set(the_chosen_ones) == Pytheusses_favorites:
                  counter += 1
          prob = counter / n
          counter = 0
          weights[:7] = [ p - p*factor1 for p in weights[:7] ]
          weights[7:] = [ p + p*factor2 for p in weights[7:] ]
          weights = [ x / sum(weights) for x in weights]
          days += 1
      print(time.clock() - start)
      print("Number of days, he has to wait: ", days)
      

      Die Werte für die Anzahl der Tage weichen ab, wenn n nicht groß genug ist.

      Der folgende Code ist die Lösung ohne Abrundungs-Fehlern. Wir verwenden Fraction aus dem Modul fractions.

      In [ ]:
      import time
      from fractions import Fraction
      amazons = ["Airla", "Barbara", "Eos",
                 "Glykeria", "Hanna", "Helen",
                 "Agathangelos", "Iokaste", 
                 "Medousa", "Sofronia", 
                 "Andromeda"]
      weights = [ Fraction(1, 11) for _ in range(len(amazons)) ]
      Pytheusses_favorites = {"Iokaste", "Medousa", 
                              "Sofronia", "Andromeda"}
      n = 1000
      counter = 0
      prob = Fraction(1, 330)
      days = 0
      factor1 = Fraction(1, 13)
      factor2 = Fraction(1, 12)
      start = time.clock()
      while prob < 0.9:
          #print(prob)
          for i in range(n):
              the_chosen_ones = weighted_sample_alternative(amazons, weights, 4)
              if set(the_chosen_ones) == Pytheusses_favorites:
                  counter += 1
          prob = Fraction(counter, n)
          counter = 0
          weights[:7] = [ p - p*factor1 for p in weights[:7] ]
          weights[7:] = [ p + p*factor2 for p in weights[7:] ]
          weights = [ x / sum(weights) for x in weights]
          days += 1
      print(time.clock() - start)
      print("Number of days, he has to wait: ", days)
      

      Wir können sehen dass die Lösung mit fractions schön aber langsam ist. Dabei spielt die Präzision in unserem Fall keine Rolle.

      Jedoch haben wir die Leistung von PYthon nicht genutzt. Das machen wir in der nächsten Implementierung:

      In [ ]:
      import time
      import numpy as np
      amazons = ["Airla", "Barbara", "Eos",
                 "Glykeria", "Hanna", "Helen",
                 "Agathangelos", "Iokaste", 
                 "Medousa", "Sofronia", 
                 "Andromeda"]
      weights = np.full(11, 1/len(amazons))
      Pytheusses_favorites = {"Iokaste", "Medousa", 
                              "Sofronia", "Andromeda"}
      n = 1000
      counter = 0
      prob = 1 / 330
      days = 0
      factor1 = 1 / 13
      factor2 = 1 / 12
      start = time.clock()
      while prob < 0.9:
          for i in range(n):
              the_chosen_ones = weighted_sample_alternative(amazons, weights, 4)
              if set(the_chosen_ones) == Pytheusses_favorites:
                  counter += 1
          prob = counter / n
          counter = 0
          weights[:7] = weights[:7] - weights[:7] * factor1
          weights[7:] = weights[7:] + weights[7:] * factor2
          weights = weights / np.sum(weights)
          #print(weights)
          days += 1
      print(time.clock() - start)
      print("Number of days, he has to wait: ", days)
      

      </li> </ol>





      Fußnoten:


      1 Ich bedanke mich bei Dr. Hanno Baehr der mich auf das Problem der "Zufalls-Extrahierung" aufmerksam gemacht hat bei der Teilnahme eine Python-Training-Kurses in Nürnberg im Januar 2014. Hanno hat ein paar Bits des theoretischen Rahmens entworfen. Während einer Nacht-Session in der Bar "Zeit & Raum" (engl. "Time & Space") habe ich ein entsprechendes Python-Programm implementiert um die theoretische Lösung empirisch zu unterstützen.

      In [ ]: