EINFÜHRUNG

Eine der wichtigsten Fragestellungen überhaupt besteht darin zu verstehen oder vorauszusagen, wie sich ein System im Laufe der Zeit verhält. Dies ist umso schwieriger, je mehr Einflussparametern das System ausgesetzt ist. Will man das System mit dem Computer simulieren, so kann dies nur in diskreten Zeitschritten geschehen, was für sich kontinuierlich verändernde Systeme immer nur eine Approximation darstellt, die zwar mit kleineren Zeitschritten immer genauer wird, wobei sich aber der Rechenaufwand dabei massiv vergrössern kann.

Es gibt viele Systeme, bei denen der Übergang vom aktuellen Systemzustand Z_i zum Nachfolgezustand Z_kdurch Wahrscheinlichkeiten geregelt ist und nur vom Zustand i abhängt (also nicht von noch weiter zurückliegenden). Willst du beispielsweise mit einem Würfel alle 6 Zahlen erwürfeln, so kannst du diesen Prozess als Durchlauf einer Reihe von Zuständen Z₀ bis Z₆ auffassen und grafisch als Markoff-Kette darstellen.

Die Übergangswahrscheinlichkeiten p_i vom Zustand i zum Nachfolgezustand kannst du im folgenden Schema eintragen:

MITTLERE WARTEZEITEN

Statt nach Wahrscheinlichkeiten zu fragen, willst du abklären, wie oft du im Mittel würfeln musst, um von einem bestimmten Anfangszustand in einen bestimmten Endzustand zu gelangen. Wenn die in gleichen Zeitschritten würfelst, so fragst du dich also, wie lange im Mittel der Prozess dauert. Diese Zeit nennst du die (mittlere) Wartezeit. Offensichtlich ist beim vorher beschriebenen Würfelprozess die Wartezeit um von Z₀ zu Z₆ die Summe der Wartezeiten für alle Übergänge. Wie gross sind aber die einzelnen Wartezeiten?

In deinem ersten Programm findest du experimentell die wichtigste Eigenschaft von Wartezeitproblemen heraus.:

Ist p die Wahrscheinlichkeit um von Z₁ nach Z₂ zu gelangen, so beträgt die Wartezeit (in einer geeigneten Einheit) u = 1/p.

In der Simulation untersuchst du die Wartezeit, um eine bestimmte Augenzahl, sagen wir eine 6, zu würfeln. Ein Simulationsversuch besteht hier nicht aus einem einzigen Würfelwurf. Vielmehr wirfst du in der Funktion sim() soviel Mal, bis du du die 6 erhalten hast und gibst die Anzahl der dazu nötigen Würfe zurück. Diesen Versuch wiederholst du 10000 Mal und bildest das Mittel der benötigten Würfe. Gleichzeitig stellst du in einem Häufigkeitsdiagramm die Zahl Versuche dar, bei denen du in k = 1, 2, 3,... Würfen die 6 erhieltest (bei k = 50 hörst du auf).

from gpanel import *
from random import randint
n = 10000
p = 1/6

def sim():
    k = 1
    r = randint(1, 6)
    while r != 6:
        r = randint(1, 6)
        k += 1
    return k

makeGPanel(-5, 55, -200, 2200)
drawGrid(0, 50, 0, 2000)
title("Waiting on a 6")
h = [0] * 51
lineWidth(5)
count = 0
repeat n:
    k = sim()
    count += k
    if k <= 50:
        h[k] += 1
        line(k, 0, k, h[k])
mean_exp = count / n

lineWidth(1)
setColor("red")
count = 0
for k in range(1, 1000):
    pk = (1 - p)**(k - 1) * p
    nk = n * pk
    count += nk * k
    if k <=50:
        line(k, 0, k, nk)
mean_theory = count / n
title("Experiment: " + str(mean_exp) + "Theory: " + str(mean_theory))

MEMO

AUSBREITUNG EINER KRANKHEIT

Du gehst von der folgenden Geschichte aus, die zwar erfunden ist, aber durchaus gewisse Parallelen zu aktuellen Lebensgemeinschaften hat.

"Auf einer einsamen Karibikinsel leben abgeschnitten von der Außenwelt 100 Menschen. Ein alter Mann infiziert sich durch den Verzerr eines Zugvogels, den er gedankenlos gegessen hat, mit einer ansteckenden Krankheit. Trifft ein kranker Einwohner mit einem gesunden zusammen, so wird dieser in kurzer Zeit krank. Jede Stunde treffen sich zwei Menschen zufällig." Mit einer Computersimulation willst du untersuchen, wie sich die Krankheit ausbreitet. Dazu bestimmst du die Zahl der Infizierten in Abhängigkeit von der Zeit.

Es ist elegant, die Population durch eine Liste mit booleschen Werten zu modellieren, wobei gesund mit False und krank mit True codiert sind. Der Vorteil dieser Datenstruktur besteht darin, dass sich die Wechselwirkung beim Zusammentreffen von zwei Personen in der Funktion pair() mit einer logischen OR-Verknüpfung ausdrücken lässt:

from random import randint

n = 10000

def sim():
    z = []
    i = 0
    while True:
        r = randint(1, 6)
        i += 1
        if not r in z:
            z.append(r)
        if len(z) == 6:
            return i

count = 0
repeat n:
    count += sim()

print("Mean waiting time:", count / n)

MEMO

Interessant ist die Frage, wie lange es im Mittel geht, bis alle Einwohner krank sind. Du kannst das Problem direkt mit einer Computersimulation lösen, in der du vielmals dieselbe Population simulierst und dabei die Schritte zählst, bis alle krank sind.

from gpanel import *
from random import randint

def pair():
    # Select two distinct inhabitants
    a = randint(0, 99)
    b = a
    while b == a:
        b = randint(0, 99)
    z[a] = z[a] or z[b]
    z[b] = z[a]

def nbInfected():
    count = 0
    for i in range(100):
        if z[i]:
            count += 1
    return count

makeGPanel(-50, 550, -10, 110)
title("The spread of an illness")
drawGrid(0, 500, 0, 100)
lineWidth(2)
setColor("blue")

z = [False] * 100
tmax =  500
t = 0
a = randint(0, 99) 
z[a] = True # random infected inhabitant
move(t, 1)

while t <= tmax:
    pair()
    infects = nbInfected()
    t += 1
    draw(t, infects)

Du kannst die Ausbreitung der Krankheit aber auch als Markoff-Kette auffassen. Ein bestimmter Zustand ist durch die Zahl der Infizierten charakterisiert. Die Wartezeit, bis alle krank sind, ist die Summe der Wartezeiten für die Übergänge von k Kranken zu k+1 Kranken für k von 1 bis 99. Dazu benötigst du die Wahrscheinlichkeit p_k für diesen Übergang. Es gilt:

pk = Summe der Wahrscheinlichkeiten als ersten einen Kranken und als zweiten einen Gesunden auszuwählen und umgekehrt: Im Programm stellst du pk auch noch gerade grafisch dar und ermittelst die Summe der Reziprokwerte.

from random import randint

n = 1000 # number experiment

def pair():
    # Select two distinct inhabitants
    a = randint(0, 99)
    b = a
    while b == a:
        b = randint(0, 99)
    z[a] = z[a] or z[b]
    z[b] = z[a]

def nbInfected():
    count = 0
    for i in range(100):
        if z[i]:
            count += 1
    return count

def sim():
    global z
    z = [False] * 100
    t = 0
    a = randint(0, 99) 
    z[a] = True # random infected inhabitant
    while True:
        pair()
        t += 1
        if nbInfected() == 100:
            return t

count = 0
for i in range(n):
    u = sim()
    print("Experiment #", i + 1, "Waiting time:", u)
    count += u
    
print( "Mean waiting time:", count / n)

MEMO

WARTEZEITPARADOXON

Das Warten an einer Haltestelle von öffentlichen Verkehrsmitteln (Bus, Tram, S-Bahn, U-Bahn) gehört zum Alltag. Wir wollen uns hier überlegen, wie gross die mittlere Wartezeit ist, falls du zu völlig zufälliger Zeit an die Haltestelle eines Busses kommst (also den Fahrplan nicht kennst). Dabei gehst du zuerst davon aus, dass die Busse exakt im Takt von 6 Minuten hintereinander fahren. Es ist einleuchtend, dass du machmal nur sehr kurz und manchmal maximal fast 6 Minuten warten musst. Im Mittel beträgt die Wartezeit also 3 Minuten. Wie steht es aber, wenn die Busse nicht mehr exakt im Takt fahren, sondern beispielsweise mit gleich verteilter Wahrscheinlichkeit im Bereich von 2 und 10 Minuten?

Da sie auch dann noch im Mittel alle 6 Minuten vorbei kommen, nimmst du vielleicht an, dass sich an der mittleren Wartezeit von 3 Minuten nichts ändert. Die erstaunliche, und darum paradoxe Antwort ist, dass die mittlere Wartezeit nun grösser als 3 Minuten wird.

In einer animierten Simulation sollst du herausfinden, wie gross die Wartezeit unter der Annahme ist, dass die Busse gleichverteilt im Bereich 2 und 10 Minuten (in der Simulation sind es Sekunden) hintereinander fahren. Dabei verwendest du mit Vorteil die Game-Library JGameGrid ein, da sich die beteiligten Objekte, wie Busse und Fahrgäste als Sprite-Objekte modellieren lassen.

Der Programmcode erfordert einige Erklärungen:

• xxx
• xxx

from gpanel import *

n  = 100

def p(k):
    return 2 * k * (n - k) / n / (n - 1)

makeGPanel(-10, 110, -0.1, 1.1)
drawGrid(0, 100, 0, 1.0)

count = 0
for k in range(1, n - 1):
    if k == 1:
        move(k, p(k))
    else:
        draw(k, p(k))
    count += 1 / p(k)

title("Time until everyone is ill: " + str(count))

MEMO

AUFGABEN