Sobol-Sequenzen: Dem Zufall auf die Sprünge helfen

Eines der typischen Probleme, die sich bei geographischen Informationssystemen stellen, ist die Frage, wie groß die Fläche ist, die sowohl in Gebiet A als auch in Gebiet B liegt, also beispielsweise in einer Gemeinde namens Zufallshausen, die ihrerseits im Vorwahlbereich 0815 liegt.

Mit solchen Verschneidungen kann man versuchen, Daten, die auf Gemeindeebene vorliegen, auf die Vorwahlbereiche zu übertragen. Das Problem lässt sich zwar exakt durchrechnen, aber die Verschneidung zweier Gebiete, die jeweils sowohl Enklaven als auch Exklaven besitzen können, ist eine ziemlich anspruchsvolle und aufwändige Aufgabe.

Die hier vorgestellte Methode hat den Vorteil, dass sie in wenigen Minuten an die unterschiedlichsten Problemstellungen angemaßt werden kann, auch wenn sie nicht die von Geodäten gewünschte Genauigkeit liefert. Doch dafür funktioniert sie selbst in hoch dimensionalen Räumen, wo andere Methoden scheitern.

Abbildung 1: Vergleich der Monte-Carlo-Integration mit Zufallszahlen und Sobol- Sequenzen in Abhängigkeit von der Anzahl der Iterationsschritte. Auch wenn es auf den ersten Blick nicht so scheint: Mit genügend vielen Iterationen nähert sich auch die Version mit Zufallszahlen dem exakten Ergebnis an.

Das Testproblem

Als Übungsaufgabe werden wir hier die Fläche berechnen, die sowohl im Einheitskreis als auch im Quadrat |x| 0,815, |y| 0,815) liegt. Eine früher übliche und durchaus legitime Lösung besteht darin, das Schnittgebiet auf Karton zu zeichnen, auszuschneiden und zu wiegen. Ein Vergleich mit dem Gewicht des Einheitsquadrats liefert das Ergebnis. Für Digitalrechner ist das aber noch zu kompliziert. Deshalb hängen wir unsere Zeichnung an die Wand und werfen mit Dartpfeilen darauf.

Nach einer stattlichen Anzahl von Würfen zählen wir die Treffer in unserem Schnittgebiet und die Treffer im umschließenden Quadrat |x|, |y| 1. Der Quotient liefert uns ein Näherung der Lösung. Die Umsetzung in C++ ist in Listing 1 zu finden, die Sourcen auf [1]. Für eine Diskussion der Problematik numerisch erzeugter Zufallszahlen verweise ich auf [2].

Zunächst definieren wir in Zeile 4 bis 8 eine Zufallszahlen-Funktion, die eine zufällige Zahl zwischen Null und Eins liefert unter Verwendung der C++ Funktion rand(). Zeilen 10 bis 14 definieren eine alternative Version mit der oft verwendeten drand48() ( 48 steht für die Anzahl der intern verwendeten Bits).

Im Hauptprogramm bestimmen wir zunächst die Anzahl der gewünschten Iterationen N. Unsere Zufallszahl skalieren wir in Zeile 31 bis 32 auf das umschließende Rechteck. Schließlich testen wir in Zeile 36, ob unser Zufallspunkt sowohl in A als auch in B liegt. Falls ja, erhöhen wir unseren Score. Die gesuchte Fläche ergibt sich nun einfach aus dem Quotienten der Treffer und der Versuche multipliziert mit der Fläche des umschließenden Rechtecks. Für komplexere Flächen braucht man nur den Hit-Test durch entsprechende Punkt-in-Polygon-Funktionen zu ersetzen.

Abbildung 2: Verteilung der ersten 10000 x/y-Zufallszahlen von rand().

1:    #include <iostream>
 2:    #include <stdlib.h>
 3:
 4:    void Rand_1( double& x, double& y)
 5:    {
 6:            x = double(rand()) / RAND_MAX;
 7:            y = double(rand()) / RAND_MAX;
 8:    }
 9:
10:    void Rand_2( double& x, double& y)
11:    {
12:            x = drand48();
13:            y = drand48();
14:    }
15:
16:    int main(int argc, char* argv[])
17:    {
18:            int N = 25;
19:            int nHit = 0;
20:            double x , y, r;
21:
22:            if( argc >1 )
23:                    N = atoi( argv[1] );
24:
25:            for( int i = 0; i < N; ++i)
26:            {
27:                    //Bestimmung der Zufallszahl
28:                    Rand_1( x, y );
29:
30:                    // Einpassen ins umschließende Rechteck
31:                    x = 2.0*x - 1.0;
32:                    y = 2.0*y - 1.0;
33:
34:                    // Hit Test
35:                    r = x*x + y*y;
36:                    if( (r <= 1) && ( fabs(x) <= 0.815 ) && ( fabs(y) <= 0.815 ) )
37:                            ++nHit;
38:
39:                    // Ausgabe
40:                    //cout << i << " " << x << `<>t' << y << `<>t'
41:                    //      << 4.0 * nHit / double(i+1) << endl;
42:            }
43:
44:            cout.precision(12);
45:            cout << N << `<>t' << 4.0 * nHit / double(N) << endl;
46:    }

Monte Carlo mit Sobol-Sequenzen

Leider hat das Verfahren einen Nachteil: Die Genauigkeit ist typischerweise proportional zu . Die gute Nachricht: Es geht besser, statt Zufallszahlen, deren Erzeugung sowieso prinzipielle Probleme aufwirft, kann man deterministische Zahlenfolgen verwenden, die das Gebiet möglichst gut abtasten, etwa Fraktale oder wie hier Sobol-Sequenzen.

Abbildung 1 vergleicht das Ergebnis der Monte-Carlo-Simulation mit Zufallszahlen mit dem bei Sobol-Sequenzen. Das Ergebnis spricht für sich, da die Erzeugung der Sobol-Sequenzen schneller ist als die der (Pseudo-)Zufallszahlen - die ersten zehn Millionen Sobol-Paare benötigen auf meinem AMD 800 weniger als zwei Sekunden.

Das Programm ist identisch mit dem vorherigen, bis auf die Verwendung eines Sobol-Zahlen-Generators statt der Zufallszahlen. Dieser Generator stammt aus der Gnu Scientific Library v. 0.8, alle Sourcen unterliegen also automatisch der GPL [1].

Im Gegensatz zu Pseudo-Zufallszahlen sind Sobol-Sequenzen deterministische Zahlenfolgen, die gleichmäßiger als Zufallszahlen verteilt sind. Der große Vorteil bei den Sobol-Sequenzen ist, dass sie sich sozusagen gegenseitig meiden: Die Wahrscheinlichkeit, in der Nähe einer Sobol-Zahl eine zweite zu finden, ist kleiner als für Zufallszahlen.

Das sieht man im Vergleich von Abbildung 2 mit Abbildung 3, wo jeweils die ersten 10000 Zahlenpunkte aufgetragen sind: Sobol-Sequenzen bedecken die Fläche gleichmäßiger als Zufallszahlen. Für eine weiter gehende Einführung verweise ich auf die Originalliteratur [3] und die berühmt berüchtigten Numerical Recipes [4].

Abbildung 3: Verteilung der ersten 10000 x/y-Zufallszahlen von Rand_Sobol().