/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ /* [ Created with wxMaxima version 0.8.4 ] */ /* [wxMaxima: input start ] */ /* Die Maxima-Pakete "distrib" und "descriptive" enthalten statistische Funktionen */ /* und werden fuer diese Simulation gebraucht */ load (distrib) $ load (descriptive); /* Die Simulation benoetigt normalverteilte Zufallszahlen */ /* Diese effizient aus rechteckvertteilten 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 Laenge von 480 Elementen erzeugt und zunaechst mit Nullen gefuellt */ t_carlo:makelist(0.0,x,-240,240)$ /* Sie soll spaeter die Treffer der MonteCarlo Simulation akkumulieren */ /* Die Liste "tru_t" wird in einer Laenge von 480 Elementen erzeugt und */ /* mit der "richtigen" Student-t-verteilung aus der distrib-Package befuellt */ /* 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 traegt die aequidistanten x-Koordinaten von -3 bis +3, die fuer 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 fuer den Vergleich im Plot benoetigt, schliesslich sollen unabhaengig von der Zahl der Schuesse */ /* die MC-Simulation im Pegel mit der echten Student-t Verteilung uebereinstimmen */ /* Spaeter(drittletzte Zeile) wird fuer 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 zurueck */ /* 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 Schaetzung des Fehlers fuer den Mittelwert der Grundgesamtheit auf Basis nur */ /* einer Stichprobe ohne Kenntnis der wahren Standardabweichung */ /* 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 fuer den Stichprobenumfang 3 (2 Freiheitsgrade)*/ /* incrementiere das vom MonteCarlo-Schuss getroffene Tabellenelement von "t_carlo" */ /* ab und zu normiere die t-carlo-Liste und ploitte 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"]))$, ); /* [wxMaxima: input end ] */ /* Maxima can't load/batch files which end with a comment! */ "Created with wxMaxima"$