Mathematischer Programmierkurs in Matlab

3. Monte-Carlo-Simulation und der Mex-Compiler

Zwar hängen Monte-Carlo-Simulation und die Handhabung des Matlab-mex-Compilers in keiner Weise voneinander ab -- aber da die Verwendung des mex gerade für zeitaufwändige Simulationen oft vorteilhaft ist soll es hier mit behandelt werden.

Zu detaillierten Grundlagen zu Monte-Carlo-Simulation sei auf die Vorlesung "Stochastische Simulation" verwiesen, hier werden wir uns vor allem um Anwendungen kümmern.

Grundlagen

Informieren Sie sich, welche Zufallszahlen-Generatoren im Standard-Matlab und in den Toolboxen zu finden sind.
Erzeugen Sie mit diesen Generatoren Zufallszahlen und stellen Sie diese in einem Histogramm dar.
Wie verhält sich der Zufallszahlengenerator beim Neustart von Matlab ?
Welche Auswirkung hat das Leeren des Workspace mit clear ?
Wie kann man den ZZ-Generator initialisieren so dass bei jedem neuen Durchlauf immer andere Folgen entstehen ?
Erzeugen Sie normalverteilte Zufallszahlen ohne Verwendung einer fertigen Matlab-Funktion. Stellen Sie diese in einem Histogramm zusammen mit der (analytischen) Dichtefunktion dar.

Dartpfeile auf Pi ;-)

Zum Start soll eine Monte-Carlo-Simulation zur Schätzung von Pi implementiert werden. Im Algorithmus werden dazu zwei gleichverteilte Zufallszahlen x,y aus [0,1] generiert und der Versuch zählt als Treffer, wenn gilt \sqrt{ 1-x^2 } > y.

Entwickeln Sie ein Simulationsprogramm, welches alle K (z.B K=100) Simulationen eine Statuszeile mit der aktuellen Schätzung von Pi ausgibt.

Ratsimu -- Die Ratte im Labyrinth

Dies ist ein Beispiel einer klassischen Markov-Kette mit 2 Absorptionszuständen.

Eine Ratte startet in der Kammer Nummer 4 und entscheidet sich völlig zufällig für eine der möglichen Richtungen. In der grünen Kammer Nummer 5 steht eine Futterkiste und die Ratte bleibt für alle Zeiten dort. In der roten Kammer Nummer 6 hingegen lauert eine Schlange und die Ratte wird gefressen.

Nun die Programmieraufgabe:

Ermitteln Sie mit einer Monte-Carlo-Simulation die Wahrscheinlichkeit, dass die Ratte die Futterkiste erreicht.
Wie viele Simulationen erscheinen für eine stabile Schätzung sinnvoll ?
Ermitteln Sie die Laufzeit für 100,1000,1 Million, 10 Millionen ... Simulationen.
Welchen Schluss ziehen Sie beim Vergleich mit den Laufzeiten eines analogen C-Programmes, welches für 100 Million Simus, ca. 10 Sekunden braucht ?

Der mex - Compiler und ein Hello World

Der sogenannte mex-Compiler ist eine Möglichkeit Funktionen für Matlab z.B. in C oder Fortran zu schreiben, diese zu compilieren und damit wirklich schnellen Code zu erreichen. Wie üblich soll zuerst ein Hello World - Programm als mex-Erweiterung von Matlab geschrieben werden. Speichern Sie dazu den Quelltext hworld.c


/**************************************************************
 *
 *  First test, compile with:
 *    mex hworld.c
 *  inside matlab, then call 
 *    hworld
 *
 *  R. Unger 10/2010
 *
 **************************************************************/ 

#include "mex.h"

void mexFunction( int nlhs, mxArray *plhs[] ,
                   int nrhs, const mxArray *prhs[] )
{
    printf("Hello World from mex - c \n");
    printf("  nlhs = %i \n",nlhs);
    printf("  nrhs = %i \n",nrhs);


    return;

}

in Ihrem Arbeitsverzeichnis ab und compilieren Sie es (im Matlab-Command-Window) mit

>> mex hworld.c

und rufen Sie es anschliessend mit verschiedenen Varianten

>>hworld
>>hworld(1.23)
>>x=hworld
>>x=hworld(1.23)
>>[x,y] = hworld();

auf. Beobachten Sie die Ausgaben und verstehen Sie alle Variablen und Parameter.

mex mymin

Analog zur ersten Übung soll noch einmal die Minimumssuche in einem Vektor, diesmal über eine mex-Funktion implementiert werden.

Dazu ist der Quelltext schon fertig: mymin.c. Speichern Sie den Quelltext in Ihrem Arbeitsverzeichnis, erstellen Sie mit dem mex-Compiler das Binary und testen Sie die Laufzeiten im Vergleich mit den beiden Varianten aus der ersten Übung.

mexrat

Bauen Sie nun obiges hworld.c aus, um die Ratte-im-Labyrinth-Simulation komplett in C auszulagern. Die Funktion soll mexrat.c heissen, einziger Eingangsparameter soll die Anzahl der Simulationen sein, einziger Ausgabeparameter soll die Wahrscheinlichkeit des Erreichens der Futterkiste ein. Der Aufruf erfolgt also als:

>> p = mexrat(1000);

Messen Sie auch hier wieder die Laufzeiten für 10,100 ... 10 Mio Simus und vergleichen Sie mit der reinen Matlab Variante.

Hinweis: Nicht zu sehr ärgern wenn das komplette Matlab beim Entwickeln ein paar mal abstürzt, manchmal hat die Verwendung des mex-Compilers ein wenig Gemeinsamkeiten mit Russischem Roulette. ;-)

Solaris

solaris.m