Ao. Prof. M. Hudec: 405 008 Statistische Software & Algorithmen

1998-03-17

1998-03-19

1998-03-24

1998-03-26

1998-03-31

1998-04-02

1998-04-21

1998-04-28

1998-05-07

1998-05-12

1998-05-19

1998-06-04

1998-06-09

1998-06-16


1998-03-17 Daten in S-PLUS

Informelle Datentypen:

Vektor: geordnete Menge von Elementen, Position eindeutig, alle Elemente haben den gleichen Typ

Matrix: Vektor, aber zweidimensional, alle Elemente haben den gleichen Typ

Atomare Datentypen: (Vektor, Matrix): können selbst keine atomaren Typen enthalten.

Liste: Elemente können unterschiedliche Typen haben, Liste kann andere Listen enthalten (rekursiv)

Funktion: Zusammenfassung von S-PLUS-Ausdrücken

Formale Klassen:

factor: Vektor, der kategorielle Werte repräsentiert (z.B. männlich/weiblich) (Klasse bekommt eine "Semantik", Bedeutung)

Objekte der Klasse lm: Datenobjekte, die das Ergebnis eines Fits eines linearen Modells repräsentieren.

print(x): verhält sich polymorph, erkennt die Eigenschaften der Objekte.

DataFrame: 2dimensionales Array, das eine klassische Datenmatrix repräsentiert (Spalten = Variablen; Zeilen = Beobachtungen), Spalten können verschiedenen mode aufweisen.

Objekte werden neben dem Typ durch ihren Modus charakterisiert.

Der mode gibt an, welche Werte von eienm Datenelement prinzipiell gespeichert werden können.

modezulässige WerteBeispiele storage mode
nullRealisierung der leeren MengeNULL null
logicalBoolesche WahrheitswerteT, F logical
numericganze Zahlen17, 4integer

Vektoren in S-PLUS

Vektoren (atomic objects):

Zusammenfassung elementarer Objekte des gleichen Modus zu einer Folge endlicher Länge, Elemente werden durch einen Index charakterisiert wird.

Konstruktion von Vektoren

elementarer_Datentyp(length=n) bzw. elementarer_Datentyp(n)

z.B. x<-numeric(10) ... x<-logical(3)

vector(mode="elementarer_Datentyp", length=n)

z.B.: x<-vector(mode="numeric", length=10)

Durch Zuweisungsausdruck ohne vorhergehende Deklaration:

z.B. x<-1:5

Prinzip: "self describing objects"

In Datenbjekten werden neben den Orthoinformationen auch Metainformationen (Attribute des Objekts) gespeichert.

mode(x)liefert Datentyp
length(x)liefert Länge des Vektors
is.vector(x)liefert booleschen Wert T, falls x ein Objekt vom Typ Vektor ist
is.vector(x, "numeric")

is.numeric(x)

liefert T, wenn x ein Objekt vom Typ Vektor und Mode numeric ist.

Alle Eigenschaften eines Objekts sind dynamisch und können im Verlauf einer Sitzung verändert werden:

> <-c("a", "b", "c") #c: concat

> y

> [1] "a" "b" "c"

> mode(x) <-mode(y)

> x

> length(x)<-length(y) #schneidet ab, wenn nötig

> as.numeric(x) #schreibt konvertierte Daten

> x<-numeric(x) #konvertiert x in numerische Daten

>y<-as.numeric(y)

Warning messages:

3 missing values generated coercing from character to numeric in: as.numeric(y)

>y

(1) NA NA NA

usw..

Subsetting

>x<-rnorm(10) #random normalverteilte Variablen

>x

(1) -1.087 .... ...

(8) -1.080 .... ...

>x[4]

(1) 1.738

>x[2:5]

(1) 1.053 ....

>x[c(2:5, 7:9)] #nicht zusammenhängender Bereich

>x[length(x)] #letzter Wert

>x[-(2:8)] #Alles, was nicht in diesem Bereich liegt

>x[c(T,F,T,F,T,F,T,F,T,F)] #Alle nehmen, die True sind

>x[c(T,F)] #Durch Recycling Rule dasselbe wie voriger Befehl

>x[x>0] #Alle, die größer als 0 sind

>min(x[x>0]) #kleinstes positives Element

>x[x>0][1:4] #die ersten vier positiven Elemente

>x<-1:3

>names(x)<-c("a", "b","c") #Variablen benennen

>x

a b c

1 2 3

>x["b"]<-5 #Variable wird über den Namen angesprochen

>x

a b c

1 5 3

Operationen mit Vektoren

Verketten

>c(x,y) #zwei Vektoren zusammenfügen

>append(x,y, StartPos) #Ab StartPos werden die Daten aus dem zweiten
Vektor eingefügt

Ersetzen

>x[3]<-9 #drittes Element ersetzen

>replace(x, 3, 9) #dasselbe wie vorher

>replace(x, c(3,5),8) #drittes und fünftes Element ersetzen

Folgen

>3.1:8

3.1 4.1 5.1 6.1 7.1

>seq(from=1, to=21, by=5)

1 6 11 16 21

>seq(from=1, by=5, length=6) #sechs Elemente

1 6 11 16 21 26

>seq(to=21, by=5, length=6) #sogar das funktioniert: from = -4

-4 1 6 11 16 21

>seq(from=1, to=21, length=5) #by=5

>seq(-5) # -5 ist to, default-from ist 1, default-by 1 oder -1

Wiederholungen

>rep(1,5) #replicate

1 1 1 1 1

>rep(1:3, 3:1) #wie oft jede Variable wiederholt wird

1 1 1 2 2 3

>rep(1:3, 3:1, 9) #length übersteuert den times-Parameter

Arithmetik

>x<-1:5

>x+2 #durch Recycling Rule: zu jedem Element 2 addieren

3 4 5 6 7

>1:3^2 #Achtung, das liefert 1 2 ... 9! Lieber Klammern verwenden.

Getelem<-function(x, from = 1, each = 1)

{
# Selektiert Elemente des Vektors x
# Beginnend von der Position from wird jedes each-te 
# Element ausgewählt
        sel(<-c(T, rep(F, each -1)) # eine Sequenz mit einem T und 
                        # der gewünschten Anzahl von F
        x[from:length(x)][sel] # Recycling ausnutzen












}

Permutationen

>x<-c("a","b","c","d")

>x

"a" "b" "c" "d"

>x[c(3,1,2,4)]

"c" "a" "b" "d"

>shift.left<-function(x, k)

{

# x ... any vector with length(x) >0
# k ... number of positions to be shifted
        k<-k%%length(x) # %%: modulo
        if(k!=0) # k ungleich 0
              x[c((k+1):length(x),1:k)]
        else x












}

Der Rückgabewert einer Funktion ist immer das letzte Statement.

Aufgaben: 2 (letter-Funktion); 3 (Funktion rle); 4; 7; 8; 9


1998-03-19

Ermitteln von Indexwerten

(An welcher Stelle kommt der ... Wert vor)

runif = random uniform
>(1:10)[x>0.5]  # indexing: ich bekomme die Positionen zurück
[1] 2 3 4 7 8 9 10
pos.max <-function(x)
{
        (1:length(x))[x==max(x)] # Vektor mit der gleichen Länge von 1 bis Anzahl;
                # alle Positionen, wo x = max(x) ist
}
pos.any<-function (x, any)
{
        (1:length(x))[any]      # Bedingung wird übergeben, z.B. x == min(x) oder x > 0.7
}













Reduktion von Vektoren

if (logische Bedingung) 
        führe X aus 
else 
        führe Y aus
>x
[1]5 9 2 7 2 7 8 9 7 ...
peaks(x, span=3, strict=T)      #lokale Maxima
[1] F T F T F F F
>any(x>5)
[1] T
>all(x>5)
[1] F
range(x) # kleinster und größter Wert
diff(range(x)) # Spannweite, Unterschied zwischen kleinstem und größtem Wert
cummax(x) # Bisheriges Maximum, monoton steigendes Ergebnis













Sortieren

>x

>sort(x) # wenn NAs existieren, werden sie weggeschnitten,

#Länge des Vektors ändert sich!

Umgekehrt sortieren

>y<-sort(x)[length(x):1]

>y<-rev(sort(x))

Sortliste

>sort.list(x) # gibt aus, wo die Daten wären, wenn sie sortiert würden

>order(x) # dasselbe, kann aber auch mehrfach sortieren (nach mehreren Spalten)

>rank(x) # inverses sort.list

Manipulationen von Listen

>x<-1:4
>y<-c("a", "b", "c")
>z<-c(T,F)
>l<-list(x,y,z) #generiert eine Liste aus 3 Komponenten
>l
[[1]]:
[1] 1 2 3 4
[[2]]:
[1] "a" "b" "c"
...
>l[[1]] #Zugriff auf einzelne Komponenten
[1] 1 2 3 4     #Ergebnis ist ein Vektor
>l[[1]][2]      #Zugriff auf Elemente der Komponenten
[1] 2
>l[2]   #Zugriff auf eine Subliste!
[[1]]:
[1] "a" "b" "c" # Eine Subliste wird zurückgeliefert, nicht der Vektor
>is.list(l[[2]]) # liefert False
>is.list(l[2]) #liefert True
is.recursive(list(1:10)) #liefert True, weil die Liste rekursiv sein könnte
is.recursive(1:10) #liefert False, weil der Vektor atomar ist













Listen als Rückgabewert von Funktionen

my.moments<-function(x=NULL, p=4)
{

if (is.null(x))
{ cat("Please enter your data \n")

x<-scan()

(....)

list(a,b)

}

else

list(NA, NA)

}


1998-03-24

Komponentenweise Operationen über Listen

>a1<-rnorm(round(runif(1)*100))
>a2<-rnorm(round(runif(1)*100))
(...)
>a<-list(a1, a2, a3, a4)
>lapply(a, length) #Befehl length auf die Liste a anwenden
[[1]]:
[1] 45
[[2]]:
[1] 95











...

>sapply(a, mean) # s für simple, Ergebnis nicht als Liste, sondern (diesmal) als Vektor

[1] -0.0592 ... 0.4092 ...

>lapply(a, range)

[[1]]:

[1] -2.948012 1.477

(...)

>sapply(a, range)

# jetzt kommt als einfachste Form eine Matrix heraus

[1][2]
[1]-2.94-2.48
[2]1.4772.249

Implementation von Datenmatrizen als Listen:

> s <- [irgendwas] usw.

> daten.matrix <- list(Jahr=j, Limit=l, Strasse=s, Unfallanz=u)

>daten.matrix

$Jahr: ....

$Limit: ....

>attach(daten.matrix, pos=1) #Bezeichner als Standardobjekt definieren

#(Nummer 1 in der Search-Liste)

>mean(Unfallanz)

[1] 62.5

>sapply(split(Unfallanz, Strasse), mean) #Split teilt die Daten nach den

#genannten Kategorien auf

Autobahn Landstr Nebenstr

30.25 63.5 93.75

Tricks mit Listen

Komplexe Folgen

>unlist(lapply(as.list(1.5), seq))

1 1 2 1 2 3 1 2 3 4 1 2 3 4 5

# seq(3): [1] 1 2 3

# as.list(1:5)

# lapply(as.list(1:5), seq): [[1]]: [1] 1 [[2]] [1] 1 2 ...

# unlist: Struktur wird aufgehoben

Runlength encoding

>x<-rep(1:3, c(8, 6, 16))

>x

[1] 1 1 1 1 ... 2 2 2 2 ... 3333 ...

>rx <-rle(x)

>rx

$lengths:

[1] 8 6 16

$values:

[1] 1 2 3

>rep(rx$values, rx$lengths)

[1] 1 1 1 ... 2 2 2 ... 3 3 3 ...

Matrizen

>A<-matrix(1:12, nrow=3, ncol=4)

>A

[,1][,2][,3][,4]
[,1]14710
[,2]25811
[,3]36912

#Die Daten werden spaltenweise eingefügt!

>mode(A)

"numeric"

>is.vector(A)

F

>is.matrix(A)

T

>dim(A) #Dimensionsattribute

3 4

>matrix(1:12, 3, 4) #dasselbe

>matrix(1:12, 3, 4, byrow = TRUE) #zeilenweise einfügen

1 2 3 4

5 6 7 8

9 10 11 12

>matrix(1:3 nrow=3, ncol=4)

1 1 1 1

2 2 2 2

3 3 3 3

>diag(3) #Diagonalmatrix mit der Dimension 3

1 0 0

0 1 0

0 0 1

# diag() mit numerischer Konstante generiert eine n*m Einheitsmatrix

# diag() mit numerischem Vektor generiert eine Diagonalmatrix aus den Werten

>diag(1:3)

1 0 0

0 2 0

0 0 3

>A<-matrix(1:9, 3, 3)

1 4 7

2 5 8

3 6 9

>diag(A) #Anwendung auf eine Matrix

[1] 1 5 9

>diag(diag(A))

1 0 0

0 5 0

0 0 9

>diag(A)<-0

>A

0 4 7

2 0 8

3 6 0

>A<-diag(2)

>B<-matrix(1:4, 2, byrow=TRUE)

>A

1 0

0 1

>B

1 2

3 4

>cbind(A, B) #column-bind

1 0 1 2

0 1 3 4

>rbind(A, B) #row-bind

1 0

0 1

1 2

3 4

>C<-B

>D<-A

>rbind(cbind(A,C), cbind(B, D))

1 0 1 2

0 1 3 4

1 2 1 0

3 4 0 1

Zugriff auf Elemente der Matrix

>A

1 2 3

2 3 1

3 1 2

>A[2, 2]

[1] 3

>A[1:2,]

1 2 3

2 3 1


1998-03-26

> A<-matrix(1:8,2,4)

> A

[,1] [,2] [,3] [,4]

[1,] 1 3 5 7

[2,] 2 4 6 8

> 2*A

[,1] [,2] [,3] [,4]

[1,] 2 6 10 14

[2,] 4 8 12 16

> 2*A #jedes Element wird mit 2 multipliziert

[,1] [,2] [,3] [,4]

[1,] 2 6 10 14

[2,] 4 8 12 16

> (1:2)*A# recycling rule schlaegt zu!

[,1] [,2] [,3] [,4]

[1,] 1 3 5 7

[2,] 4 8 12 16

> B<-matrix(1:4, 2,2)

> B

[,1] [,2]

[1,] 1 3

[2,] 2 4

> B*B

[,1] [,2]

[1,] 1 9

[2,] 4 16

> A*B # keine Recycling Rule für Matrizen!

Error in A * B: Dimension attributes do not match

Dumped

> B%*%A#echte Matrixmultiplikation

[,1] [,2] [,3] [,4]

[1,] 7 15 23 31

[2,] 10 22 34 46

> A%*%B

Error in "%*%.default"(A, B): Number of columns of x should be the same as number of rows of y

Dumped

> B%*%(1:2)

[,1]

[1,] 7

[2,] 10

> (1:2)%*%B

[,1] [,2]

[1,] 5 11

> t(A) # transponiert die Matrix

[,1] [,2]

[1,] 1 2

[2,] 3 4

[3,] 5 6

[4,] 7 8

> crossprod(A,B) #Matrixprodukt der Form A'B

[,1] [,2]

[1,] 5 11

[2,] 11 25

[3,] 17 39

[4,] 23 53

Was bedeutet sum(diag(crossprod(A))?

> tapply(x, col(x), sum) #col(x) gruppiert die Elemente des Objekts

1 2 3 4 5

10 26 42 58 74

> apply(x, 2, sum)

[1] 10 26 42 58 74

> t(A)

[,1] [,2] [,3] [,4] [,5]

[1,] 1 2 3 4 5

[2,] 6 7 8 9 10

[3,] 11 12 13 14 15

> A

[,1] [,2] [,3]

[1,] 1 6 11

[2,] 2 7 12

[3,] 3 8 13

[4,] 4 9 14

[5,] 5 10 15

> t(t(A)-apply(A,2,mean)) #elegant, verwendet Recycling Rule nach der Transformation

[,1] [,2] [,3]

[1,] -2 -2 -2

[2,] -1 -1 -1

[3,] 0 0 0

[4,] 1 1 1

[5,] 2 2 2

> sweep(A, 2, apply(A,2,mean))

[,1] [,2] [,3]

[1,] -2 -2 -2

[2,] -1 -1 -1

[3,] 0 0 0

[4,] 1 1 1

[5,] 2 2 2

sweep subtrahiert den Vektor im 3. Parameter aus der Matrix im ersten Parameter

Arrays

Arrays sind so etwas wie Matrizen, aber sie können höhere Dimensionsattribute haben.


1998-03-31

Bei Bereichsüberschreitungen beim Indizieren von Arrays warnt S-Plus nicht, sondern löst sie auf.

Graphiken mit S-Plus

Grafisches Ausgabefenster notwendig, Befehl: win.graph() oder win.printer()

Der plot()-Befehl verwendet immer das aktuelle Grafik-Gerät und löscht, was vorher darauf war. Es ist aber möglich, neue Grafik-Geräte zu öffnen.

Inkrementelle Grafikbefehle

Sie zeichnen in ein bestehendes Grafikgerät hinein, z.B. Achsen (axis()), Titel (title())

> plot(x,y,type="p")

> title(main="Hallo Welt")

> axis(1, pos=0, labels=F)

Zusätzliche Parameterübergabe

function(x, y, ...)

#

plot(x, y, ...)

Die zusätzlichen Parameter im Funktionsaufruf werden an die aufgerufene Funktion weitergereicht.

(Erben der Funktionalität von plot())

Die Funktion parse()

Mit parse() kann eine Zeichenkette daraufhin überprüft werden, ob sie den Regeln der Sprache entspricht. Das Ergebnis von parse() ist dann eine evaluierbare Funktion (mit eval()).

Zusätzliche Grafik auf einem Grafikdevice

Nach plot() können z.B. mit lines() zusätzliche Grafikelemente in einen Plot aufgenommen werden.


1998-04-02 Nachbesprechung der Übungsaufgaben

Aufpassen bei runif(): wenn der Output einfach mit einer Zahl multipliziert wird, kommt keine Gleichverteilung, weil der erste und der letzte Wert nur die halbe Wahrscheinlichkeit haben.

Denksport: 7. Nur einmal sortieren

2. Übung:

7. Funktion outer()

9. rekursiv lösen


1998-04-21

Die Funktion matplot() wird verwendet, um zweidimensionale Tabellen zu plotten. Sie paßt auch selbst die Grenzen an, was ein zusätzliches lines() nicht tut.

Klassen

summary(), print(), plot(): generische Funktionen; sie ermitteln die aufzurufende Funktion (z.B. summary.list(), summary.vector()) aus der Klasse der Parameter.

class(x.list) <-"konz" #der Variablen eine "Klasse" zuweisen

summary.konz() # klassenspezifische Implementation. Der Name muß zur Klasse passen, sonst wird die Standard-Summary-Funktion verwendet

Teilung der Zeichenfläche

par(mfrow=c(2,2)): Zeichenfläche aufteilen, jeder neue Grafikbefehl gilt für eine neue Teilfläche

In S-PLUS gibt es zwei Grafiklibraries, die "classic" und die "Trellis-Library". Trellis kann in der Regel deutlich mehr, die Syntax ist aber etwas anders.

Mehrdimensionale Darstellung

persp() #Perspektivansicht

contour() # "Höhenlinien"


1998-04-28 Höherdimensionale Grafiken

Möglichkeiten, höherdimensionale Daten darzustellen:

1. Daten als Funktionsparameter (z.B. faces())

2. Mehrere Scatterplots (z.B. pairs())

3. Dreidimensionale Darstellung (z.B. cloud() aus der Trellis-Library)

Rekursive Funktionen

Mit recall() kann die aktive Funktion rekursiv aufgerufen werden.

Dialogboxen

dialog.main<-function()

{

graph <- dialogbox(size=c(120,200), title="Fenstertitel", controls = list(

listbox.control(label="&Verteilungen",

listvalues=c("Normalverteilung", "Chiquadrat"),

multiplesel=F, selected=1, size=c(100,40))))

dialog.back<-dialogbox.display(graph)

usw.

}

Datenmatrix in S-PLUS

Interne Organisation sollte "Liste von Vektoren" sein. Damit können die einzelnen Variablen unterschiedliche Typen (z.B. m/w, numerisch) haben.

Eigenschaften einer Matrix in bezug auf Indexierung: "data frames"

Generierung:

attach(x)


1998-05-07

browser(): zum Debuggen in S-PLUS

Beispiel. 3: ML-Schätzer: x quer


1998-05-12

scan(dateiname, what="") #einen Wert aus einer Textdatei einlesen, herauskommt ein Vektor oder eine Matrix (je nach Anzahl der what-Parameter).

identify(): zeigt die Grafik zum Anklicken an, zeigt die Namen der Werte und liefert, welche Punkte ausgewählt wurden. Ausstieg mit rechter Maustaste.

abline(): Legt eine Gerade über eine Grafik; entweder aus Werten, oder aus einer Funktion wie lsfit (Kleinste Quadrate), l1fit (least-absolute-deviation) usw.


1998-05-19

1. Modellannahme: Linearität

y=Xb+e

(x*1): Daten

B.L.U.E. = best linear unbiased estimator

2. Modellannahme: Homoskealastizität: Varianz konstant

3. Modellannahme: Unkorreliert

Autokorrelation: z.B. saisonale Daten

Modelldiagnose

Punkte mit großen Residuen: "outlier" (Nichtlinearität)

Leverage points: influential data points (ob das Ergebnis von Einzeldatenpunkten abhängig ist)

Output von lsfit()

coef

residuals: Residuen der einzelnen Datenpunkte

wt: Gewichtsvektor

intercept: boolesch, ob die Regressionsgerade durch den Nullpunkt geht

qr: numerische Zerlegung

ls.print: Methode, um den Output von lsfit() "schön" auszugeben.

Hat-matrix: X(X'X)-1X'

ratio estimator = Verhältnisschätzer

mean(yi/xi) oder mean(xi)/mean(yi)


1998-06-04

y ~ x*A ^= y ~ x+A+x*A

"Reparametrisieren": contrasts()

contrast treatment: 0 und 1. Parameter: Abweichungen von der Nullstufe

contrast sum: -1 und +1.

Faktoren mit k Stufen generieren k - 1 Dummies (=Differenzen)


1998-06-09

Mallows CP-Statistik: Least Square-Analogon (Spezialfall von Akaike IC fürs lineare Modell)

Akaike IC: maximum likelihood-Prinzip

Suche nach dem besten Modell

1. All subsets (leaps())

2. stepwise

anova()

step() Automatische Berechnung des besten Modells

Richtungen des Modellbaus: forward (Anzahl der Variablen wird größer); backward (#Variablen wird kleiner); combined (alternierend)

stepwise(): überprüft selbst, ob eine "All subsets"-Analyse noch funktioniert ("exhaustive search") oder sucht das beste Modell.


1998-06-16

© Balázs Bárány
zuletzt geändert (JMT): 1999-10-01