Notebook zu: Mathematik für Informatiker
Band 2

©2001-2018 Gerald Teschl (http://www.mat.univie.ac.at/~gerald/) und Susanne Teschl  (http://staff.technikum-wien.at/~teschl/)

Dieses Notebook ist als Ergänzung zu unserem Buch Mathematik für Informatiker, Band 2, 3. Aufl., Springer 2014, gedacht. Dort finden Sie auch weitere Erklärungen.

Elementare Funktionen

Nullstellen

Der Solve -Befehl in Mathematica versucht, eine Gleichung nach der angegebenen Variablen aufzulösen. Gleichungen werden mit einem doppelten Gleichheitszeichen eingegeben :

Solve[(x^3)-7(x^2)+7x-1==0,x]

Mathematica hat also drei Nullstellen gefunden.

Faktorisierung

Polynome können mit dem Befehl

Factor[(x^3)-1]

faktorisiert und mit dem Befehl

Expand[%]

ausmultipliziert werden. Eine Faktorisierung erfolgt nur, wenn die Nullstellen rationale Zahlen sind.

Polynomdivison

Den Quotienten s(x) der Polynomdivision p(x) = s(x) q(x) + r(x) können wir mit PolynomialQuotient[p (x), q (x), x]

PolynomialQuotient[3(x^4)+(x^3)-2x,(x^2)+1,x]

und den Rest r(x) mit PolynomialRemainder[p(x), q(x), x]

PolynomialRemainder[3(x^4)+(x^3)-2x,(x^2)+1,x]

berechnen.

Interpolationspolynome

Interpolationspolynome können mittels InterpolatingPolynomial[{{(Subscript[x, 0]), (Subscript[y, 0])}, . . .}, x] berechnet werden, wobei ((Subscript[x, 0]), (Subscript[y, 0])), ... die gewünschten Stützpunkte sind. Um uns die Eingabe der Stützpunkte zu ersparen, können wir leicht mit dem Table-Befehl eine Liste gleichmäßig verteilter Stützpunkte erzeugen. Zum Beispiel :

f[x_]=(1/(1+25 (x^2)));

Stuetzpunkte[n_]:=Table[{(j/n),f[((j )/n)]},{j,-n,n}]

P[x_]=InterpolatingPolynomial[Stuetzpunkte[6],x]//Expand

Stuetzpunkte[n] erzeugt eine Liste von 2 n + 1 Stützpunkten, die gleichmäßig im Intervall[−1, 1] verteilt sind. Damit wurden hier 13 Stützpunkte erzeugt und dann
das Interpolationspolynom P(x) von f(x) = (1/(1+25 (x^2))) berechnet.

Plot[{f[x],P[x]},{x,-1,1},PlotRange->All,PlotStyle->{Black,GrayLevel[0.5`]}]

Clear[f,P];

Graphische Darstellung von Funktionen

Erinnern Sie sich an den Plot-Befehl, mit dem Funktionen gezeichnet werden können. Es ist dabei sogar möglich mehrere Funktionen gleichzeitig zu zeichnen, indem man statt einer Funktion eine Liste von Funktionen angibt :

Plot[Sin[x],{x,0,2 π},Ticks->{{0,(π/2),π,((3 π)/2),2 π},Automatic}]

Mit der Option Ticks wurde hier festgelegt, dass die x - Achse an den Positionen 0, (π/2),  π, etc.  zu beschriften ist (die Beschriftung der y - Achse erfolgt hier automatisch).

Es ist dabei sogar möglich mehrere Funktionen gleichzeitig zu zeichnen, indem man statt einer Funktion eine Liste von Funktionen angibt :

Plot[{Sin[x],Sin[2 x]},{x,0,2 π},Ticks->{{0,π,2 π},Automatic},PlotStyle->{{},Dashing[{0.05`,0.02`}]}]

Mit der Option PlotStyle wurde für die zweite Funktion durch Dashing[{0.05, 0.02}] eine strichlierte Linie gewählt. Als PlotStyle kann durch Angabe von RGBColor[r, g, b] auch ein Farbwert gewählt werden (z.B. RGBColor[1, 0, 0] für rot)

Plot[{Sin[x],Sin[2 x]},{x,0,2 π},Ticks->{{0,π,2 π},Automatic},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}}]

Bogenmaß und Gradmaß

Winkelwerte werden von Mathematica defaultmäßig im Bogenmaß interpretiert:

Sin[(π/2)]

Haben Sie den Winkel im Gradmaß gegeben, so muss er mit dem Faktor (π/180)  ins Bogenmaß umgerechnet werden. Dafür gibt es die eingebaute Konstante Degree, die auch als ◦ über die Palette eingegeben werden kann :

Sin[90 Degree]

Funktionen mit komplexen Argumenten, Polardarstellung

Natürlich können wir in Mathematica die trigonometrischen (bzw. alle anderen) Funktionen auch mit komplexem Argument aufrufen (Eingabe der imaginären Einheit mit dem Großbuchstaben “I” oder über die Palette).

Sin[I]

Oder:

Exp[2 + I (π/3)]//ComplexExpand

Hier mussten wir noch mit ComplexExpand nachhelfen, damit das Ergebnis auch in Real - und Imaginärteil aufgespalten wird.

Der Betrag einer komplexen Zahl kann mit Abs[z] und der Winkel mit Arg[z] berechnet werden. Der Winkel wird im Bogenmaß ausgegeben. Die Polarkoordinaten von 1 + i (Sqrt[3]) sind zum Beispiel

{Abs[1+I (Sqrt[3])],Arg[1+I (Sqrt[3])]}

Die n-te Wurzel einer komplexen Zahl z wird in der Form (z^n) oder über die Palette eingegeben :

(Sqrt[1 + I(Sqrt[3])])//ComplexExpand

Oder:

((8Exp[I π/2])^(1/3))//ComplexExpand

Den komplexen Logarithmus erhält man mit

Log[1 + I(Sqrt[3])]//ComplexExpand

und analog komplexe Potenzen

{(3^(-I)),(((1+I))^I),(((1-I))^(2+3I))}//ComplexExpand//Simplify

Differentialrechnung  I

Grenzwert

Grenzwerte erhalten wir wie bei Folgen mit dem Befehl Limit:

Limit[(((x^2)-1)/(x-1)),x->1]

Ableitungen

Ableitungen werden mit dem Befehl D berechnet :

D[2 (x^3)-4x+1,x]

Die n-te Ableitung nach x bekommen wir mit D[f(x), {x, n}] :

D[2 (x^3)-4x+1,{x,2}]

Eine Ableitung kann in Mathematica auch einfach als f'[x] anstelle von D[f[x], x] eingegeben werden

f[x_]:=2 (x^3)-4x+1; f'[x]

Clear[f]

Splines

Mathematica stellt leider keinen Befehl zur Berechnung von Splinefunktionen zur Verfügung, man muss also selbst zur Tat schreiten und den Algorithmus aus dem Buch implementieren. (Im Paket NumericalMath'SplineFit' gibt es allerdings eine Funktion SplineFit, die aus gegebenen Stützpunkten eine Splinekurve berechnet.)

NaturalSpline[xy_,t_]:=Module[{x,y,z,n,u,v,d,k,Fh,F,spl},
{x,y}=Transpose[xy];
n=Length[x];u=Table[0,{j,n}];v=Table[0,{j,n}];
d=Table[x[[j+1]]-x[[j]],{j,n-1}];
k=Table[((y[[j+1]]-y[[j]])/(d[[j]])),{j,n-1}];
u[[1]]=0;u[[2]]=3 k[[1]];v[[1]]=1;v[[2]]=-2;
Fh[j_,z_]:=((-1)/(d[[j-2]]))(2(d[[j-1]]+d[[j-2]])z[[j-1]]+ d[[j-1]] z[[j-2]]);
F[j_,z_]:=Fh[j,z]+(3/(d[[j-2]]))(k[[j-2]]d[[j-1]]+k[[j-1]]d[[j-2]]);
Do[u[[j]]=F[j,u];v[[j]]=Fh[j,v];,{j,3,n}];
z=u+ ((3k[[n-1]]-u[[n-1]]-2u[[n]])/(2v[[n]]+v[[n-1]])) v;
spl=Flatten[Table[{t<x[[j+1]],Expand[y[[j]]+k[[j]](t-x[[j]])+((z[[j]]-k[[j]])/((d[[j]])^2))(t-x[[j]])(((t-x[[j+1]]))^2)+((z[[j+1]]-k[[j]])/((d[[j]])^2))(((t-x[[j]]))^2)(t-x[[j+1]])]},{j,n-1}]] ;
spl[[-2]]=True;spl [[0]]=Which;
spl
];

Ein Vergleich mit dem Interpolationspolynom aus dem letzten Abschnitt:

f[x_]=(1/(1+25 (x^2)));

Stuetzpunkte[n_]:=Table[{(j/n),f[(j/n)]},{j,-n,n}];
XY=Stuetzpunkte[6];
S[x_]=NaturalSpline[XY,x];
P[x_]=InterpolatingPolynomial[XY,x];
dots=ListPlot[XY,PlotStyle->{PointSize[0.02`]}];
gr1=Plot[f[x],{x,-1,1},PlotRange->{-1,1}];
gr2=Plot[S[x],{x,-1,1},PlotRange->{-1,1}];
gr3=Plot[P[x],{x,-1,1},PlotRange->{-1,1}];
Show[gr1,dots,AxesLabel->{"x","f[x]"}]
Show[gr2,dots,AxesLabel->{"x","Spline"}]
Show[gr3,dots,AxesLabel->{"x","Interpolationspolynom"}]

Clear[f,S,P,XY]

Differentialrechnung II

Taylorpolynome

Taylorpolynome werden mit dem Befehl Series berechnet. Anzugeben ist dabei das Argument der Funktion (hier x), der Entwicklungspunkt (hier 0), sowie der gewünschte Grad des Polynoms (hier 3) :

Series[Sin[x],{x,0,3}]

Den “O-Term” können wir mit

Normal[%]

entfernen.

Gauß'sche Glockenkurve

Um die Extrema und Wendepunkte der Gauß'schen Glockenkurve zu berechnen, definieren wir zunächst

f[x_]:=(1/((Sqrt[2π])σ))(E^(-((((x-μ))^2)/(2(σ^2)))));

und bestimmen dann die Nullstellen der Ableitung:

f'[x]

Es gibt also eine mögliche Extremstelle, x = μ. Da die zweite Ableitung an dieser Stelle

Solve[f'[x]==0,x]

f''[μ]

negativ ist, liegt hier ein Maximum  vor. Für die Wendepunkte ermitteln wir die Nullstellen der zweiten Ableitung,

Solve[f''[x]==0,x]

und werten die dritte Ableitung an diesen Stellen aus:

{f'''[μ-σ],f'''[μ+σ]}

Da das Ergebnis jeweils ungleich null ist, liegen hier tatsächlich Wendepunkte vor.

Und so sieht die Glockenkurve aus:

μ=0;σ=1;Plot[f[x],{x,-4,4}]

Sekantenverfahren

Wir definieren die Funktion, deren Nullstellen wir numerisch mit dem Sekantenverfahren bestimmen möchten:

f[x_]:=x^3+2x^2+10x-20;

{f[1],f[2]}

Dann geben wir die Formel für die rekursive Berechnung der Näherungswerte ein:

x[n_]:=x[n-1]-f[x[n-1]]((x[n-1]-x[n-2])/(f[x[n-1]]-f[x[n-2]]))

Nun geben wir die Startwerte vor und lassen uns eine Liste der nächsten Näherungswerte ausgeben:

x[0]=1.;x[1]=2.;{x[2],x[3],x[4],x[5]}

Clear[f,x]

Newtonverfahren

Wir geben die Funktion ein, deren Nullstellen wir bestimmen möchten, und geben die rekursive Formel für das Newton-Verfahren ein :

f[x_]:=x^2-2;
x[n_]:=x[n-1]-((f[x[n-1]])/(f'[x[n-1]]))

Mit dem Startwert (Subscript[x, 0]) = 1 lassen wir uns eine Liste der nächsten fünf Näherungswerte ausgeben :

x[0]=1.;{x[1],x[2],x[3],x[4],x[5]}

Natürlich können Sie auch einfach den Befehl FindRoot[f[x], {x, x0}] verwenden, um mit dem Startwert (Subscript[0, □]) einen Näherungswert für eine Nullstelle von f(x) zu finden.

FindRoot[f[x],{x,1}]

Clear[f,x]

Konvergenzbereich des Newton-Verfahrens

Anbei eine Funktion die mit dem Newton-Verfahren die Nullstellen von (z^3)-1  sucht und jedem komplexen Startwert eine Farbe zuordnet die der Nullstelle entspricht zu der der Startwert konvergiert.

WhichZero=Compile[{{z0,_Complex}},
Block[{z=z0,k},
For[k=0,k<200&&Abs[z^3-1]>0.5, k++,z=(1/(3 (z^2)))+((2 z)/3)];
(Round[Im[z]]+1)/2]
];
MyColor[x_]:=Which[x==0,RGBColor[1,0,0],x==1/2,RGBColor[0,1,0],x==1,RGBColor[0,0,1],True,RGBColor[1,1,1]];

DensityPlot[WhichZero[x+I y],{x,-2,2},{y,-2,2},PlotPoints->400,Mesh->False,ColorFunction->MyColor]

DensityPlot[WhichZero[x+I y],{x,-2,-1},{y,-0.5`,0.5`},PlotPoints->400,Mesh->False,ColorFunction->MyColor]

Dioden-Logik

Wir definieren zuerst (Subscript[I, D])(U)

c=39.6;
ID[U_]=10^-9(Exp[c  U] -1);

Plot[ID[U],{U,-0.5,0.5},PlotRange->All]

und unsere Werte für die Parameter (die Widerstände)

R0=10000; R1=4700;R=((R0 R1)/(R0+R1));

Nun lösen wir mit FindRoot die Gleichung für (Subscript[U, 2]) und geben eine Liste mit den Werten für (Subscript[U, 1]), (Subscript[U, 2]), (Subscript[U, 0]) aus:

U0[E1_,E2_]:={E1-E2+U2,U2,E2-U2}/.
FindRoot[U2==E2 - R (ID[E1-E2+U2]+ID[U2]),{U2,-4}]

(Der Startwert − 4 wurde gewählt, damit das Newton - Verfahren in allen vier untenstehenden Fällen konvergiert; -) Damit erhalten wir

{U0[0,0],U0[5,0],U0[0,5],U0[5,5]}  //Chop//TableForm

Der Befehl Chop ersetzt Werte, die innerhalb der Rechengenauigkeit null sind, durch 0. Das macht die Ausgabe oft leichter lesbar.

Clear[T,ID,R,R0,R1,c,U0]

Physikalische Konstanten

Mathematica kennt viele wichtige Konstanten und kann auch mit Einheiten rechnen.

Quantity["BoltzmannConstant"]

UnitConvert[%]

UnitSimplify[%]

Zum Beispiel können wir den Wert (q/(k T)) der in der Kennlinie einer Diode auftritt ausrechnen:

UnitConvert[((Quantity["ElectronCharge"])/(Quantity["BoltzmannConstant"]Quantity[293.,"Kelvin"]))]

UnitSimplify[%]

Übungen

f[x_]:=(x/(Sqrt[1+(x^2)]));
x[n_]:=x[n-1]-((f[x[n-1]])/(f'[x[n-1]]))

x[0]=0.5;{x[1],x[2],x[3],x[4]}

x[0]=1;{x[1],x[2],x[3],x[4]}

x[0]=2;{x[1],x[2],x[3],x[4]}

Clear[f,x]

Integralrechnung

Unbestimmtes Integral

Das Integralzeichen geben wir mithilfe der Palette ein.Alternativ kann auch der Befehl Integrate[f, x] verwendet werden.

∫(Sin[x]+(x^2))\[DifferentialD]x

Mathematica gibt keine Integrationskonstante aus.

Bestimmtes Integral

WiederverwendenwirdiePaletteoderdenBefehl Integrate[f, {x, xmin, xmax}]:

(∫_0^1)(2(E^x)+x)\[DifferentialD]x

Kann das Integral nicht analytisch berechnet werden, so kann ein numerischer Wert mit dem Befehl NIntegrate[f, {x, xmin, xmax}] erhalten werden.

Uneigentliches Integral

Analog berechnen wir ein uneigentliches Integral:

(∫_1^∞)(1/(t^3))\[DifferentialD]t

Gammafunktion

Die Gammafunktion ist eingebaut:

Gamma[5]

Plot[Gamma[x],{x,0,6}]

Fourierreihen

Fourierpolynom

Natürlich gibt es auch einen Befehl in Mathematica, um ein Fourierpolynom zu berechnen:

FourierTrigSeries[Abs[t],t,3,FourierParameters->{1,2π}]

Mathematica nimmt als Intervall automatisch[−π, π] an, den allgemeinen Fall [−T, T] erhält man mit der Option FourierParameters → {1, ((2 π)/T)}.

Fourierreihe

Wir definieren zunächst die Funktion

f[t_]:=If[t>0,1,0];Plot[f[t],{t,-π,π},Ticks->{{-π,-(π/2),0,(π/2),π},Automatic}]

und berechnen

a[0]=(1/π)(∫_(-π)^π)f[t]\[DifferentialD]t

Weiters gilt

a[k_]=(1/π)(∫_(-π)^π)f[t]Cos[k t]\[DifferentialD]t

Das ist aber für ganzzahliges k gleich Null, wie auch mit Mathematica leicht zu sehen ist:

Simplify[%,k∈Integers]

(Die Option k ∈ Integers legt fest, dass k ganzzahlig ist. Eingabe entweder über die Palette oder als Element[k,Integers].)

Analog berechnen wir

b[k_]=(1/π)(∫_(-π)^π)f[t]Sin[k t]\[DifferentialD]t

was sich ebenfalls noch etwas vereinfachen lässt:

Simplify[%,k∈Integers]

Berechnen wir noch explizit die ersten sechs Koeffizienten (Subscript[b, k]):

{b[1],b[2],b[3],b[4],b[5],b[6]}

Zum Schluss noch ein Vergleich zwischen der Funktion und zweier Fourierpolynome:

F[n_,t_]:=FourierTrigSeries[f[t],t,n];
Plot[Evaluate[{f[t],F[3,t],F[12,t]}],{t,-π,π},PlotStyle->{{},{Dashing[{0.04`,0.02`,0.01`,0.02`}]},{}},Ticks->{{-π,-(π/2),0,(π/2),π},Automatic}]

Clear[a,b,f,F];

Differentialrechnung in mehreren Variablen

Graphische Darstellung

Zur Veranschaulichung einer reellwertigen Funktion von zwei Variablen kann der Befehl Plot3D[f[x, y], {x, xmin, xmax}, {y, ymin, ymax}] verwendet werden.

Plot3D[2-(x^2)-(y^2),{x,-2,2},{y,-2,2}]

Plot3D[2-(x^2)-(y^2),{x,-2,2},{y,-2,2},ClippingStyle->None,PlotPoints->50,Mesh->False,PlotRange->{-1,2}]

Partielle Ableitung

Mit Mathematica können wir partielle Ableitungen wie folgt berechnen:

D[(Subscript[x, 1])((Subscript[x, 2])^3),(Subscript[x, 1]),(Subscript[x, 1])]

oder

D[(Subscript[x, 1])((Subscript[x, 2])^3),(Subscript[x, 1]),(Subscript[x, 2])]

Hesse-Matrix

Die Hesse-Matrix erhalten wir mit

Table[D[(Subscript[x, 1])( (Subscript[x, 2])-1) + ((Subscript[x, 1])^3),(Subscript[x, k]),(Subscript[x, j])],{k,2},{j,2}]//MatrixForm

Differentialgleichungen

Differentialgleichungen

Mit Mathematica können wir Differentialgleichungen mit dem Befehl DSolve lösen:

DSolve[x'[t]==2(Sqrt[x[t]]),x[t],t]

Die Parameter in der allgemeinen Lösung werden mit C[1], C[2], etc.bezeichnet.Da unsere Gleichung von der Ordnung eins ist, gibt es genau eine Konstante.

Auch Anfangsbedingungen können angegeben werden :

DSolve[{x'[t]==2(Sqrt[x[t]]),x[0]==x0},x[t],t]

Beachten Sie, dass die erste von Mathematica ausgegebene Lösung keine Lösung der gegebenen Differentialgleichung ist (sie erfüllt x′ (t) = −2 (Sqrt[x(t)])). Außerdem findet das Programm für (Subscript[x, 0]) = 0 nicht alle Lösungen, die wir im Buch gefunden haben. (Es schadet also nie, die Ergebnisse des Computers kritisch zu überprüfen; -)

Auch Differentialgleichungen höherer Ordnung sind kein Problem:

DSolve[{x''[t]==x'[t]+6x[t],x[0]==1,x'[0]==8},x[t],t]

Parabolspiegel

Die Lösung der Differentialgleichung des Parabolspiegels:

DSolve[y'[x]==(((y'[x]*y[x]+x))/(Sqrt[((y[x])^2)+(x^2)])),y[x],x]

Systeme von Differentialgleichungen

Der DSolve-Befehl kann zwar theoretisch auch Systeme lösen, meistens ist aber eine analytische Lösung nicht möglich. Man muss dann das System numerisch mit NDSolve lösen.

Eine Lösungskurve x(t) = (x(t), y(t)) des Volterra-Lotka-Systems zu den Anfangsbedingungen (x(0), y(0)) = (0.7, 0.7) im Intervall t ∈ [0, 7] erhalten wir mit

loesung= {x[t],y[t]}/.NDSolve[{x'[t]==(1-y[t])x[t],y'[t]==(x[t]-1)y[t],x[0]==0.7,y[0]==0.7},{x,y},{t,0,7}][[1]];

ParametricPlot[loesung,{t,0,7}]

Folgendes Programm berechnet die Lösungskurven für mehrere Anfangswerte und stellt sie gemeinsam dar.

PhasePlot[f_,ic_,tmax_,opts___]:=Block[{i,n=Length[ic],ff,ivp,sol,phaseplot},
ff=f /. {x->x[t],y->y[t]};
Off[ParametricPlot::"ppcom"];
Do[
ivp={x'[t]==ff[[1]],y'[t]==ff[[2]],x[0]==ic[[i,1]],y[0]==ic[[i,2]]};
sol=NDSolve[ivp,{x[t],y[t]},{t,-tmax,tmax}];
phaseplot[i]=ParametricPlot[{x[t],y[t]} /. sol,{t,-tmax,tmax}]
,{i,1,n}];
On[ParametricPlot::"ppcom"];
Show[Table[phaseplot[i],{i,1,n}],opts]
];

Und so wird der Befehl verwendet (die Anfangsbedingung (x(0), y(0)) = (1, 1) entspricht übrigens dem Fixpunkt):

PhasePlot[{(1-y)x,(x-1)y},{{0.3,0.3},{0.5,0.5},{0.7,0.7},{0.9,0.9},{1,1}},3.7,PlotRange->{{0,3.5},{0,3.5}},Ticks->False]

Analog kann die Lösung der Lorenz-Gleichung mit dem Befehl

σ=10; ρ=28; β=8/3;loesung={x[t],y[t],z[t]} /.NDSolve[{x'[t]==-σ (x[t]-y[t]),y'[t]==-x[t] z[t]+ρ x[t]-y[t],z'[t]==x[t] y[t]-β z[t],x[0]==30,y[0]==10,z[0]==40},{x,y,z},{t,0,20}][[1]];

berechnet werden.

EIne grafische Darstellung erhält man mit

ParametricPlot3D[loesung,{t,0,20},PlotPoints->2000,Axes->False]

Beschreibende Statistik und Zusammenhangsanalysen

Stabdiagramm

Ein Stabdiagramm kann mit dem Befehl

BarChart[{1/20,3/20,7/20,1/4,1/10,1/10},BarSpacing->0.5,ChartLabels->{"397","398","399","400","401","402"},AxesLabel->{"(Subscript[a, i])","(Subscript[f, i])"}]

gezeichnet werden.

Verteilungsfunktion

Eine Verteilungsfunktion kann wie folgt dargestellt werden. Zuerst summieren wir die Werte der Reihe nach auf:

F=Accumulate[{1/20,3/20,7/20,1/4,1/10,1/10}]

Das sind die Funktionswerte an den einzelnen Punkten. Nun können wir sie wie folgt darstellen:

DiscretePlot[F[[n-396]],{n,397,402},ExtentSize->{None,Full},PlotMarkers->"Point"]

Histogramm

Ein Histogramm erhalten wir so:

Histogram[{3,7,12,18,19,20,25,25,27,28,29,31,32,34,37,38,40,41,45,47},{{0,10,20,30,40,50}},Ticks->{{5,15,25,35,45},{2,4,6} }]

Das erste Argument enthält die Urliste und das zweite (optional) die Klassengrenzen.

Arithmetische Mittel

Das arithmetische Mittel kann mit dem Befehl Mean berechnet werden :

Mean[{1,2,2,2,3,3,3,3,3,4,4}]

Median und Quantile

Den Median berechnen wir so :

data={1,2,2,2,3,3,3,3,3,4,4};
Median[data]

Ein p-Quantil (bzw. mehrere gleichzeitig, wenn man die verschiedenen p-Werte als Liste übergibt) erhalten wir mit :

Quantile[data,{0.25,0.5,0.75}]

(Achtung : Mathematica verwendet eine leicht abweichende Definition für das p-Quantil.)

Eine Box Plot erhält man mit :

BoxWhiskerChart[{3,7,12,18,19,20,25,25,27,28,29,31,32,34,37,38,40,41,45,47}]

Standardabweichung

Varianz und Standardabweichung können mit Variance bzw. StandardDeviation berechnet werden:

data = {1,1,2,2,3,3,3,3,4,4,7};
{Variance[data], StandardDeviation[data]}

Lorenz-Kurve und Gini-Koeffizient

Lorenz-Kurve für die häufigsten Automarken in Österreich 2012

data={187244,204044,210476,222903,225007,227863,275077,323840,344505,942268};
n=Length[data];
adata=Join[{{0,0}},Transpose[{Table[i,{i,n}],Accumulate[data]}]];
f=Interpolation[adata,InterpolationOrder->1];

Show[Plot[{((f[n])/n)x,f[x]},{x,0,n}],ListPlot[adata]]

Gini=((n f[n]/2-NIntegrate[f[x],{x,0,n}])/(n f[n]/2))//N

Lorenz - Kurve für die Einkommensverteilung in Österreich 2011

data={370,667,1049,1317,1540,1733,1933,2182,2537,3200,5444};
data=Join[{0},Table[((data[[i]]+data[[i+1]]-1)/2.),{i,Length[data]-1}]]

Mean[%[[2;;-1]]]

n=Length[data];
adata=Transpose[{Table[10(i-1),{i,n}],Accumulate[data]}];
f=Interpolation[adata,InterpolationOrder->1];

Show[Plot[{((f[100])/100)x,f[x]},{x,0,100}],ListPlot[adata]]

Gini=((NIntegrate[((f[100])/100)x-f[x],{x,0,100}])/(f[100]100/2))

Streudiagramm

Ein Streudiagramm erhält man leicht auf folgende Weise: Die Stichprobendaten werden als Liste eingegeben (z.B. mit dem Variablennamen data), danach wird mit dem Befehl ListPlot die Graphik erzeugt:

data={{163, 59},{165, 62},{166,65},{169,69},{170, 65}, {171, 69}, {171, 76}, {173, 73}, {174, 75},{175, 73} , {177, 80} , {177, 71},  {179, 82}, {180, 84}, {185, 81}};
ListPlot[data, PlotRange->{{160, 187},{55, 85}},AxesLabel->{Groesse x, Gewicht y},PlotStyle->PointSize[0.016]]

(Hier wurde – was nicht unbedingt notwendig ist – der dargestellte x,y-Bereich mit PlotRange gewählt, die Achsen wurden beschriftet und die Punkte vergrößert.)

Korrelationskoeffizient

Nach dem Laden des Zusatzpakets

Needs["MultivariateStatistics`"]

kann der empirische Korrelationskoeffizient mit dem Befehl Correlation berechnet werden:

xdata={30,45,45,60,75,80,90,100,110,120};
ydata={6,5,3,4,3,5,2,4,2,3};
Correlation[xdata,ydata]//N

xdata={163,165,166,169,170,171,171,173,174,175,177,177,179,180,185.};
ydata={59,62,65,69,65,69,76,73,75,73,80,71,82,84,81.};
Correlation[xdata,ydata]//N

xdata={163,165,166,169,170,171,171,173,174,175,177,177,179,180,185};ydata={2900,1100,3600,2300,4000,5600,2100,5100,3400,1800,2100,2600,4600,3600,2300};Correlation[xdata,ydata]//N

Regressionsgerade

Der Befehl Fit gibt uns eine Regressionsgerade aus :

data={{163, 59},{165, 62},{166,65},{169,69},{170, 65}, {171, 69}, {171, 76}, {173, 73}, {174, 75},{175, 73} , {177, 80} , {177, 71},  {179, 82}, {180, 84}, {185, 81}};
g=Fit[data, {1,x},x]

Die Syntax ist Fit[Daten, Funktionen, Variablen]. Dadurch werden Daten mithilfe der Gauß'schen Methode der kleinsten Quadrate durch eine Linearkombination der Funktionen genähert. Hier möchten wir eine Gerade, also eine Linearkombination der Funktionen 1 und x.

Nun können wir die Regressionsgerade zusammen mit den Daten darstellen:

g1=ListPlot[data,PlotRange->{{160, 187},{55, 85}},PlotStyle->PointSize[0.016]];
g2=Plot[g,{x,160,190}];Show[g1,g2, AxesLabel->{"Gewicht x", "Größe y"}]

Natürlich können wir auch direkt die Formeln aus dem Buch verwenden:

xdata=Transpose[data][[1]];
ydata=Transpose[data][[2]];
{mx,my}=Mean/@{xdata,ydata};
{sx,sy}=StandardDeviation/@{xdata,ydata};
r=Correlation[xdata,ydata];
k=r (sy/sx);
d=my-k mx;
{mx,my,sx,sy,r,k,d}//N

Clear[Gini, F,mx,my,sx,sy,r,k,d,f,g,g1,g2,data,xdata,ydata]

Anscombe-Quartett

Eine Darstellung der Anscomb-Quartett:

AnscombeA={{10,8.04},{8,6.95},{13,7.58},{9,8.81},{11,8.33},{14,9.96},{6,7.24},{4,4.26},{12,10.84},{7,4.82},{5,5.68}};AnscombeB={{10,9.14},{8,8.14},{13,8.74},{9,8.77},{11,9.26},{14,8.1},{6,6.13},{4,3.1},{12,9.13},{7,7.26},{5,4.74}};AnscombeC={{10,7.46},{8,6.77},{13,12.74},{9,7.11},{11,7.81},{14,8.84},{6,6.08},{4,5.39},{12,8.15},{7,6.42},{5,5.73}};AnscombeD={{8,6.58},{8,5.76},{8,7.71},{8,8.84},{8,8.47},{8,7.04},{8,5.25},{8,5.56},{8,7.91},{8,6.89},{19,12.5}};

Anscombe[data_]:=Block[{g},
g[x_]=Fit[data, {1,x},x];
Show[ListPlot[data,PlotRange->{{2,20},{0,14}},PlotStyle->{Orange,PointSize[Large]},AspectRatio->1,Frame->True,Axes->False],Plot[g[x],{x,2,20}]]
]

GraphicsGrid[{{Anscombe[AnscombeA],Anscombe[AnscombeB]},{Anscombe[AnscombeC],Anscombe[AnscombeD]}}]

Elementare Wahrscheinlichkeitsrechnung

Optimale Stoppstrategie

Die Wahrscheinlichkeiten für die optimale Stoppstrategie

P[n_,j_]:=(j/n)(∑_(k=j)^(n-1))(1/k);

können mit Mathematica leicht berechnet werden:

Table[N[P[12,j],2],{j,11}]

Clear[P];

Zufallsvariablen

FilledPlot

Ein Bereich einer Verteilung kann wie folgt veranschaulicht werden:

f[x_]:=(2/(Sqrt[π]))(x^2) Exp[-(x^2)];

gr1=Plot[f[x],{x,-4,1},Filling->Axis];
gr2=Plot[f[x],{x,-4,3}];
Show[gr1,gr2,AxesLabel->{"t","f(t)"},PlotRange->All]

Zufallszahlen

Zur Simulation von komplexen Zufallsvorgängen (z.B. bei der Planung von Großprojekten) benötigt man Zufallszahlen, die bestimmten Verteilungen genügen. In Mathematica können im Intervall [0, 1] gleichverteilte Zufallszahlen mit dem Befehl Random[] (ohne Argument) erhalten werden.

Mit dem Befehl

RandomReal[ExponentialDistribution[1.5]]

können exponentialverteilte (k = 1.5) Zufallszahlen erhalten werden. Es stehen noch weitere Verteilungen zur Verfügung, die wir später noch kennen lernen werden.

Eine Veranschaulichung des Hauptsatz der Statstik anhand eines Computerexperiments mit einer exponentialverteilten (k = 1) Zufallsvariablen erhalten wir wie folgt:

n=4;xlist=RandomReal[ExponentialDistribution[1],n];
F[x_]:=1-(E^(-x));
Fn[x_]:=Count[xlist,y_/;(y≤ x)]/Length[xlist];
Plot[{Fn[x],F[x]},{x,0,10}]

Der Befehl Count zählt die Anzahl der Elemente y in der Liste, die die Bedingung y <= x erfüllen. Division durch die gesamte Anzahl (= Länge der Liste) ergibt die empirische Verteilungsfunktion.

Clear[f,F,Fn,n]

Diskrete Wahrscheinlichkeitsverteilungen

Allgemein gilt: Der Erwartungswert wird mit Mean, die Varianz mit Variance berechnet. Die Wahrscheinlichkeit P(X=x) erhalten wir mit PDF („ probability density function ") und die Wahrscheinlichkeit P(X<=x) wird mit dem Befehl CDF („cumulative distribution function ") aufgerufen.

Hypergeometrische Verteilung

Mit HypergeometricDistribution[n, M, N] erhalten wir (symbolisch) die hypergeometrische Verteilung mit den Parametern n, M, N. Wir können nun z.B. den Erwartungswert und die Varianz berechnen :

hdist=HypergeometricDistribution[5,4,20];
{Mean[hdist], Variance[hdist]}

Hier wurde als Abkürzung für die Verteilung die Variable hdist eingeführt. Die Wahrscheinlichkeit P(X=2) bzw.die Verteilungsfunktion F(2) = P(X<=2) erhalten wir mit :

{PDF[hdist,2],CDF[hdist,2]}//N

Binomialverteilung

Mit BinomialDistribution[n, p] kann die Binomialverteilung aufgerufen werden:

bdist=BinomialDistribution[50,0.03];PDF[bdist,2]

Poisson-Verteilung

Der Befehl für die Poisson - Verteilung ist PoissonDistribution[λ] :

pdist=PoissonDistribution[1.6];
CDF[pdist,2]

Stetige Wahrscheinlichkeitsverteilungen

Wie schon bei den diskreten Verteilungen kann der Erwartungswert mit Mean und die Varianz mit Variance berechnet werden. Die Anweisung CDF („cumulative distribution function“) gibt die Verteilungsfunktion F und der Befehl PDF („probability density function“) liefert die Dichtefunktion f.

Normalverteilung

Mit NormalDistribution[μ, σ] erhält man (symbolisch) die Normalverteilung mit den Parametern μ, σ :

ndist = NormalDistribution[4,2];
CDF[ndist,6]//N

(Wir haben diese Normalverteilung mit ndist abgekürzt, um sie z.B. handlicher als Argument des CDF-Befehls schreiben zu können.) Wenn Sie nicht numerisch rechnen (also // N weglassen), so bekommen Sie einen symbolischen Output, der hier die Errorfunktion erf(x) enthält. Sie ist nichts anderes als eine Abkürzung für das (nur numerisch berechenbare) Integral erf(x) = (2/(Sqrt[π]))(∫_0^x)(e^(-(t^2)))dt.

p-Quantile

Das p - Quantil einer Verteilung erhalten wir mit dem Befehl Quantile[verteilung, p].Für eine Normalverteilung ergibt sich zum Beispiel :

Quantile[NormalDistribution[100,5],0.05]

Normalverteilung als Näherung

Näherung der Binomialverteilung:

n=(10^7);p=(10^(-5));μ=n p; σ=(Sqrt[n p (1-p)]);x=90;
1-{CDF[BinomialDistribution[n,p],x],CDF[NormalDistribution[μ,σ],x+0.5]}//N

Näherung der Poisson-Verteilung:

λ=20;x=14;{CDF[PoissonDistribution[λ],x],CDF[NormalDistribution[λ,(Sqrt[λ])],x+0.5]}//N

Prüfverteilungen

Das 0.9-Quantil der Chi-Quadrat-Verteilung mit m = 12 Freiheitsgraden ist

Quantile[ChiSquareDistribution[12], 0.9]

Das 0.9-Quantil der t-Verteilung mit m = 12 Freiheitsgraden erhalten wir mit

Quantile[StudentTDistribution[12],0.9]

Das 0.95-Quantil der F-Verteilung mit (Subscript[m, 1]) = 2 und (Subscript[m, 2]) = 12 Freiheitsgraden berechnet sich so:

Quantile[FRatioDistribution[2,12], 0.95]

Clear[ndist, n,p,μ, σ, x]

Schließende Statistik

Funktionen zur Berechnung von Konfidenzintervallen stehen nach dem Laden des Zusatzpakets

Needs["HypothesisTesting`"]

zur Verfügung. Die Beispieldaten aus dem Buch sind

data={100.,97.,101.,96.,98.,102.,96.,100.,101.,98.};
{Mean[data],Variance[data],StandardDeviation[data]}

Konfidenzintervall für den Erwartungswert

Das Konfidenzintervall einer Stichprobe aus einer normalverteilten Grundgesamtheit mit bekannter Varianz erhalten wir mit:

MeanCI[data,KnownVariance->4,ConfidenceLevel->.9]

Ohne Angabe des Konfidenzniveaus (ConfidenceLevel) wird ein Niveau von 1 − α = 0.95 verwendet :

MeanCI[data,KnownVariance->4]

Wird die Varianz nicht angegeben, so wird sie automatisch als unbekannt angenommen:

MeanCI[data]

Konfidenzintervall für den Vergleich zweier Erwartungswerte

Das Konfidenzintervall für den Vergleich zweier Erwartungswerte bekommen wir so:

data1={152,148,149,146,146,152,149,150};
data2={153,150,149,149,152,151,154,152};
MeanDifferenceCI[data1,data2,EqualVariances->True,ConfidenceLevel->0.9]

Konfidenzintervall für die Varianz

Das Konfidenzintervall für die Varianz bekommen wir mit

VarianceCI[data]

Konfidenzintervall für den Anteilswert

Ein Befehl zur Berechnung des Konfidenzintervalls für einen Anteilswert steht zwar nicht zur Verfügung, wir können aber leicht einen definieren. Da Mathematica die F-Verteilung kennt, können wir die exakte Formel implementieren :

ProportionCI[n_,x_,opts___]:=Block[{p,Fp},
p =1-((1-ConfidenceLevel)/2) /. {opts,ConfidenceLevel->.95};
Fp[m1_,m2_]:=Quantile[FRatioDistribution[m1,m2],p];
{(x/(x+(n-x+1)Fp[2(n-x+1),2x])),(((x+1)Fp[2(x+1),2(n-x)])/(n-x+(x+1)Fp[2(x+1),2(n-x)]))}];

ProportionCI[500,60]

Hypothesentests

Funktionen zur Durchführung von Hypothesentests stehen nach dem Laden des Zusatzpakets

Needs["HypothesisTesting`"]

zur Verfügung. Den Gauß-Test für H0 : μ >= 100 bekommen wir so :

MeanTest[data,100,KnownVariance->4,SignificanceLevel->0.05]

MeanTest[data,100,TwoSided->True,FullReport->True,SignificanceLevel->0.05]

Mathematica führt automatisch einen einseitigen Test durch (je nachdem, ob x > μ0 oder x < μ0 gilt) und gibt den zugehörigen p-Wert aus. Wird die Varianz nicht angegeben, so wird automatisch ein t-Test durchgeführt. Einen zweiseitigen Test können wir durch Angabe der Option TwoSided → True erhalten. Mit der Option FullReport → True werden zusätzlich Stichproben-Mittelwert, Prüfwert und die verwendete Verteilung ausgegeben.

(χ^2)-Anpassungstest auf Normalverteilung

Wir beginnen mit den gegeben Füllmengen:

xlist={102.5,115.7,93.8,102.6,110.,111.5,98.6,96.8,101.7,110.9,103.,104.4,108.4,104.2,97.7,105.5,92.1,100.3,97.9,105.5,106.1,90.3,108.9,96.4,90.4,91.5,94.9,99.9,84.9,102.,101.5,96.8,99.9,104.6,92.7,87.9,104.,108.6,94.7,107.3,98.6,96.6,105.4,101.4,104.,94.2,108.3,106.2,101.1,109.1,94.,95.6,100.1,89.5,101.2,94.1,92.,100.1,105.5,105.1,94.1,113.1,101.6,86.,92.1,91.5,98.6,90.6,101.4,93.6,88.3,88.5,88.6,95.9,108.2,101.2,101.1,96.9,100.2,104.7,96.6,109.2,108.5,108.4,111.6,99.2,90.8,111.7,99.7,100.4,96.,95.7,90.9,95.5,106.6,100.2,114.1,101.6,113.6,98.5};
n=Length[xlist]

Die gesuchten Schätzwerte sind

{m,s}={Mean[xlist],((n-1)/n)(Sqrt[Variance[xlist]])}

und die zugehörige Normalverteilung ist

ndist = NormalDistribution[m,s];

Nun berechnen wir für gegebene Klassengrenzen alle benötigten Werte:

k=16;
{min,max}={Floor[Min[xlist]],Ceiling[Max[xlist]]};
clist=Table[If[j≤0,-∞,If[j≥16,∞,min+(max-min)(j/k)]],{j,0,k}]

tab=Table[{{clist[[i]],clist[[i+1]]},Count[xlist,y_/;(clist[[i]]<y≤ clist[[i+1]])],n(CDF[ndist,clist[[i+1]]]-CDF[ndist,clist[[i]]])},{i,1,Length[clist]-1}];
tab//MatrixForm

Zusammenlegung einiger Klassen:

clist={-∞,90,92,94,96,98,100,102,104,106,108,110,∞};

tab=Table[{{clist[[i]],clist[[i+1]]},Count[xlist,y_/;(clist[[i]]<y≤ clist[[i+1]])],n(CDF[ndist,clist[[i+1]]]-CDF[ndist,clist[[i]]])},{i,1,Length[clist]-1}];
tab//MatrixForm

Daraus folgt der gesuchte Prüfwert und das gesuchte Quantil:

{Sum[(((tab[[i,2]])^2)/(tab[[i,3]])),{i,1,Length[clist]-1}]-n,Quantile[ChiSquareDistribution[Length[clist]-4],0.95]}