zurück zur Homepage     zurück zur Kurshomepage Statistik

Monte Carlo Simulation zur Generierung der Student-t Verteilung mit Maxima

MC-Simulation-Grafikoutputstud_t_explicit

Die Herleitung der expliziten Formel für die Student-t Verteilung ist schwierig, würde eine normale Statistikvorlesung sprengen und wird deshalb in aller Regel nicht präsentiert.
Man kann die Verteilung und deren Anwendung aber trotzdem verstehen, wenn man eine Monte Carlo Simulation versteht, die sie generiert.
Sie können das auf Ihrem Computer durchführen und der Verfeinerung bei größer werdender "Schusszahl" der Simulation zuschauen.

Sie brauchen dafür lediglich die kostenfreie und sowieso überaus nützliche Software(Open Source/Public License) Maxima für symbolische (und auch numerische) Berechnungen:

1. Ich empfehle wxMaxima, das Sie von hier kostenfrei für alle Betriebssysteme herunterladen und dann installieren können.
(Hinweis: Linux-user finden das Paket in der Regel auch in Ihren Repositories.)

2. Folgende kleine wxMaxima-Zelle "maxima_student_t.wxm" bekommen Sie von mir. Laden Sie sie herunter und legen Sie sie an einem Ort Ihrer Wahl auf dem Computer ab. Mit dieser Zelle erzeugen Sie eine Monte Carlo Simulation zur Student-t Verteilung zu 2 Freiheitsgraden.
(Hinweis: Wegen verschiedener Zeichencodierung in verschiedenen Umgebungen können die Umlaute in Ihrem Editor eventuell etwas wunderlich dargestellt werden - das ist aber kein Problem, da das nur in Kommentarzeilen und auch nur außerhalb von wxMaxima auftritt.)

3. Öffnen Sie einfach die Zelle "maxima_student_t.wxm" mit wxMaxima.
(Hinweis: Sicherheitsbewusste Besucher können auch den unten auf dieser Seite wiedergegeben Inhalt selbst abtippen oder per Copy and Paste in eine leere wxMaxima-Zelle einfügen.)
Starten Sie die Berechnung der Maxima-Zelle(Gehen Sie mit dem Cursor in das Programm und drücken Sie "shift return").
Die Berechnung kann durch Schließen des Maxima-Haupt-Fensters vorzeitig abgebrochen werden.

Hinweis:
Unter Windows (XP) werden die Iterationsblöcke jeweils durch Schließen des Plotfensters getriggert.
Unter Linux (Ubuntu Lucid) erfolgen die Iterationsblöcke mit Grafikaktualisierung automatisch und werden bei gleicher Hardware etwa in doppelter Geschwindigkeit ausgeführt.

Kleine Knobeleien zum Schluss gefällig ?
A: Welche (zwei !) Zahlen müssen geändert werden, um die Verteilung für 3,4,5,..n Freiheitsgrade zu generieren ?
B: Wie könnte man die Maxima-Zelle abwandeln, um einen größeren Wertebereich als -3  - +3 zu erfassen ?
C: Wie könnte man die Maxima-Zelle abwandeln(vereinfachen), um die Standardnormalverteilung zu generieren ?
Der Inhalt der Maxima-Zelle ist hier wiedergegeben:

/* [ Created with wxMaxima version 0.8.4 ] */

/* Die Maxima-Pakete "distrib" und "descriptive" enthalten statistische Funktionen */
/* und werden für diese Simulation gebraucht */
load (distrib) $
load (descriptive);

/* Die Simulation benötigt normalverteilte Zufallszahlen */
/* Diese effizient aus rechteckverteilten Zahlen zu erzeugen ist nicht trivial */
/* Unter Maxima scheint das Verfahren nach Marsaglia am besten zu laufen*/
/* die jetzt definierte Funktion "gaussdm" erzeugt bei Aufruf eine standardnormalvt. Zufallszahl */
gaussdm():= block(s:2,for im:1 thru 900000 while s >= 1.0
do block(x:random(2.0)-1.0, y:random(2.0)-1.0, s:(x^2 +y^2)),ret:y*sqrt(-2.0*log(s)/s)/*,if x>y then ret:-ret ,ret*/);

/* Die Liste "t_carlo" wird in einer Länge von 480 Elementen erzeugt und zunächst mit Nullen gefüllt */
t_carlo:makelist(0.0,x,-240,240)$ /* Sie soll später die Treffer der MonteCarlo Simulation akkumulieren */

/* Die Liste "tru_t" wird in einer Länge von 480 Elementen erzeugt und */
/* mit der "richtigen" Student-t-verteilung aus der distrib-Package befüllt */
/* Die Liste gibt die Werte von -3 bis +3 wieder ->(240/80 = 3) */
tru_t:makelist(pdf_student_t(1/160+x/80,2),x,-240,240)$

/* Die Liste "xval" wird erzeugt. Sie trägt die äquidistanten x-Koordinaten von -3 bis +3, die für den Plot */
/* jeweils der echten Student-t und der MonteCarlo-Simulation zugeordnet werden */
xval:makelist(x/80,x,-240,240)$

/* Die Funktion normlist wird definiert, die Liste "lst" auf eine bestimmte "norm" normiert */
/* Wird für den Vergleich im Plot benötigt, schliesslich soll unabhängig von der Zahl der Schüsse */
/* die MC-Simulation im Pegel mit der echten Student-t Verteilung übereinstimmen */
/* Später(drittletzte Zeile) wird für das Intervall -3 - +3 auf die Norm 0.904 normiert, der Wert kann student-t Tabellen entnommen werden */
normlist(lst,norm) :=(block(a:0.0, in:0), for in: 1 step 1 thru length(lst) do (a:a +lst[in]),
for in: 1 step 1 thru length(lst) do (lst[in]: norm*lst[in]/a));

/* Diese Funktion zieht "anz" standardnormalverteilte Zahlen, und gibt eine liste mit zwei Elementen zurück */
/* ertes Listenelement der Returnliste ist der Mittelwert, das zweite die empirische Standardabw. */
ggs_m(anz):= block(x:0, a:0, ims:0, mitt:0.0, loclist:makelist(0.0,x,1,anz),
retlist:makelist(0.0,x,1,2), for ims:1 step 1 thru anz do(loclist[ims]: gaussdm()),
for ims:1 step 1 thru anz do(a: a+loclist[ims]), retlist[1]: a/anz,
retlist[2]: std1(loclist),retlist );

/* Diese Funktion ermittelt den Abstand des wahren Mittelwerts der Standardnormalvt. */
/* IN EINHEITEN DER EMPIRISCHEN Standardabweichung(nicht der wahren Standardabweichung!) */
/* vom Mittelwert der Stichprobe. */
/* Dies bildet die Schätzung des Mittelwertfehlers gegenüber der Grundgesamtheit aufgrund nur einer Stichprobe ab */
/* und soll damit Student-t verteilt sein. */
gme_m(anz):= block(l1:makelist(0,x,1,2), l1:ggs_m(anz), a:-l1[1]*sqrt(3)/l1[2]);

/* Mach das viele, viele, viele Male für den Stichprobenumfang 3 (2 Freiheitsgrade)*/
/* incrementiere das vom MonteCarlo-Schuss getroffene Tabellenelement von "t_carlo" */
/* bilde ab und zu eine normierte Version der t-carlo-Liste (n_tcarlo)und plotte sie zusammen mit der echten Student-t-Verteilung */
for ii:1 thru 400 do
block( for iz: 1 step 1 thru (1000 + 6*ii) do
( a: gme_m(3), indx: 240 + floor(0.5+80*a), if ((indx > 0) and (indx <= 481)) then t_carlo[indx]:t_carlo[indx] +1.0 ),
n_tcarlo: copylist(t_carlo),
normlist(n_tcarlo,80*0.904),
plot2d([[discrete,xval,n_tcarlo], [discrete,xval,tru_t]],[legend,"MC_Sim","student_t"]))$,
);