Statistik mit Mathematica

Copyright

Copyright 1994,  Claudia Funke
TU-Berlin, FB 13, Fachgebiet Ökonometrie und Statistik
Dieses Notebook darf ausschließlich als Unterrichtsmaterial und für private 
Zwecke verwendet und nicht ohne Zustimmung der Autorin verändert werden.

Ausgangsituation

Nach den jüngsten Veröffentlichungen über die Schädlichkeit des Rauchens
macht sich die Firmenleitung einer größeren Gesellschaft Gedanken, wie sie 
die bei ihr beschäftigten Raucher und Raucherinnen positiv beeinflußen kann.
Dazu läßt sie sich von ihrer statistischen Abteilung mittels einer 
repräsentativen Stichprobe das Rauchverhalten ihrer Beschäftigten ermitteln.

Diese wurden dementsprechend in fünf von der Firmenleitung vorgegebenen
Beschäftigungsgruppen aufgeteilt befragt. Die Werte wurden anschließend
in zuvor definierte Raucherklassen eingeordnet:

nicht = 	         0 Zigaretten / Tag
leicht =  	   1 - 3 Zigaretten / Tag
mittel = 	 4 - 19 Zigaretten / Tag 
stark =      19 - 35 Zigaretten / Tag

Es ergab sich danach folgende Gesamtübersicht:

RAUCHERTYPUS
	               	nicht	leicht 	mittel	stark	
---------------------------
Jüng.Arb.		36	48	66	26
Ält. Arb.		50	20	24	0
Jüng. Ang.	0	6	14	0
Ält. Ang. 		8	4	6	4
Management	20	12	14	4
------------------------------------------------------------------------------------

Die MitarbeiterInnen der statistischen Abteilung haben nun die Aufgabe, 
diese Daten statisch aufzubereiten.

Deskriptive Statistik

Einlesen und Darstellen eines Datensatzes

Die als ASCII - File vorliegenden numerischen Daten werden mit

ReadList[ "file", Number, RecordLists-> True  ]

in Mathematica eingelesen. Diese spezielle Form der Funktion  ReadList[ ]
liest eine Folge von Zahlen  (Number)  ein, die durch ein Leerzeichen 
voneinander getrennt sind. Die Option    RecordLists ->  True  bewirkt, 
daß ein als zweidimensionale Matrix in genau diesem Format eingelesen
wird.

xKlass = ReadList["rauch.dat", Number, RecordLists -> True];
MatrixForm[xKlass]

Im Gegensatz zur Matrixdarstellung läßt die Darstellung einer 
mehrdimensionalen Liste mit

TableForm[ Liste ]

mehr Gestaltungsmöglichkeiten zu, z.B. Bezeichnungen für die Zeilen
und Spalten einer Matrix:

TableForm[ Liste, TableHeadings -> {"Text1", "Text2", ... } ]

TableForm[xKlass, TableHeadings ->         { {"jung. Arb.","alt. Arb.","jung. Angst.",
          "alt. Angst.","Manag."},
         {"nicht","leicht","mittel", "stark"} },        TableDirections -> {Column, Row},        TableAlignments -> {Center,Right}   ]

Analyse des Datensatzes

Einen repräsentativen Überblick über die Situation der Nichtraucher und
Nichtraucherinnen erhält man, indem man die Daten für die Gesamt-
belegschaft betrachtet.

Umstrukturierung des Datensatzes

Ddie Daten der gesamten Belegschaft (für die einzelnen Rauchertypen)
ergeben sich aus der Randhäufigkeit der Spalten:

n  
h_j = Sum   h_ij
i=1

mit h_ij =  Häufigkeit einer Merkmalkombination in der  i - ten Zeile und
der j - ten Spalte.

Die Summen aller Spalten der Datenmatrix erhält man in Mathematica
am elegantesten mit einer Funktion, die die benötigte Operation auf
mehrere Listen anwendbar macht. Mit

Apply[ Plus, Liste ]

oder

Plus @@ Liste

werden die Elemente jeder Liste in  list  aufaddiert.

belegKlass = Apply[Plus, xKlass]

Der Umfang der Stichprobe wiederum ist die Summe der Randhäufigkeiten.

k
N  = Sum h_j
   j=1

spumfang = Sum[belegKlass[[i]], 
               {i, Length[belegKlass]}
           ]

Die neue Zeile   belegKlass   läßt sich mit der Funktion   AppendRows[ ]
aus dem Package     LinearAlgebra`MatrixManipulation`  in die 
Datenmatrix   X  einfügen.

<<LinearAlgebra`MatrixManipulation`

xKlassCols = Append[xKlass, belegKlass]; TableForm[xKlassCols, TableHeadings ->       { {"jung. Arb.","alt. Arb.","jung. Angst.",
        "alt. Angst.","Manag.","Gesamt"},        {"nicht","leicht","mittel","stark"} },
     TableDirections -> {Column, Row},      TableAlignments -> {Center,Right} ]

Graphische Darstellungsmöglichkeiten

Die Werte für die Gesamtbelegschaft lassen sich anschaulich
in einem Balkendiagramm (absolute Häufigkeit) und in einem
Tortendiagramm (relative Häufigkeit) darstellen. Die entsprechenden
Funktionen sind im Package    Graphics`Graphics`   enthalten.

<<Graphics`Graphics`

Mit

BarChart[ {y1, y2, ... } ]

wird ein Balkendiagramm der Beobachtungswerte y_i,  i = 1, 2, ... , n erzeugt.

Da   BarChart[ ]  ein Objekt der Art  Graphics  erzeugt, können alle dieser
Objektart zugeordneten Optionen bei der Gestaltung der Graphik verwendet
werden.  Die Option   Ticks -> {xticks, yticks}   ermöglicht z.B. die Bennenung
einzelner Balken, die Option   PlotLabel -> "String"  die Beschriftung der
Graphik.

barpic = 
BarChart[belegKlass, Ticks -> { {{1,"nicht"}, {2,"leicht"}, {3,"mittel"}, {4,"stark"}},    Automatic}, PlotLabel -> "Rauchertypen" ]

Mit der Funktion InputForm [ ] kann man hier sehr schön sehen, wie 
Mathematica solch eine Graphik aus einzelnen Graphikelementen 
aufbaut:

InputForm[barpic]

Mit

PieChart[ {y1, y2, ... } ]

wird ein Tortendiagramm für die Raucherklassen bezüglich der Gesamt-
belegschaft erzeugt.

tortpic =
PieChart[{ {belegKlass[[1]], "nicht"}, 
           {belegKlass[[2]], "leicht"}, 
           {belegKlass[[3]], "mittel"},  
           {belegKlass[[4]], "stark"} } 
        ]

Statistische Kenngrößen

Der erste graphische Eindruck soll nun durch statistische Kenngrößen 
quantifiziert werden.

Doch dazu werden statt der klassierten Daten die ursprünglichen 
Stichprobenwerte verwandt (Zigaretten / Tag).  Dieser Datensatz
wird wiederum mit   ReadList[ ]  eingelesen:

x = ReadList["rauchorg.dat", Number, 
              RecordLists -> True]

Dimensions[x]

Die Werte für die Gesamtbelegschaft lauten damit

beleg = Apply[Plus, x]

Zur Probe wird noch einmal der Stichprobenumfang berechnet.

spumfang = Sum[beleg[[i]], {i, Length[beleg]}]

Die Liste der Ursprungswerte enthält noch jeweils in der ersten Spalte 
die Zahl der Nichtraucher und Nichtraucherinnen, die für die folgenden 
Betrachtungen nicht mehr relevant sind. Mit

Drop[ Liste, n ]

werden die ersten n-Elemente aus der Liste  entfernt. Damit ergibt sich 
der Ursprungsdatensatz für die rauchende Belegschaft mit

belegRauch = Drop[Apply[Plus, x], 1]

Damit ändert sich auch der Länge des nun zu bearbeitenden 
Datensatzes auf

spumfang = Sum[belegRauch[[i]], 
               {i, Length[belegRauch]}
           ]

Für die folgenden Berechnungen wird auf verschiedene Packages aus 
dem Bereich Statistik zurückgegriffen, die alle unter dem Verzeichnis
....\statistics\    liegen. Durch das Laden des zugehörigen Package
Master.m  sind alle Funktionen der in diesem Verzeichnis enthaltenen
Packages verfügbar.

<< Statistics`Master`

Lageparameter

Die Lageparameter Mittelwert und Median werden mit den nun verfügbaren 
Funktionen

Mean[ Daten ]  ,    Median[ Daten ]

berechnet.

{mittel, med} = {Mean[belegRauch], Median[belegRauch]}; 
TableForm[N[{mittel, med}],
          TableHeadings -> { {"Mittelwert:","Median:"}, 
          {"Rauchende Belegschaft"}}]

Streuungsparameter

Für die wichtigsten Streuungsparameter wird in Mathematica ebenfalls eine
Funktion angeboten:

Dispersion[ Daten ]

Damit werden die Einzelwerte

	Variance[ ]
	StandardDeviation[ ]
	SampleRange[ ]
	MeanDeviation[ ]
	MedianDeviation[ ]
	QuartileDeviation[ ]

in einem Aufruf berechnet.

streupar = DispersionReport[beleg]; N[streupar]

Den nicht implementierten Variationskoeffizienten

v = s / x_bar

berechnet man mit

varkoef = (Variance /. streupar) / (Mean[beleg]);
N[varkoef]

Lineare Regression

Ausgangssituation

Der Firmenleitung ist aus verschiedenen Publikationen bekannt, daß
Lärm als ein Streß verursachender Faktor einen Einfluß auf das Rauch-
verhalten haben kann. Sie beauftragt daher die statistische Abteilung,
diesen Einflußfaktor quantitativ zu bewerten. Die statistische Abteilung
erhebt daraufhin bei 50 repräsentativ ausgesuchten Beschäftigten
(aus dem Produktionsbereich) folgende Daten:

- täglicher Zigarettenkonsum am Arbeitsplatz [Stück/Tag]
- durchschnittlicher täglicher Lärmpegel am Arbeitsplatz [dB]

Regressionsanalyse

Um einen ersten Überblick über die Art des funktionalen Zusammenhangs 
zwischen den Variablen Lärm und Zigarettenkonsum zu erhalten, fertigt
man einen Scatterplot an.

regdat = ReadList["reg.dat", Number, 
         RecordLists -> True]

reglistplot = 
ListPlot[regdat, Prolog -> AbsolutePointSize[2],
                 PlotRange -> {{55, 90}, {0,20} },
                 AxesLabel -> {"Laerm[dB]","Zig./Tag"},
                 PlotLabel -> "Scatterplot"
        ]

Ein nach dem Scatterplot zu vermutender linearer Zusammenhang 
zwischen dem Lärm (erklärende Variable, Regressor) und dem Zigaretten-
konsum (zu erklärende Variabe) wird in folgendem Ansatz modelliert:

y_i = beta0 + beta1 x1 + eps_i   mit e_i ~ N(0,1)

Die Funktionen zur Regressionsanalyse sind in Mathematica in dem
Package   Statistics`LinearRegression`    enthalten. Diese Funktionen
stehen durch das Einlesen des Packages   Master.m   aus dem
Verzeichnis   ...\statistics\       schon zur Verfügung. Mit

Regress[ Daten, {1, Funktionsform}, x ]

wird ein Modell der Form  Funktionsform  mit der Variablen  x  den Daten 
angepaßt. Für die Datenmatrix  Daten   wird angenommen, daß die 
letzte Spalte (das letzte Element jeder Liste)  die Werte der abhängigen
Variablen enthält.

Das Symbol x bekam weiter oben schon einen Wert zugewiesen. Diese 
Zuweisung und x an sich ohne Zuweisung wird mit ClearAll[ ]  gelöscht.

ClearAll[x]

Regress[regdat, {1, x}, x]

Mit der Option    OutputControl   ->  NoPrint    und eventuell darauffolgende 
Variablennamen kann die Ausgabe der Funktion    Regress[ ]  gezielt
gesteuert werden.

output = 
Regress[regdat, {1, ex}, ex, OutputControl -> NoPrint, OutputList -> {FitResiduals,BestFitCoefficients}];

Die Funktionsform erhält man mit

Fit[ Daten, Funktionsform, x  ]

func[x_] = Fit[regdat, {1, x}, x]

Die graphische Darstellung der Regressionsgeraden  ergibt eine
erst einmal optisch nachvollziehbare Modellwahl.

regplot = 
Plot[func[x],                                                       
     {x, Min[Transpose[regdat] [[1]] ], 
         Max[Transpose[regdat] [[1]] ]},
     DisplayFunction -> Identity                                     
    ]

regall = 
Show[reglistplot, regplot, 
     PlotLabel -> "Beob., Regression",
     AxesLabel -> {"Laerm[dB]","Zig./Tag"}
    ]

Die Gültigkeit des geschätzten Regressionsmodells hängt im Wesentlichen
von der Übereinstimmung der theoretischen Vorgaben (Modellanahmen) mit 
den vorliegenden Daten ab. 

Als Beispiel sei dazu die Analyse der Residuen betrachtet. Ein Plot der 
Residuen übepr