Höchstwahrscheinlich (Linux-Magazin, Juli 2014)

Mathematische Rätsel mit bedingten Wahrscheinlichkeiten knackt der Fachmann mit der Bayes-Formel oder auch mit diskreten Verteilungen, erzeugt von kurzen Perl-Skripts.

Mit dem in Statistikerkreisen bekannten Ziegenproblem [6] und seiner verblüffenden Lösung kann der Rätselfreund reihenweise mathematische Laien hereinlegen, aber auch so mancher Fachmann hat sich daran schon damit blamiert, weil er leidenschaftlich gegen die richtige Lösung gewettert hat. Wer hätte gedacht, dass sich Wahrscheinlichkeiten in offensichtlich festgezimmerten Fernsehstudios dramatisch ändern, nur weil der Moderator eine weitere Tür ohne Preis öffnet?

Das menschliche Gehirn scheint sich allgemein schwer zu tun mit sogenannten bedingten Wahrscheinlichkeiten, die zufällige Ereignisse mit Vorbedingungen beschreiben. So auch bei der vor 15 Jahren im Perl-Snapshot ([4]) schon einmal vorgestellten Aufgabe: "Gegeben ist ein Zylinder, in dem drei Karten liegen: Die erste ist vorne und hinten schwarz, die zweite vorne und hinten rot, und die dritte auf der einen Seite schwarz und auf der anderen rot. Eine Karte wird gezogen. Man sieht die Vorderseite, und die ist rot. Wie hoch ist die Wahrscheinlichkeit daß auch die Rückseite der gezogenen Karte rot ist?"

Abbildung 1: Der Satz von Bayes (von der Wikipedia)

Intuitiv falsch

Die meisten Leute antworten darauf mit "50%", denn scheinbar gibt es zwei gleichwahrscheinliche Fälle, die rot-schwarze und die rot-rote Karte, und einmal ist die Kehrseite schwarz und einmal rot. Steckt der Rätselfreund dem erstaunten Zuhörer aber dann, dass die korrekte Lösung "66%" lautet, macht sich oft Entrüstung breit. Der Kniff liegt an den Vorbedingungen des Experiments: Die rot-rote Karte wird einfach doppelt so häufig gezogen wie die rot-schwarze, weil das Experiment als Vorbedingung "eine Seite ist rot" stellt.

Abbildung 2: Aus einem Zylinder mit einer schwarz/roten, einer schwarz/schwarzen und einer rot/roten Karte zieht ein Proband Karten.

Schon vor 250 Jahren hat sich der Mathematiker Thomas Bayes daran gemacht, solche gefühlsmäßig schwer zu verstehenden Probleme in eine Formel zu fassen, an der sich nicht rütteln lässt. Der "Satz von Bayes" ([5]) beschreibt in seiner "diachronischen Interpretierung", wie sich die Wahrscheinlichkeiten von Hypothesen im Lauf der Zeit verändern, falls in einem Experiment neue Daten auftauchen. Die Formel

    P(H|D) = P(H) * P(D|H) / P(D)

definiert die Wahrscheinlichkeit einer Hypothese H (zum Beispiel "ich habe die rot-rote Karte gezogen") in Abhängigkeit von neu auftauchenden Daten D ("die Vorderseite der gezogenen Karte ist rot"). Auf der rechten Seite ist P(H) dabei die Ausgangswahrscheinlichkeit der Hypothese, noch bevor neue Daten vorliegen. Multipliziert wird sie mit P(D|H), also der Wahrscheinlichkeit, dass die Daten unter der Hypothese vorliegen. Wie wahrscheinlich ist es, dass der Kartenfreund auf eine rote Vorderseite starrt, falls er die rot-rote Karte gezogen hat? Und im Nenner der rechten Seite steht schließlich mit P(D) die Wahrscheinlichkeit der vorliegenden Daten, unabhängig von irgendeiner Hypothese. Wie wahrscheinlich ist es, dass der Proband nach dem Ziehen irgendeiner Karte auf eine rote Vorderseite starrt?

Bayes schlägt Schnippchen

Alle drei möglichen Hypothesen (zufällig gezogene Karten) sind in diesem Experiment offensichtlich gleich wahrscheinlich. Also ist P(H), also die Wahrscheinlichkeit, dass die rot-rote Karte gezogen wird, gleich 1/3, denn sie wird in einem Drittel aller Fälle gezogen. Und bei Annahme dieser Hypothese ist die Kehrseite der gezogenen rot-roten Karte trivialerweise in exakt 100% aller Fälle ebenfalls rot, also ist P(D|H) = 1. Unabhängig von der Hypothese einer gezogenen Karte ist die Wahrscheinlichkeit, nach dem Ziehen auf eine rote Kartenoberseite zu blicken 50%, also ist P(D) = 1/2. Schließlich liegen im Zylinder 6 Kartenhälften, von denen drei rot und drei schwarz sind. In die Bayes-Formel eingesetzt ergibt sich für P(H|D) entsprechend 1/3 * 1 / (1/2) = 2/3. Die Bayes-Formel sagt also mit 66% das im Experiment nachweisbar korrekte Ergebnis voraus, und schlägt damit der falsch liegenden Intuition ein Schnippchen.

Diskret programmiert

Das etwa ein halbes Jahr alte Buch "Think Bayes" [2] beschreibt nun zusätzlich zur Bayes-Formel ein numerisches Verfahren, das sich zum Experimentieren leicht in handliche Skripts fassen lässt. Dazu legt man die Hypothesen in einer statistischen Verteilung ab, speichert sie also in einem Python-Dictionary oder Perl-Hash mit ihren Ausgangswahrscheinlichkeiten ab, multipliziert dann die Einzelwerte mit den bedingten Wahrscheinlichkeiten eintreffender Daten unter Annahme der jeweiligen Hypothese, normalisiert dann die Werte für alle Hypothesen, und erhält die Wahrscheinlichkeit der gesuchten Hypothese im Lichte neu eingetroffener Daten. Dieses Verfahren, und speziell die abschließende Normalisierung funktioniert allerdings nur, falls:

    1) höchstens eine der Hypothesen wahr ist und
    2) es keine weiteren Möglichkeiten gibt, also
       mindestens eine der Hypothesen wahr ist

Listing 1 zeigt die Implementierung des Skripts cards, das am Ende die gesuchte Lösung "0.666.." ausgibt. Es nutzt das in Listing 2 vorgestellte Modul Distrib.pm, das die Hilfsfunktionen zur Berechnung der statistischen Verteilung bereitstellt. Die Wahrscheinlichkeiten für die Hypothesen der Ziehung verschiedener Karten ("BB", "BR", "RR" für black-black, black-red, red-red) füttern die Zeilen 7-9 mit jeweils dem Wert 1/3 in das Modul. Die hinzugekommenen Daten, also die Oberseite der tatsächlich gezogenen Karte, geben die Zeilen 11-13 für jede Hypothese vor und multiplizieren die eingespeisten Hypothesenwerte damit. So beträgt die Wahrscheinlichkeit, dass sich dem Proband eine rote Oberseite präsentiert bei der schwarz-schwarzen Karte gleich Null. Bei der schwarz-roten Karte beträgt sie 50%, je nachdem wie herum die Karte aus dem Zylinder herauskommt. Und die rot-rote Karte wiederum präsentiert dem Laborassistenten in 100% aller Fälle eine rote Seite, also setzt Zeile 13 den Multiplikator auf den Wert 1. Abschließend gibt der Aufruf der Methode prob("RR") die normalisierte Wahrscheinlichkeit für die Hypothese "Rot-rote Karte" zurück, die 2/3 beträgt.

Listing 1: cards

    01 #!/usr/local/bin/perl -w
    02 use strict;
    03 use Distrib;
    04 
    05 my $distrib = Distrib->new();
    06 
    07 $distrib->set( "BB", 0.33 );
    08 $distrib->set( "BR", 0.33 );
    09 $distrib->set( "RR", 0.33 );
    10 
    11 $distrib->mult( "BB", 0 );
    12 $distrib->mult( "BR", 0.5 );
    13 $distrib->mult( "RR", 1 );
    14 
    15 $distrib->normalize();
    16 
    17 print $distrib->prob( "RR" ), "\n";

Moose als Gimmick

Listing 3 nutzt das CPAN-Modul Moose, um sich die Definition des Distrib-Konstruktors in Perl zu ersparen, allerdings muss der Code statt dessen mit "has" die Klassenattribute deklarieren und initialisieren. Im vorliegenden Fall nur ein netter Gimmick, bei weiteren Attributen wäre der Code deutlich schlanker als bei manueller Klassenfahrt.

Listing 2: Distrib.pm

    01 use Moose;
    02 
    03 has 'values' => (is => 'rw', isa => 'HashRef', 
    04                  default => sub { {} } );
    05 
    06 sub set {
    07     my( $self, $hypo, $prob ) = @_;
    08     
    09     $self->values()->{ $hypo } = $prob;
    10 }
    11 
    12 sub mult {
    13     my( $self, $hypo, $prob ) = @_;
    14 
    15     $self->values()->{ $hypo } *= $prob;
    16 }
    17 
    18 sub normalize {
    19     my( $self ) = @_;
    20 
    21     my $values = $self->values();
    22     my $sum    = 0;
    23 
    24     for my $hypo ( keys %$values ) {
    25         $sum += $values->{ $hypo };
    26     }
    27     for my $hypo ( keys %$values ) {
    28         $values->{ $hypo } /= $sum;
    29     }
    30 }
    31 
    32 sub prob {
    33     my( $self, $hypo ) = @_;
    34 
    35     return $self->values()->{ $hypo };
    36 }
    37 
    38 1;

Die Methode set() nimmt den Namen einer Hypothese (zum Beispiel "BB") und deren A-priori-Wahrscheinlichkeit entgegen, und speichert sie in dem Objekt-internen Array values. Zum Multiplizieren eines Hypothesenwerts mit einem konstanten Wert nimmt mult() beide entgegen, findet den bislang gespeicherten Wert in values und multipliziert ihn mit dem hereingereichten Wert für $prob. Die Methode normalize() iteriert über alle bislang eingespeisten Werte im Hash, summiert sie in $sum auf und dividiert dann alle Werte durch diesen Wert. So beträgt die neue Summe aller Wahrscheinlichkeitswerte zum Beispiel nach einer Multiplikation wieder 1, jeder Wert kann also als Wahrscheinlichkeit zwischen 0 und 1 interpretiert werden. Am Ende findet prob() den Wert für eine gesuchte Hypothese, indem es den Werte-Hash mittels der von Moose automatisch erzeugten Methode values() als Referenz hervorholt und unter dem Schlüssel der gesuchten Hypothese den dort abgelegten Wert extrahiert.

Abstrakter geht's auch

Führt man mehrere solcher Tests für verschiedenartigste Probleme durch, zeigt sich ein wiederkehrendes Muster: Nach dem Aufsetzen der Hypothesen werden die Wahrscheinlichkeiten zusätzliche Daten jeweils mit allen definierten Hypothesen multipliziert. Es bietet sich also wie in [2] vorgemacht an, eine von Distrib abgeleitete Musterklasse HypoTest wie in Listing 2 zu definieren, die mittels einer Methode update() die Werte für alle Hypothesen auffrischt. Weiter verlässt sie sich darauf, dass wiederum von ihr abgeleitete Klassen (wie z.B. CardHypoTest in Listing 4) die abstrakte Methode likelihood() überladen und basierend auf der übergebenen Hypothese und den zusätzlich verfügbaren Daten die Wahrscheinlichkeit P(D|H) zurückgeben. Das Framework HypoTest ruft dann diese Methode wiederholt auf, um die Werte für einzelne Hypothesen einzuholen, bevor es diese in der Verteilung Distrib setzt. Das Framework bietet weiterhin eine Methode print(), mit der es die Werte aller aktualisierten Wahrscheinlichkeiten für die jeweilige Hypothese ausgibt.

Listing 3: HypoTest.pm

    01 use Moose;
    02 use Distrib;
    03 
    04 has 'distrib' => ( 
    05     is => 'rw', 
    06     isa => 'Distrib', 
    07     default => sub { Distrib->new() } );
    08 
    09 sub hypo_add {
    10   my( $self, $hypo ) = @_;
    11 
    12   $self->distrib()->set( $hypo, 1 );
    13 }
    14 
    15 sub update {
    16   my( $self, $data ) = @_;
    17 
    18   my $hypos = $self->distrib()->values();
    19 
    20   for my $hypo ( keys %$hypos ) {
    21     $self->distrib()->mult( $hypo,
    22       $self->likelihood( $data, $hypo ) 
    23     );
    24   }
    25 
    26   $self->distrib()->normalize();
    27 }
    28 
    29 sub print {
    30   my( $self ) = @_;
    31 
    32   my $values = $self->distrib()->values();
    33 
    34   for my $hypo ( keys %$values ) {
    35     print "$hypo: $values->{ $hypo }\n";
    36   }
    37 }
    38 
    39 sub likelihood {
    40     die "Subclass needs to override this";
    41 }
    42 
    43 1;

So nimmt likelihood() in Listing 4 in $data vom Hauptprogramm den Buchstaben "R" entgegen (um eine gezogene Karte mit roter Vorderseite als Zusatzbedingung anzumelden) und rechnet dann basierend auf der ebenfalls hereingereichten Hypothese ("RR", "RB", "BB") aus, wie wahrscheinlich es ist, dass der Proband auf eine rote Vorderseite blickt: Eine glatte 1 (also 100%) für "RR", 0.5 für "RB" und 0 für "BB". Die Funktion setzt dazu einen regulären Ausdruck ab, der die Anzahl der "R"s in der Hypothese zählt, und dividiert das Ergebnis anschließend durch zwei als Fließkommawert, damit Perl keine Integerdivision vornimmt und den Divisionsrest unterschlägt. Am Ende gibt hypotest die Wahrscheinlichkeiten für alle Hypothesen in der Verteilung mittels der Methode print() aus dem Modul HypoTest an und berichtet zurecht, dass die rot-rote Karte in 2/3 aller Fälle erscheint und damit der Proband, der auf eine rote Vorderseite blickt, nach dem Umdrehen der Karte auch eine rote Rückseite sieht:

    $ ./hypotest 
    RB: 0.333333333333333
    RR: 0.666666666666667
    BB: 0

Listing 4: hypotest

    01 #!/usr/local/bin/perl -w
    02 use strict;
    03 package CardHypoTest;
    04 use base qw( HypoTest );
    05 
    06 sub likelihood {
    07   my( $self, $data, $hypo ) = @_;
    08 
    09     # count the number or Rs in the string
    10   my @matches = ( $hypo =~ /($data)/g );
    11   return @matches / 2.0;
    12 }
    13 
    14 package main;
    15 my $suite = CardHypoTest->new();
    16 
    17 $suite->hypo_add( "RR" );
    18 $suite->hypo_add( "BB" );
    19 $suite->hypo_add( "RB" );
    20 
    21 $suite->update( "R" );
    22 $suite->print();

Mit einfachen Skripts lassen sich so mathematisch fundierte Folgerungen aufstellen, und auch wer schon beim Anblick einer Formel in Ohnmacht fällt, kommt so auf seine Kosten und lernt, wie weitere Beispiele in [2] oder [3] (zwar in Python und nicht in Perl) zeigen, dass statistische Alltagsprobleme oft erstaunliche Lösungen haben, die auf den ersten Blick gar nicht intuitiv zu erfassen sind.

Infos

[1]

Listings zu diesem Artikel: ftp://www.linux-magazin.de/pub/listings/magazin/2014/07/Perl

[2]

"Think Bayes", Allen B. Downey, O'Reilly, 2013

[3]

"Probabilistic-Programming-and-Bayesian-Methods-for-Hackers", Cameron Davidson-Pilon, https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers

[4]

"Ziegenproblem" - Statistische Probleme -- geknackt mit Perl, Michael Schilli, Linux-Magazin 09/1999

[5]

Satz von Bayes, http://de.wikipedia.org/wiki/Satz_von_Bayes

[6]

Ziegenproblem, http://de.wikipedia.org/wiki/Ziegenproblem

Michael Schilli

arbeitet als Software-Engineer bei Yahoo in Sunnyvale, Kalifornien. In seiner seit 1997 laufenden Kolumne forscht er jeden Monat nach praktischen Anwendungen der Skriptsprache Perl. Unter mschilli@perlmeister.com beantwortet er gerne Ihre Fragen.