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.9227335362865872

Nun wollen wir eine kryptografisch starke Zufallszahl erzeugen:

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

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(3, secure=False))
[0.2372999369862978, 0.399237260441624, 0.4165971640424626]

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)
323 µs ± 5.58 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
random_list(100, secure=False)
7.67 µs ± 72.3 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)
Wir können die folgende Ausgabe erwarten, wenn wir den obigen Python-Code ausführen:
array([ 0.19713999,  0.08812961,  0.64416305,  0.94776954,
        0.00616254,  0.34216703,  0.29579779,  0.60599821,
        0.96039013,  0.31544036])
%%timeit
np.random.random(100)
1.59 µs ± 13.8 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, also [start, 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)]
Führt man obigen Code aus, erhält man folgendes Ergebnis:
[4, 4, 6, 3, 5, 3, 3, 3, 4, 5, 3, 5, 6, 6, 3, 5, 3, 4, 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)
4
import secrets
c = secrets.SystemRandom()
c.randint(1, 6)
Führt man obigen Code aus, erhält man folgendes Ergebnis:
4

Würfeln wir sehr oft, so sollte jede Augenzahl ungefähr in einem Sechstel der Fälle 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({2: 16927, 4: 16753, 5: 16743, 3: 16546, 1: 16531, 6: 16500})
Augenzahl: 4,           relative Häufigkeit: 0.16753
Augenzahl: 1,           relative Häufigkeit: 0.16531
Augenzahl: 5,           relative Häufigkeit: 0.16743
Augenzahl: 3,           relative Häufigkeit: 0.16546
Augenzahl: 2,           relative Häufigkeit: 0.16927
Augenzahl: 6,           relative Häufigkeit: 0.165

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.16735
Augenzahl: 5,           relative Häufigkeit: 0.16816
Augenzahl: 6,           relative Häufigkeit: 0.16743
Augenzahl: 3,           relative Häufigkeit: 0.16415
Augenzahl: 2,           relative Häufigkeit: 0.1668
Augenzahl: 1,           relative Häufigkeit: 0.16611

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, das ein geschlossenes Intervall erwartet.

Die formale Definition:

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

Diese Funktion liefert zufällige ganze Zahlen 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 für high der Wert None oder nichts übergeben wird, liegen die Ergebnisse in der Range von [0, "low").

Der Parameter "size" definiert die Shape des Ergebnisses. Wenn für size None oder nichts übergeben wird, generiert die Funktion einen Integer-Wert. 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("Eine Integer-Zahl: ", 
      np.random.randint(1, 7))
print("\nEin eindimensionales Array mit einem Element:\n",
      np.random.randint(1, 7, size=1))
print("\nEin eindimensionales Array mit zehn Elementen:\n", 
      np.random.randint(1, 7, size=10))
print("\nWie eben, aber alternative size-Definition:\n", 
      np.random.randint(1, 7, size=(10,))) 
print("\nZweidimensionales Array:\n", 
      np.random.randint(1, 7, size=(5, 4)))
Eine Integer-Zahl:  6
Ein eindimensionales Array mit einem Element:
 [1]
Ein eindimensionales Array mit zehn Elementen:
 [4 2 2 1 4 2 1 4 3 1]
Wie eben, aber alternative size-Definition:
 [2 3 3 4 4 5 2 3 2 4]
Zweidimensionales Array:
 [[5 4 3 2]
 [1 5 2 2]
 [2 4 2 6]
 [3 3 2 4]
 [4 5 6 1]]

Wir können das Würfeln 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 bei 1/6 liegt.

Wie können wir einen gezinkten Würfen 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.



Stichproben / Auswahlen


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.

Sequenzen können beispielsweise Listen, Strings und Tupels sein, aber auch Iteratoren. Das bedeutet, wir sind in der Lage aus einem String ein zufälliges Zeichen zu holen oder ein zufälliges Element aus einer Liste oder einem Tupel:

from random import choice
professions = ["scientist", "philosopher", "engineer", "priest"]
print(choice("abcdefghij"))
print(choice(professions))
print(choice(("apples", "bananas", "cherries")))
print(choice(range(10)))
f
philosopher
apples
8

"choice" liefert zufällig ein Objekt aus einer nicht-leeren Sequenz, wobei die Chancen für die Elemente ausgewählt zu werden gleichmäßig verteilt sind. So liegt die Chance für die Rückgabe von "scientist" beim Aufruf von choice(profession) bei 1/4. Das hat allerdings wenig mit der Realität zu tun. Um die Realität nachzubilden benötigen wir auch hier - wie beim gezinkten Würfel - eine gewichtete Auswahl.

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


Zufallsintervalle

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
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=", ")
-1, 0, 0, 1, 2, 2, 3, 4, 6, 6, -1, -1, 



Gewichtete Zufalls-Auswahl

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

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

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, auf der weighted_choice basiert:

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: Stellen wir uns vor, dass wir einen gezinkten Würfel haben, sodass für die Wahrscheinlichkeiten des Auftretens der Augenzahlen 6 und 1 gilt: P(6)=3/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 mit den Augenzahlen des Würfels und der Liste der korrespondierenden Gewichte auf. Jeder Aufruf entspricht einem Wurf des gezinkten Würfels.

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

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()))
[0.0772, 0.1615, 0.1669, 0.1674, 0.1734, 0.2536]

Die Werte der Aufteilungs-Liste definieren die Intervalle, 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 Unterteilung als ein Intervall von

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:

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=", ")
Endpoints are included:
-1, 0, 0, 1, 2, 2, 3, 6, 6, -1, -1, 
Endpoints are not included:
0, 1, 1, 2, 3, 3, 4, 7, 7, 8, 8, 



Zufalls-Samples mit Python


Eine Sample oder Stichprobe 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 im halb-offenen Intervall [0.0, 1.0) liefert. 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 die size zum Beispiel mit (3, 4) angeben, erhalten wir ein Array mit der Shape (3, 4), das mit Zufalls-Werten gefüllt sind:

import numpy as np
x = np.random.random_sample((3, 4))
print(x)
[[ 0.8264495   0.6310656   0.20705883  0.12143256]
 [ 0.68678488  0.57697705  0.89330325  0.76748047]
 [ 0.81193228  0.31446597  0.21202368  0.40172387]]

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:

x = np.random.random_sample(7)
print(x)
y = np.random.random_sample((7,))
print(y)
[ 0.90397582  0.74251759  0.85259247  0.21515893  0.1109198
  0.44684381  0.76229266]
[ 0.51943675  0.72439752  0.61401382  0.64369677  0.45759271
  0.75285897  0.60699391]

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

(b - a) * random_sample() + a

Beispiel:

a = -3.4
b = 2
A = (b - a) * np.random.random_sample((3, 4)) + a
print(A)
[[-3.36441019 -3.04912123 -3.26656291  1.47122998]
 [ 0.19256774  0.53046134  0.79345875 -2.73161851]
 [-0.42905873  1.28669361 -3.07284688  0.90472292]]

Das Standard-Modul "random" von Python hat eine allgemeinere Funktion "sample", die Samples einer Population produziert. Die Population ist eine Sequenz also beispielsweise eine Liste, Menge oder ein String.

Die Syntax von "sample":

sample(population, k)

Die Funktion erstellt eine Liste, die k Elemente aus der population beinhaltet. 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 einer Lotto-Ziehung in Deutschland:

import random
print(random.sample(range(1, 50), 6))
[27, 29, 43, 46, 3, 36]
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)
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:

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))
['yellow', 'blue', 'orange']
['blue', 'orange', 'white']
['green', 'orange', 'pink']
['blue', 'orange', 'pink']
['green', 'blue', 'black']
['green', 'orange', 'black']
['yellow', 'green', 'orange']
['green', 'blue', 'orange']
['yellow', 'blue', 'orange']
['yellow', 'green', 'orange']

Im folgenden vergleichen wir die beiden Varianten gewichtete Stichproben zu berechnen:

n = 10000
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)
0.7092
0.7107



Kartesische Auswahl

Die Funktion cartesian_choice, die wir im Folgenden definieren werden, ist nach dem Kartesischen Produkt aus der Mengenlehre benannt.

Kartesisches Produkt

Das kartesische Produkt, auch Mengenprodukt genannt, ist ein Konstruktionsprinzip um aus Mengen eine neue Menge zu konstruieren.

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

A x B = { (a, b) | a &isin; A and b &isin; B }

Ganz allgemein besteht das kartesische Produkt mehrerer Mengen aus der Menge aller Tupel von Elementen der Mengen, wobei die Reihenfolge der Mengen und damit der entsprechenden Elemente fest vorgegeben ist. Die Ergebnismenge des kartesischen Produkts wird auch Produktmenge, Kreuzmenge oder Verbindungsmenge genannt.

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.

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"])
Der obige Code führt zu folgendem Ergebnis:
['A', 'blue', 'light', 'blinks']

Wir definieren nun eine gewichtete Version der vorher definierten Funktion:

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 = (["door", "elephant", "fish", "light", 
          "programming language", "Python"], 
         [0.2, 0.2, 0.1, 0.1, 0.3, 0.1])
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", "consists"], 
         [0.1, 0.3, 0.3, 0.2, 0.1])
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) + ".")
It may or may not be true:
No yellow Python is made of happiness.
A blue fish dreams of happiness.
Every red elephant thinks of challenges.
No yellow Python dreams of happiness.
Each yellow programming language thinks of happiness.
The grey light smells of happiness.
A grey fish is made of wisdom.
A green programming language dreams of happiness.
No blue elephant dreams of happiness.
A green programming language dreams of happiness.

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

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)
It may or may not be true:
is made of 0.1007
smells 0.3998
dreams 0.3011
thinks 0.1984

Echte Zufallszahlen

Zur Erzeugung echter Zufallszahlen kann man physikalische Phänomene nutzen. MAn kann sie aus radioaktiven Zerfallsprozessen oder dem Rauschen elektronischer Bauelemente gewinnen. Aber es geht auch einfacher zum Beispiel mit dem Werfen einer oder mehrerer Münzen oder unter Benutzung von Würfeln. Diese Verfahren sind jedoch einerseite zeitaufwendig oder technisch kompliziert zu bewerkstelligen.

Wir wir bereits geschrieben haben, liefert das random-Modul von Python keine "echten" oder "wahren" Zufallszahlen. Die meisten Programme liefern nur Zufallszahlen die pseudozufällig sind. Die Zahlen werden in einer vorhersehbaren Art generiert, weil der verwendete 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.



Seed Key

Random Seed

Unter einem "Seed Key" oder einfach nur Seed (deutsch wörtlich "Saatschlüssel") versteht man einen Wert, mit dem ein Zufallszahlengenerator initialisiert wird. Deshalb bezeichnet man ihn häufig auch als Startwert einer Zufallsfolge. Der Zufallszahlengenerator erzeugt mit der Seed als Startwert eine Folge von Zufallszahlen bzw. Pseudozufallszahlen. Startet man einen Algorithmus mit dem gleichen Seed-Wert, erhält man auch exakt die gleichen Pseudozufallszahlen.

Verwendet man Zufallszahlen in der Kryptographie, um damit verschlüsselte Daten zu erzeugen, so möchte man normalerweise nicht, dass die Seed erraten werden kann. In diesem Fall sollte die Seed ein Zufallswert sein. Diesen kann man beispielsweise in einem Programm durch das Hin- und Herbewegen einer Maus erzeugen.

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. Wie sieht es aber aus, wenn wir diese Funktion zum ersten Mal in unserem Programm verwenden? Genau, dann gibt es noch keinen vorangegangenen Zufalls-Wert. Wenn ein Zufalls-Generator zum ersten Mal aufgerufen wird, so muss der erste "Zufalls"-Wert irgendwie erzeugt werden.

Setzt man den Startwert, also die Seed, eines Pseudo-Zufallszahlen-Generators, so stellt man einen ersten Zufallswert zur Verfügung. Aus diesem wird dann deterministisch eine Sequenz von weiteren Zufallswerten generiert. Startet man wieder eine Zufallsfolge mit dem gleichen Startwert, so erhält man wieder exakt die selbe Folge von Werten. Kennt man also den Startwert, kann man auch die aus diesem erzeugten verschlüsselten Werte berechnen.

Geht es also um Sicherheit, braucht man einen Weg, einen echten Zufallswert als Startwert zu verwenden.

Zufällige Startwerte werden in vielen Programmiersprachen anhand des Systemstatus generiert, der meistens der Systemzeit entspricht.

Für Python trifft dies ebenfalls zu. help(random.seed) sagt, dass wenn die Funktion mit None oder keinem Argument aufgerufen wird der Seedwert aus der aktuellen Systemzeit oder einer anderen System-Spezifischen Zufalls-Quelle" generiert wird.

import random
help(random.seed)
Help on method seed in module random:
seed(a=None, version=2) method of random.Random instance
    Initialize internal state from hashable object.
    
    None or no argument seeds from current time or from an operating
    system specific randomness source if available.
    
    If *a* is an int, all bits are used.
    
    For version 2 (the default), all of the bits are used if *a* is a str,
    bytes, or bytearray.  For version 1 (provided for reproducing random
    sequences from older versions of Python), the algorithm for str and
    bytes generates a narrower range of seeds.

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

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=", ")
2, 1, 5, 4, 4, 3, 2, 9, 2, 10, 
Let's create the same random numbers again:
2, 1, 5, 4, 4, 3, 2, 9, 2, 10, 



random.gauss und random.normalvariate

Bei random.gauss und random.normalvariate handelt es sich um zwei verschiedene Implementierungen von Funktionen, die jeweils normalverteilte Zufallszahlen zurückliefern.

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.

Wenn wir uns die help-Information anschauen, sehen wir, dass der wesentliche Unterschied darin besteht, dass gauss etwas schneller ist und nicht Thread-sicher!

import random
help(random.normalvariate)
help(random.gauss)
Help on method normalvariate in module random:
normalvariate(mu, sigma) method of random.Random instance
    Normal distribution.
    
    mu is the mean, and sigma is the standard deviation.
Help on method gauss in module random:
gauss(mu, sigma) method of random.Random instance
    Gaussian distribution.
    
    mu is the mean, and sigma is the standard deviation.  This is
    slightly faster than the normalvariate() function.
    
    Not thread-safe without a lock around calls.

Wir berechnen uns nun mittels der Funktion gaus 1000 Zufallzahlen einem Mittelwert von 180 und einer Standardabweichung von 30:

from random import gauss
n = 1000
values = []
frequencies = {}
while len(values) < n:
    value = int(gauss(180, 30))
    if 130 < value < 230:
        frequencies[value] = frequencies.get(value, 0) + 1
        values.append(value)
values[:7]    
Führt man obigen Code aus, erhält man folgende Ausgabe:
[173, 183, 186, 214, 199, 183, 157]

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 unbedingt notwendig:

%matplotlib inline
import matplotlib.pyplot as plt
freq = list(frequencies.items())
freq.sort()
plt.plot(*list(zip(*freq)))
Führt man obigen Code aus, erhält man folgendes Ergebnis:
[<matplotlib.lines.Line2D at 0x7f3dda110748>]

Wir können dies natürlich auch mit der Funktion normalvariate durchführen:

from random import normalvariate
n = 1000
values = []
frequencies = {}
while len(values) < n:
    value = int(normalvariate(180, 30))
    if 130 < value < 230:
        frequencies[value] = frequencies.get(value, 0) + 1
        values.append(value)
freq = list(frequencies.items())
freq.sort()
plt.plot(*list(zip(*freq)))
Führt man obigen Code aus, erhält man folgende Ausgabe:
[<matplotlib.lines.Line2D at 0x7f3dda041828>]



Übung mit Binärsender

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:

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:

n = 1000000
sum(random_ones_and_zeros(0.8) for i in range(n)) / n
Der obige Python-Code führt zu folgender Ausgabe:
0.800606

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.

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)
n = 1000000
firstn_values = firstn(random_ones_and_zeros(0.8), n)
sum(x for x in firstn_values) / n
Der obige Code führt zu folgendem Ergebnis:
0.799768

Bitstreams Lighthouses

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

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
        
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)
n = 1000000
sum(x for x in firstn(ebitter(random_ones_and_zeros(0.8)), n)) / n
Der obige Code führt zu folgendem Ergebnis:
0.499749
n = 1000000
sum(x for x in firstn(ebitter2(random_ones_and_zeros(0.8)), n)) / n
Führt man obigen Code aus, erhält man folgendes Ergebnis:
0.500011

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 geht es um die synthetische Erzeugung von Daten. Dazu erzeugen wir eine Datei mit Verkaufszahlen einer fiktiven Kette von Geschäften in diversen europäischen Städten.

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

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 fiktive Verkaufszahlen für unsere nicht-existierenden Ladenlokale für die Jahre von 1997 bis 2017.

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):

min_percent = 0.98  # corresponds to -1.5 %
max_percent = 1.06   # 6 %
sample = np.random.random_sample(12)
print(sample)
growthrates = (max_percent - min_percent) * sample + min_percent
print(growthrates)
[ 0.19199923  0.74882856  0.32861167  0.92521988  0.67732252
  0.13142505  0.62833576  0.52683706  0.53815011  0.13553688
  0.37325259  0.52726136]
[ 0.99535994  1.03990628  1.00628893  1.05401759  1.0341858
  0.990514    1.03026686  1.02214697  1.02305201  0.99084295
  1.00986021  1.02218091]

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

sales * growthrates
Führt man obigen Code aus, erhält man folgende Ausgabe:
array([ 1240.10899394,  2308.59195249,  1646.05724873,
        2040.84155935,  1036.28519895,  2079.21766154,
         745.90290466,  1012.30368982,   553.9212797 ,
        1748.83780752,  1820.61637616,  2043.33963726])

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

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

import numpy as np
fh = open("sales_figures.csv", "w")
fh.write("Year, Frankfurt, Munich, Berlin, Zurich, Hamburg, \
         London, Paris, Luxembourg, Wien, 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, 2018):
    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 %
         sample = np.random.random_sample(12)
         growthrates = (max_percent - min_percent) * sample + min_percent         
    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.

Lösung 2. Aufgabe

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

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

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:

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

Lösung 3. Aufgabe

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.perf_counter()
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.perf_counter() - 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 Rundungsfehler. 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.perf_counter()
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.perf_counter() - 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.perf_counter()
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.perf_counter() - start)
print("Number of days, he has to wait: ", days)





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 [ ]: