Python, Numpy and Wahrscheinlichkeiten


Einführung

dice

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

Diese Aussage ist auf Bill Clinton zurückzuführen über diverse Webseiten. Sie steht normalerweise in keinem Kontext und somit ist nicht klar, ob es nicht möglicherweise als "Witz" gemeint ist. Was auch immer seine Intention war, wir haben ihn zitiert um ein Statistik-Beispiel aus dem "wahren" Leben zu zeigen. Statistik und Wahrscheinlichkeites-Berechnung ist überall in unserem Leben. Wir müssen damit umgehen, immer wenn wir eine Entscheidung treffen müssen mit verschiedenen Möglichkeiten. Können wir heute Abend wandern gehen oder wird es regnen? Die Wettervorhersage sagt uns, dass die Niederschlags-Wahrscheinlichkeit bei 30% liegt. Was jetzt? Werden wir wandern gehen?
Andere Situation: Sie sind in den Ferien auf einer paradiesischen Insel, weit weg von zu Hause. Plötzlich treffen Sie ihren Nachbarn! Wie hoch standen die Chancen?
Sie spielen jede Woche Lotto und träumen von der weit entfernten Insel. Wie hoch ist die Wahrscheinlichkeit dei 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

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

Wie können wir ein Rollen des Würfels in Python simulieren? Dafür brauchen wir Numpy nicht. Das "reine" Python und die random-Module reichen aus.

import random
outcome = random.randint(1,6)
print(outcome)
1

Lass uns unseren virtuellen Würfel 10 mal werfen:

import random
[ random.randint(1, 6) for _ in range(10) ]
Der obige Code führt zu folgendem Ergebnis:
[4, 3, 2, 1, 5, 5, 6, 6, 5, 3]

Einfacher geht es mit der random-Funktion aus dem Numpy-Modul.

import numpy as np
outcome = np.random.randint(1, 7, size=10)
print(outcome)
[3 6 5 1 1 4 4 3 3 1]

Sie werden bemerkt haben, dass statt der 6 eine 7 verwendet als zweiter Parameter wurde. randint von numpy.random nutzt ein "halb-offenes" Intervall im Gegensatz zu Pythons random-Modul, dass ein geschlossenes Intervall nutzt.

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
[6]
[3 4 1 4 1 4 3 4 4 2]
[4 6 2 5 6 6 4 2 3 5]
[[4 4 5 4]
 [4 4 5 2]
 [3 3 3 4]
 [3 5 4 5]
 [6 1 2 6]]

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:

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

"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
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 3 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, wie weighted_choice funktioniert:

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:

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.0883, 0.1625, 0.1626, 0.165, 0.1677, 0.2539]

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:

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

import numpy as np
x = np.random.random_sample((3, 4))
print(x)
[[ 0.79766358  0.1071769   0.95808703  0.23360052]
 [ 0.08147494  0.11206388  0.53209231  0.22282572]
 [ 0.90766065  0.35183364  0.6290983   0.45014664]]

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.0306952   0.40171768  0.51830863  0.69956378  0.70096778  0.36403344
  0.04985113]
[ 0.1085724   0.44630135  0.48393855  0.7821603   0.29624174  0.4241912
  0.62388892]

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:

a = -3.4
b = 5.9
A = (b - a) * np.random.random_sample((3, 4)) + a
print(A)
[[  3.59028683e+00  -2.22520276e-01  -2.41319577e+00   4.49734751e+00]
 [  1.23262351e+00  -3.24489578e+00  -1.21536358e+00   5.39019638e+00]
 [  5.31013142e+00  -1.39074768e+00   2.62175282e-03  -7.61831153e-02]]

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:

import random
print(random.sample(range(1, 50), 6))
[4, 24, 40, 11, 1, 44]
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))
['blue', 'pink', 'orange']
['blue', 'green', 'black']
['blue', 'green', 'orange']
['yellow', 'black', 'red']
['pink', 'black', 'orange']
['blue', 'green', 'orange']
['blue', 'black', 'orange']
['yellow', 'white', 'orange']
['yellow', 'pink', 'red']
['blue', 'yellow', 'green']

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

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)
0.71147
0.71398



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.

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 erhalten die folgende Ergebnisse:
['The', 'grey', 'house', '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 = (["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) + ".")
It may or may not be true:
No blue elephant thinks of challenges.
The blue programming language thinks of happiness.
The grey elephant smells of happiness.
The grey light thinks of chocolate.
A yellow programming language dreams of wisdom.
The green elephant dreams of chocolate.
A red programming language smells of wisdom.
The green water smells of happiness.
A yellow fish thinks of chocolate.
Every blue water smells 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.0991
thinks 0.1996
smells 0.4043
dreams 0.297

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.

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.
    
    For version 2 (the default), all of the bits are used if *a* is a str,
    bytes, or bytearray.  For version 1, the hash() of *a* is used instead.
    
    If *a* is an int, all bits are used.

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, 



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.

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])
[173.49123947564414, 183.47654360102564, 186.96893210720162, 214.90676059797428, 199.69909520396007, 183.31521532331496, 157.85035192965537, 149.56012897536849, 187.39026585633607, 219.33242481612143]

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:

%matplotlib inline
import matplotlib.pyplot as plt
freq = list(frequencies.items())
freq.sort()
plt.plot(*list(zip(*freq)))
Wir können die folgenden Ergebnisse erwarten, wenn wir den obigen Python-Code ausführen:
[<matplotlib.lines.Line2D at 0x7f3d9b720cc0>]

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

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)))
Führt man obigen Code aus, erhält man folgendes Ergebnis:
[<matplotlib.lines.Line2D at 0x7f3d9b67d048>]



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

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
Führt man obigen Code aus, erhält man Folgendes:
0.800609

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
sum(x for x in firstn(random_ones_and_zeros(0.8), n)) / n
Wir können die folgenden Ergebnisse erwarten, wenn wir den obigen Python-Code ausführen:
0.799762

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
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 Python-Code führt zu folgender Ausgabe:
0.49975
n = 1000000
sum(x for x in firstn(ebitter2(random_ones_and_zeros(0.8)), n)) / n
Wir können die folgende Ausgabe erwarten, wenn wir den obigen Python-Code ausführen:
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 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:

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

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)
[ 1.01492188  1.01480476  1.05093035  1.05994539  1.00390292  1.02079311
  1.03427277  1.00683137  0.98921442  0.99205428  1.03891395  1.04466069]

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

sales * growthrates
Der obige Code liefert folgendes Ergebnis:
array([ 1264.48102414,  2252.86656873,  1719.08033715,  2052.31926251,
        1005.94084529,  2142.77743504,   748.80314228,   997.13558072,
         535.60025636,  1750.97580145,  1872.99562131,  2088.276718  ])

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:

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. 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))
      
      0.5061
      0.6719
      0.3402
      

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

      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:

      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)
      
      Open University in England: 123490
      University of Manchester: 37925
      University of Nottingham: 33270
      Sheffield Hallam University: 33100
      University of Birmingham: 32335
      Manchester Metropolitan University: 32160
      University of Leeds: 30975
      Cardiff University: 30180
      University of South Wales: 29195
      University College London: 28430
      King's College London: 27645
      University of Edinburgh: 27625
      Northumbria University: 27565
      University of Glasgow: 27390
      Total number of students onrolled in the UK:  2299380
      

      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:

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

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

      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))
      
      [('Open University in England', 5316), ('University of Manchester', 1705), ('Sheffield Hallam University', 1450), ('University of Leeds', 1409), ('University of Birmingham', 1383), ('Manchester Metropolitan University', 1380), ('University of Nottingham', 1352), ('Cardiff University', 1335), ("King's College London", 1257), ('University College London', 1242), ('University of Glasgow', 1206), ('University of Edinburgh', 1201), ('University of South Wales', 1199), ('Northumbria University', 1196), ('University of Sheffield', 1174), ('University of Plymouth', 1172), ('University of Central Lancashire', 1170), ('University of the West of England', 1167), ('Coventry University', 1159), ('University of Hertfordshire', 1146)]
      

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

      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)
      
      1.9674369999999897
      Number of days, he has to wait:  32
      

      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:

      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)
      
      1.6973969999999952
      Number of days, he has to wait:  33
      

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