Hier findet sich der gesamte Inhalt des Handbuchs zum Durchsuchen. Diese Seite wird nur unregelmäßig geupdated. Aktuelle Skripte finden sich auf der
Rclick Hauptseite!
Tipp: mit Pfeiltaste links/rechts kann horizontal gescrollt werden.
Syntax highlighting mit [*sourcecode language=“r“ gutter=“false“], dann den Codeblock, dann [*/sourcecode] (natürlich ohne Sternchen *), wie von WordPress vorgeschlagen.
# Komplettes Handbuch zum Durchsuchen. Im TinnR: STRG+F, weiter mit F3
# Stand: 2016-07-12 07:32:12
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 1: Grundlagen ######################################################
# 1.1 Hilfe aufrufen
# 1.2 Rechenoperatoren
# 1.3 Zuordnungen
# 1.4 Objektarten
# 1.5 Grund-Datentypen (Übersicht)
# 1.6 Objekttypen (Übersicht)
# Kommentare mit einem Kreuz (#) beginnen, werden von R nicht ausgewertet.
# Leerzeichen und -zeilen werden nicht ausgewertet; sie dienen der Übersicht.
# Mehrere Befehle in einer Zeile können mit ; voneinander getrennt werden.
# schon früh sollte man sich angewöhnen, bequem zu programmieren:
browseURL("http://RclickHandbuch.wordpress.com/install-r/tinn-r/")
# Besonders: TinnR Options > Shortcuts sowie Return focus to editor
# STRG + F : Find. Suche nach Wörtern. Weiter mit F3
# STRG + 9 schreibt Klammern und setzt den Cursor in deren Mitte
# Wer gerne mit VideoTutorials lernt, könnte hier schauen:
browseURL("https://www.youtube.com/user/Tutorlol")
# Ansonsten sei hier nochmals auf die Literaturliste im Anhang verwiesen...
# 1.1 Hilfe aufrufen -----------------------------------------------------------
# ? Befehl() oder help.search("Befehl")
?mean() # Befehlsteil in Klammern wird ignoriert
help.search("mean") # kürzer: ??mean
help.start() # HTML Hilfe (Automatischer Start des Browsers)
# am schnellsten: cursor auf Befehl im TinnR setzen, F1 drücken
# bei Befehlen mit Punkt: den ganzen Befehl zuerst markieren.
help.search("Kolmogorov Smirnov")
?ks.test
# 1.2 Rechenoperatoren ---------------------------------------------------------
7+11
3*9/5
3^2
16:21 # ganze Zahlen (Integer) von 16 bis 21 (als Vektor)
sqrt(81) # Wurzel (gleiches Ergebnis bei 81^0.5 )
abs(-12) # absolutwert (Betrag)
sin(2.4*pi) # ACHTUNG: keine Kommas (,) als Dezimaltrennzeichen verwenden!
exp(3) # e^3
log(100) # log(x) ist im R ln(x)
log10(100) # log10(x) ist der Logarithmus mit Basis 10
factorial(5) # 5! (5 Fakultät)
0.0002 # xx e yy bedeutet: xx * 10 ^ yy (Exponentialschreibweise,
3e5 # auch Scientific Notation genannt)
1.9343E+03 + 2
print(pi, digits=16) # print für eine benutzerdefinierte Ausgabe
options("digits"=16) # Ab jetzt jede Ausgabe mit 16 Nachkommastellen
sin(0.8)
options("digits"=7) # wieder zurück auf die DEFAULT = Standard Einstellung
-1/0 ; 0/0 ; 1/0 # Inf = Infinity = unendlich; NaN = Not a number
# 1.3 Zuordnungen (Englisch: Assignment) --------------------------------------
# Werte, Zeichenketten, Aussagen... als Variablen speichern (Objekten zuordnen)
# R unterscheidet Großschreibung!
probe1 <- 17.5 # Punkt als Dezimaltrennzeichen verwenden
# Anzeigen der zugeordneten Daten
probe1 # Das ist jetzt ein Objekt im Workspace, nur in dieser R-Sitzung.
probe1 + 12 # Jetzt kann damit auch weitergerechnet werden.
p2 <- 1 > 2 # Ist 1 größer als 2? Ergebnis wird als p2 gespeichert.
p2 # Symbole möglichst kurz aber wiedererkennbar wählen
p3 <- sin (5) # Werte in Radiant (pi = 180°) x = alpha*pi/180°
sin(45*pi/180) # sinus von 45°
p4 <- "ein Text" # Anführungsstriche nicht vergessen
p4
p5 <- TRUE # oder kürzer: p5 <- T (bzw F für FALSE)
sqrt(9030) -> p6 # Zuordnungen in beiden Richtungen möglich, aber "->"
# sollte vermieden werden, da es sehr ungewöhnlich, daher schlecht lesbar ist.
p7 = 7:15 # Zuordnen mit "=" ist möglich, aber ungünstig:
p7 # http://r.789695.n4.nabble.com/assign-td2291499.html
# Funktioniert als Zuordnung nur direkt im Skript, nicht in Funktionen.
# Eine Zuordnung zu bereits belegten Namen und Funktionen ist NICHT sinnvoll!
# Dazu gehören zB c, t, data, F, T, sin, ls. Data ist OK, weil großgeschrieben
# Außer bei FALSE, NA, etc (?reserved) ist es möglich, aber nicht empfehlenswert
# Wird die Zuordnung in Klammern gesetzt, wird sie gleichzeitig aufgerufen
(p8 <- 2854) # Das ist recht ungewöhnlich, man sollte es daher meiden.
?"<-" # ruft die Hilfe zur Zuordnung auf
# 1.4 Objektarten: -------------------------------------------------------------
# Abfragen des Datentyps:
mode(probe1) # numeric - Zahl
str(probe1) # STRucture # num
mode(p4) # character - Zeichenkette (Buchchstaben)
mode(p7) # numeric
str(p7) # int - Integer, ganze Zahl, Unterform von num (-> 1.5)
# Abfrage, welche Zuordnungen bestehen
ls()
ls.str()
# Zuordnungen löschen (Objekt aus dem Workspace löschen)
rm(p7)
p7 # Fehler: Objekt 'p7' nicht gefunden
rm(list=ls()) # entfernt alle Objekte, räumt also auf. Im TinnR F11 drücken.
p9 <- "1.3" # passiert manchmal beim Einlesen von Dateien (Kap 3)
p9
2*p9 # geht nicht
as.numeric(p9) # jetzt ist aus dem Text eine Zahl geworden
2*as.numeric(p9) # und kann damit entsprechend gerechnet werden
p9 # das Objekt selbst wurde dabei nicht verändert
is.numeric(p9) # is.___ fragt ab, ob ein Objekt ein bestimmter Datentyp ist
as.character(6:8) # as.___ erzwingt eine Objektart-änderung.
as.numeric(TRUE)
as.logical(0:1) # FALSE = 0, TRUE = 1
# 1.5 Grund-Datentypen: --------------------------------------------------------
?numeric # Zahlen
? integer # ganze Zahlen / werden von mode als numeric ausgegeben,
? double # KommaZahlen, NaN, Inf / str unterscheidet genauer
?logical # Wahrheitswerte: T=TRUE, F=FALSE
?character # Zeichenketten wie "Ueberschrift"
?complex # Komplexe Zahlen wie 7+3i
?factor # Kategorie, kategorielle Variablen
# 1.6 Objekttypen --------------------------------------------------------------
?vector # Vektor
?data.frame # gängigste Form einer Tabelle: verschiedene Datentypen in Spalten
?matrix # Jedes Element gleichen Datentyps! zweidimensionaler array
?array # mehrdimensionale Matrix (zB für räumliche Daten über Zeit)
?list # mehrere Modi in einem Objekt möglich, Indzierung mit [[n]]
?ts # Zeitreihe
# und einige mehr, die hier noch nicht wichtig sind.
mode(p9) # Daten- oder Objekt-typ, einiges ist dabei zusammengefasst
typeof(p9) # R-intern, vom Nutzer fast nie verwendet, genauer als mode
str(p9) # Gesamtstruktur vor allem komplizierter Objekte
class(p9) # Objekttyp
class(7:8) # zeigt bei Vektoren den Datentyp an
is.vector(7:8)# für Abfragen, ob ein Objekttyp "vector" ist.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 2: Vektoren und Tabellen ###########################################
# 2.1 Vektoren erstellen + Indizieren
# 2.2 rep, seq, :
# 2.3 Vektoren mit character
# 2.4 Basic Statistics
# 2.5 oft verwendete Operatoren
# 2.6 Tabellen
# 2.7 Tabellenindizierung
# 2.8 Tabellen zusammenführen
# 2.9 Zeilen und Spaltennamen
# 2.10 Matrizen
# 2.1 Vektoren erstellen: c + Indizieren: [] -----------------------------------
# Vektoren erstellen, Elemente durch Kommas trennen
vekt1 <- c(1 , -2 , 3, 7:11) # c steht für combine oder concatenate
vekt1
mode(vekt1)
# Auf Vektorelement zugreifen mit eckigen Klammern
vekt1 [5]
vekt1[5:6]
# Vektor mit einer entfernten Stelle erstellen/anzeigen: [-Element]
vekt2 <- vekt1 [-5]
vekt2
# Manipulation einzelner Elemente
vekt2[5] <- 28
vekt2
# Hinzufügen von Elementen. Zwischenkomponenten werden als NA eingefügt.
vekt2[10] <- -13
vekt2
vekt2[7]*3 # Zeigt nur das Ergebnis einer Berechnung an, ändert nicht den Vektor
vekt2
( k <- c(2,4,5) ) # Zuordnung in Klammern zeigt entstehendes Objekt gleich an
vekt2[k] # zugegriffen wird immer vektoriell
# 2.2 rep, seq, : --------------------------------------------------------------
# Mehrfaches Wiederholen einer Zahl im Vektor: rep(Zahl, AnzahlWiederholungen)
vekt3 <- rep(2, 10) # repeat
vekt3
# Sequenz: seq(von, bis, Schrittweite)
( vekt4 <- seq(1, 5, 0.5) )
# Sequenz: seq(von, bis, len = Anzahl Elemente)
( vekt5 <- seq(1, 3, len=15) ) # len kurz für length.out
seq(7, 6, -1/3) # Bei absteigenden Sequenzen negativen "by"-wert nutzen!
8:13 # Ganze Zahlen (integer) von 8 bis 13
-3:6
# folglich beim Entfernen beim Indizieren Klammern verwenden:
vekt1
vekt1[-(3:6)] # und nicht etwa so:
vekt1[-3:6] # Fehler: nur Nullen dürfen mit negativen Indizes gemischt werden
# 2.3 Vektoren mit character ---------------------------------------------------
# Zahlen als Text (R kann damit nicht rechnen, ist als Überschrift sinnvoll)
( vekt6 <- as.character(1:3) )
# Zusammengesetzte Vektoren
( vekt7 <- c(vekt1, -7:10, vekt2) )
# Namen zuordnen
( vekt8 <- c(2:4) )
vekt9 <- c("a","b","c")
vekt9
names(vekt8) <- vekt9
vekt8 # ist nach wie vor ein Vektor, und keine Tabelle!
vekt8["b"]
seq(7, 8, .3)
rep(1:4, 2)
rep(1:4, each=2)
paste("WunderbareSache", 3:5) # macht einen Vektor mit Modus character
paste(2, "hot", 4, "U")
paste("SeeYouL", vekt7[5], "ter", sep="") # Auf bestehende Objekte zugreifen
# sep steht für Separation, dem Trennzeichen zwischen den Zeichenketten.
# 2.4 Basic Statistics ---------------------------------------------------------
Noten <- c(3.0, 2.3, 5, 1.3, 4, 1.7)
Noten # Objekt anzeigen lassen
mode(Noten) # numeric also Zahlen
mean(Noten) # Mittelwert
var(Noten) # Varianz
sd(Noten) # Standardabweichung
median(Noten)# Median: Zahl an mittlerer Stelle der nach Größe sortierten Werten
# 2.5 oft verwendete Operatoren ------------------------------------------------
length(vekt1)# Anzahl Komponenten (Elemente) eines Vektors
length(vekt6)
max(vekt1) # Wert des größten Elements in vekt1
max(vekt2)
max(vekt2, na.rm=T) # NA Werte nicht berücksichtigen (steht für NA.ReMove)
vekt2 <- replace(vekt2, 8:9, 5) # ersetzt die Elemente 8:9 zum Wert 5
vekt2
max(vekt1, vekt2)
range(x) # entspricht c(min(x), max(x))
which.max(x) # Index des größten Elements in x
which.min(x) # Index des kleinsten Elements in x
sum(vekt1) # Summe der einzelnen Elemente
vekt2
sum(vekt2)
cumsum(vekt2) # kumulierte Summe bis zum jeweiligen Element
sum(vekt1, vekt2) # sum nimmt mehrere Daten, ohne dass diese erst
# zu einem Vektor zusammengefügt werden müssen. Bei mean ist das zB anders:
mean(vekt1, vekt2) # geht nicht
mean( c(vekt1, vekt2) )
# prod() # Produkt
# union() # Teilmenge
# intersect() # Schnittmenge
# rev() # reverse, kehrt Reihenfolge um
# sort() # sortiert Vektoren
# Ist ein bestimmter Wert (5) in einem Datensatz (y) enthalten?
y <- 3:17
is.element(5, y) # identisch zu x %in% y
2 %in% y
# ><= Abfragen
( k <- c(1,5,3,2,7,6,1,2,8) )
which(k==2) # Welche Elemente sind = 2 ?
k[k==2] <- 2.5 # ändere alle 2 zu 2.5
k
subset(k, k>3) # Welche Werte haben die Elemente mit Wert > 3 ?
# 2.6 Tabellen -----------------------------------------------------------------
# Dataframes = Tabellen mit verschiedenen Datentypen
# data.frame( Spaltenname = Vektor, zweiteSpalteÜberschrift = andererVektor)
tab1 <- data.frame(N=11:17, K=letters[1:7], Logisch=(1:7)>2)
tab1 # data.frame ist das gängigste Tabellenformat.
# Anzahl Reihen(Zeilen) / Spalten im Datensatz
nrow(tab1)
ncol(tab1)
dim(tab1)
# häufige Fehlermeldung: "Fehler in ____nrow(xx)____ : Argument hat Länge 0 "
# evtl ist xx ein Vektor --> nehme length(xx) statt nrow(xx) / dim(xx)
# Namen der einzelnen Datenspalten aufrufen
names(tab1)
str(tab1) # STRucture: Zeigt Objekttyp sowie Datentyp einzelner Spalten an
class(tab1) # data.frame
mode(tab1) # Obertyp list, hier erstmal irrelevant
# 2.7 Tabellenindizierung ------------------------------------------------------
# Auf Elemente der Tabelle zugreifen: Objekt[Zeile,Spalte]
tab1[2,3] # Ausgabe des zweiten Elementes der dritten Spalte
# Beim weglassen der Angabe von Zeile oder Spalte werden alle genommen:
tab1[ ,2:3] # Ausgabe einzelner Spalten als Spalte (Dataframe)
tab1[3,] # Ausgabe einer Zeile
tab1["K"] # Ausgabe einer Spalte als Dataframe
tab1 [3] # andere Methode für's gleiche
tab1[,3] # Ausgabe der Spalte als Vektor
tab1$K # dito
tab1[1:4,3] # Ausgabe eines Teils der Spalte als Vektor
tab1[1:4,3, drop=F] # Bleibt ein data.frame
# Noch mehr Beispiele in Kapitel 3.3
as.character(tab1) # geht manchmal nicht für Tabellen oder macht Vektor draus
# Zeilen entfernen
d <- data.frame( a=1:10 , b=11:20 )
d
( e <- d[-(3:5),] )
# Wertebedingt indizieren
blub <- data.frame(w=c("a001","b476","b342","a677"), A=c(4,5,3,2),
B=c(2,4,3,6), C=c(2,8,4,1) )
blub
blub[,2] <5
blub[blub$A<5,]
blub[blub$A<5 & blub$B<5,]
blub[blub[,2]<5 & blub[,3]<5,]
# 2.8 Tabellen zusammenführen --------------------------------------------------
# Mit rbind(a,b,d) lassen sich die Vektoren a, b und d reihenweise zu einer
# Matrix zusammenführen.
# Im Gegensatz dazu verbindet cbind() die Vektoren spaltenweise zu einer Matrix.
# data.frame() macht das auch so wie cbind.
f <- data.frame(a=11:14, b=58:55, d=236:233)
g <- data.frame(a=33:36, b=76:79, d=-28:-31)
f
g
rbind(f,g) # Spaltenanzahl ncol der beiden Matrizen muss gleich sein
cbind(f,g) # Zeilenanzahl nrow der beiden Matrizen muss gleich sein
# 2.9 Zeilen- und Spaltennamen zuordnen ----------------------------------------
a <- data.frame( 5:-1, c(3,5,1,9,6,4,8) ) ; a
colnames(a) <- c("alles","klar") ; a
rownames(a) <- LETTERS[a$klar] ; a
# 2.10 Matrizen ----------------------------------------------------------------
# matrix ist ein weiterer Objekttyp. Hat anders als Dataframe nur 1 Datentyp!
# Matrix erstellen und füllen
(m <- matrix(NA, nrow=5, ncol=3))
m[4,2] <- 17 ; m
m[1,] <- 8 ; m
m[,3] <- 3:7 ; m
( m[5,] <- LETTERS[4:6] )
m # matrix kann nur 1 Modus haben, darum wird jetzt alles in char umgewandelt!
noquote(m) # Gibt die Character ohne Anführungsstriche aus
str(m); class(m); typeof(m) # bleibt eine matrix mit charactern
# Lineare Gleichungssysteme lösen x1 x2 b
A <- matrix(c(16,4,74,16), ncol=2, byrow=T) ; A # 16 4 | 25
b <- c(25,113) # 74 16 | 113
solve(A, b) # 1.30 1.05
1.3 * 16 + 1.05 * 4 # 25
1.3 * 74 + 1.05 * 16 # 113 Die Resultate sind richtig...
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 3: Dateien und Arbeitsverzeichnis ##################################
# 3.1 Arbeitsverzeichnis
# 3.2 Tabellen einlesen
# 3.3 Anzeigen ausgewählter Elemente
# 3.4 attach
# 3.5 Tabelle exportieren (schreiben)
# 3.6 Oft verwendeter Operatoren
# 3.7 Excel-Dateien lesen und schreiben
# 3.1 Arbeitsverzeichnis (Working Directory) -----------------------------------
# An diesem Ort (Ordner) werden Dateien gelesen und gespeichert, sofern nicht
# der ganze Pfad angegeben wird. Dadurch kann ein Ordner mit Dateien und Skript
# einem Kollegen gegeben werden, ohne das der alle Pfade ändern muss.
# Schnelles und reproduzierbares Einstellen mit dem Befehl setwd:
setwd("D:/Ordner/BerryIstToll/Rclick_Handbuch")
# alle \ umwandeln in / oder \\. Bei langen Pfaden im TinnR STRG + R nutzen!
# ( \ ist in R die EscapeSequenz aus einer Zeichenkette heraus)
# Tipp: Im Windows Explorer Datei markieren, F4 drücken um den Pfad zu kopieren.
# Tipp: Wähle sinnvollerweise gut zugängliche Ordner aus.
# Alternativ per Maus, aber anderen Skriptlesern dann nicht ersichtlich:
# In der R Console "Datei" klicken, dann "Verzeichnis wechseln".
# Überprüfen des aktuellen Verzeichnisses mit
getwd()
# Wenn ein R-Skript in Tinn-R geöffnet wird, und danach R selbst mit Tinn-R
# geöffnet wird, ist das Arbeitsverzeichnis automatisch der Ordner, in dem die
# R-Datei sich befindet. Wenn also die Daten im gleichen Ordner wie das R-Skript
# vorliegen, entfällt die meiste Arbeit mit dem Working Directory.
# Anzeigen der Datei(en) im Arbeitsverzeichnis
dir() # Übersicht der Dateien und Ordner (keine Unterordner)
dir("S:/Dropbox/Public/Rclick", recursive=T) # mit Dateien in Unterordnern
file.info(dir()) # Detaillierte Übersicht
file.show("Daten/3_Einlesen.txt" ) # eine Datei anzeigen
# Dateinamen sind, anders als Zuordnungen, nicht case-sensitive (GROSS/klein).
getwd()
setwd("..") # 1 Ordner aufwärts. mit setwd("/") zur Partition, z.B. "D:/"
getwd()
# 3.2 Tabellen einlesen --------------------------------------------------------
MeiTab <- read.table("Daten/3_Einlesen.txt", header = T)
# Der Dateiname kann auch als ganzer Pfad angegeben werden, bei Weitergabe von
# Skripten und Daten sind aber relative Pfade ganz toll!
# read.table wichtigsten Argumente
# header = TRUE übernimmt Spalten-Überschriften
# dec = "," liest Kommas als Dezimaltrennstellen
# sep = "-" Nimmt Striche als Spaltentrennung an
# fill = T füllt Lücken in den Daten mit NA
# (bei fill=F (Default!) müssen alle Zeilen und Spalten gleich lang sein)
# skip = 12 überspringt die ersten 12 Zeilen (zB mit Metadaten)
# comment.char = "%" überspringt Zeilen, die mit % anfangen (wie # in R)
# read.csv() für csv-tabellen und read.csv2() für deutsche Tabellen.
# Excel-Tabellen möglichst als txt-Dateien exportieren ("Speichern unter")
?stackloss # Wärme-Abgasverlust - Effizienz industrieller Ammoniumoxidation
# Spaltennamen anzeigen
names(MeiTab)
# Grundlegende statistische Auswertung (Zusammenfassung) der Daten
summary(MeiTab)
# häufige Fehlermeldung beim Laden von Tabellen
# "Kann Datei ___ nicht öffnen: No such file or directory"
# --> Verzeichnis oder Pfad kontrollieren: getwd()
# --> Dateiname kontrollieren
# --> Wird die Datei mit Anführungsstrichen aufgerufen?
# Wer viele Excel-Dateien einlesen muss: Ganz unten Abschnitt 3.7 lesen
# 3.3 Anzeigen ausgewählter Elemente -------------------------------------------
# (können durch Zuordnung auch als neue Tabelle/Vektor gespeichert werden)
MeiTab["Air.Flow"] # Spalte "Air.Flow" als Spalte eines Dataframes
MeiTab[1] # Spalte 1 als Spalte
MeiTab$Air.Flow # Spalte "Air.Flow" als Vektor
MeiTab[,"Air.Flow"] # Spalte "Air.Flow" als Vektor
MeiTab[,1] # Spalte 1 als Vektor
MeiTab[1:3] # Spalten 1 bis 3
MeiTab[3:7,2:3] # Teil der Tabelle
MeiTab[1,] # Zeile 1 als Tabelle
MeiTab[1:3,] # Zeilen 1 bis 3
MeiTab[ 6:12 , 2:3 ] # Zeilen 6 bis 12, Spalten 2 und 3
MeiTab[1,"Air.Flow"] # Element der Zeile 1 in der Spalte "Air.Flow"
MeiTab[ 1:3 , 4 ] # Zeilen 1 bis 3 in Spalte 4 als Vektor
MeiTab[MeiTab$Air.Flow<60,] # Tabellenzeilen mit Air.Flow < 50
MeiTab[MeiTab$Air.Flow>=70,] # Tabellenzeilen mit Air.Flow grösser/gleich 70
MeiTab[MeiTab$Air.Flow==62,] # Tabellenzeilen mit Air.Flow = 62
MeiTab[MeiTab$Air.Flow!=62,] # Tabellenzeilen mit Air.Flow ungleich 62
MeiTab$Water.Temp[MeiTab$Air.Flow>70] # Water.Temp mit Air.Flow > 70
# Ausgaben als Spalte (Data.Frame) sind übersichtlich
# Mit Vektoren kann aber oft bessser gerechnet werden
mode(MeiTab[1]) ; mode(MeiTab[,1])
# häufige Fehlermeldung: Fehler in "__": undefined colums selected.
# --> beachte das Komma! Objekt[zeilen , spalten]
# --> dim(Objekt) zeigt dir, ob es überhaupt soviele Spalten hat,
# wie du ansprechen willst.
# 3.4 schlechtes Programmieren: attach -----------------------------------------
# Anhängen einer Tabelle
# mit attach() wird ein data.frame in den R Suchpfad eingebunden
# die einzelnen Spalten können jetzt mit ihrem Namen direkt angesprochen werden.
# Sind aber keine Objekte im normalen Environment (Vgl ls() vor und nach attach)
ls()
Air.Flow # noch nicht bekannt
attach(MeiTab)
Air.Flow # jetzt direkt abrufbar
ls()
log(Air.Flow) # als Zeile (Vektor, "numeric")
log(MeiTab["Air.Flow"]) # als Spalte (Tabelle, "list")
MeiTab[Water.Temp>22,]
# attach ist generell eine Unsitte, weil man später nicht weiß, woher ein objekt
# kommt: aus dem normalen Workspace oder aus einer Dataframe.
# Und es kann Objekte und Funktionen im Base "überschreiben".
search() # zeigt Suchpfad, wo die angehängten Sachen sind
# Abhängen einer Tabelle (Die Tabelle selbst ist weiterhin im Workspace)
detach(MeiTab)
# 3.5 Tabelle exportieren (schreiben) ------------------------------------------
# Wenn kein Pfad angegeben, wird im working directory abgespeichert
write.table(tab1, "MeineErsteTabelle.txt", row.names=F, sep="\t", quote=F)
# row.names = FALSE : Abspeichern ohne Zeilennummern links
# col.names = FALSE : Spaltenüberschriften nicht mitschreiben
# quote = FALSE : Keine Anführungsstriche in die Textdatei schreiben
# sep = "\t" : Daten tabstopp-getrennt speichern (ab 7 Zeichen pro Spalte
# ist die Darstellung im Editor verschoben). Zum Copypasten in Excel hinein.
# Escape-Sequenz "\" im Character belegt, Siehe Kap 17.1
# Kann daher für Verzeichnisse nicht genutzt werden (siehe oben unter 3.1).
# 3.6 Oft verwendeter Operatoren (to be continued) -----------------------------
MeiTab
order(MeiTab$Water.Temp) # Reihenfolge
MeiTab[ order(MeiTab$Water.Temp) , ] # Tabelle sortiert nach Wassertemperatur
nrow(MeiTab) # Anzahl Zeilen (rows)
ncol(MeiTab) # Anzahl Spalten (columns)
dim(MeiTab) # entspricht c(nrow(),ncol()), macht das aber auch für arrays
t(MeiTab) # transponieren: Zeilen zu Spalten, Spalten zu Zeilen
?split() # Vektoren zerlegen
data() # bestehende (R- oder Package-integrierte) Datensätze anzeigen
# Und für Matrizenoperationen:
?det() # Determinante einer Matrix
?solve(A,b) # Lösung des Gleichungssystems A x = b
?crossprod(X,Y) # Kreuzprodukt zweier Matrizen
?qr() # QR-Matrixzerlegung
?eigen() # Eigenwerte und -Vektoren
?"%*%" # Produkt zweier Matrizen / Skalarprodukt zweier Vektoren
# 3.7 Excel-Dateien lesen und schreiben ----------------------------------------
# http://cran.r-project.org/doc/manuals/R-data.html#Reading-Excel-spreadsheets
# Für Excel 2007/2010/2013 siehe besonders das Paket xlsx
# Einzelne Dateien lassen sich schneller als Text- oder CSV-dateien speichern.
# Zunächst Java installieren, falls nicht vorhanden:
browseURL("http://www.java.com/en/download/manual.jsp")
# ASK Toolbar beim Installieren ausschalten
# Beim Fehler JAVA_HOME cannot be determined from the Registry:
Sys.setenv(JAVA_HOME='C:\\Program Files\\Java\\jre7') # for 64-bit version
Sys.setenv(JAVA_HOME='C:\\Program Files (x86)\\Java\\jre7') # for 32-bit version
library(rJava)
if(!require(xlsx)) install.packages("xlsx")
library(xlsx)
read.xlsx("filename.xlsx", sheetIndex=2)
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 4: Graphische Datenauswertung ######################################
# 4.1 Basics: Beschriftungen, Symbole
# 4.2 Linien statt Punkte, Überplotten
# 4.3 Graphikausschnitte, Ränder, Barplots
# 4.4 par
# 4.5 Farben, Legende
# 4.6 Linien, Text, Gitter
# 4.7 Achsenbeschriftungen 2.0
# 4.8 Logarithmierte Achsen mit Linien im Plot
# 4.9 Anordnung mehrerer Graphiken in einem Bild
# 4.10 Graphiken mit Umrandung
# 4.11 Speichern von Bildern: bmp("dateiname.bmp")
# 4.12 Beidseitige y-Achsenbeschriftung
# 4.13 Achsen unterbrechen, Pakete installieren und laden
# 4.14 mehrdimensionale Info in Graphiken
# 4.15 Dreidimensionales Plotten
# 4.16 Fehlerbalken (Error bars)
# 4.17 Fläche unter Graphen füllen
# 4.18 Funktionen zeichnen
# 4.19 zooming in graphics
# 4.20 und ein würdiger Abschluss des Kapitels
# 4.21 Übersicht der wichtigsten Graphikbefehle
# integrierter Datensatz (Beispieldaten aus Kapitel 3)
Daten <- stackloss ; names(Daten) <- c("Air", "Water", "Acid", "StackLoss")
Daten # Zuordnungen können unter x-beliebigem Namen laufen, eben zB "Daten".
?stackloss # aber nicht "data" weil das schon eine Funktion ist
# Grundlegende Graphiken (x-y-Plot)
?plot
plot(Daten$Air, Daten$StackLoss)
# manchmal ist auch pairs interessant:
pairs(Daten)
# Intro in R's Graphikkapazitäten (Befehle sind jetzt egal, später anschauen)
demo(graphics) # in der Console ENTER drücken, dann auf die Bilder klicken
# Evtl. kommt zuerst ein leeres Grafikfenster, dann einmal darin klicken.
dev.off() # Schließt Graphikfenster
#http://blog.revolutionanalytics.com/2015/01/some-basics-for-base-graphics.html
# 4.1 Basics: Beschriftungen, Symbole ------------------------------------------
# Beschriftungen: main = "Überschrift", sub = "Unterschrift", xlab = "X-Achse"
plot(Daten$Water, Daten$StackLoss, main = "StackLoss abh. der Wassertemperatur",
sub = "Analyse", xlab = "Wassertemperatur", ylab = "StackLoss")
# Immer wenn ein + in R erscheint, ist die Eingabe noch nicht abgeschlossen.
# später hinzufügen mit title()
plot( c(8,4,14,-3,4,15) , xlab = "")
# Wenn nur ein Vektor angegeben wird, wird über ganze Zahlen geplottet (Index)
title(main = "Hallo", xlab = "Dies ist die X-Achse")
# Zeilenumbruch im Titel: \n ohne Leerzeichen im Text einfügen
plot(Daten$Water, Daten$StackLoss, main = "tolle Graphik\nStackLoss")
# Punkt- und Schriftgröße verändern: cex = ___ ( Character EXpansion)
plot(Daten$Water, Daten$StackLoss, cex = 2)
plot(Daten$Water, Daten$StackLoss, cex = 2, main = "Schuhgröße - StackLoss",
cex.main = 2) # Default: cex.main=1.2
# Achsenwerte ("label" = Zahlen / Einheiten an Achsen) skalieren: cex.axis = ___
plot(Daten$Water, Daten$StackLoss, cex.axis = 0.5)
# Ausrichtung der Achsenwerte: las = ___ 0 parallel zur Achse, 1 horizontal
plot(Daten$StackLoss, las = 1) # 2 quer zur Achse, 3 vertikal
# (Label Axis Style)
# Symbole ändern: pch = ___ oder pch="Buchstabe" zB. ausgefüllter Punkt: pch=16
plot(Daten$Water, pch = 16) # (Point Character) Siehe Übersicht im Anhang
# Achsentitel ausrichten (adjust); Werte von 0 (Text links) bis 1(Text rechts)
plot(Daten$Acid, adj = 0.82)
# Kombination von Parameter-settings (Reihenfolge ist egal)
plot(Daten$Air, main = "N in Pflanzen", xlab = "WasserTemp",
ylab = "Stack Loss", las=3, pch = 16, cex = 1.5)
# Abstand Beschriftung: mgp=c( Achsentitel, Achsenwerte, Achsen)
# Voreinstellung: mgp=c(3,1,0), somit Verschiebung nach innen und außen möglich
plot(Daten$Water, mgp=c(3.5, 1, 0))
plot(Daten$Water, mgp=c(2, 1, 0)) # könnte für MarGinPlacement stehen...
plot(Daten$Water, mgp=c(3, -1, 0))
plot(Daten$Water, mgp=c(3, -3,-2))
# 4.2 Linien statt Punkte, Überplotten -----------------------------------------
k <- c(3,0,-10,8,6,3) ; l <- c(-4,2,-6,-12,7,-3)
plot(k) # standard (default) type: "p" (points)
plot(k, type="l") # Linien
plot(k, type="b", lwd=4) # Liniendicke einstellen
plot(k, type="o", lty=2) # Linientyp einstellen
# Siehe Übersichten im Anhang zu type und lty
# mit dem Befehl points() kann etwas in die Graphik eingefügt werden.
# Die Argumente sind gleich wie bei plot(), allerdings ohne Titel.
plot(k)
points(l, pch=24, bg=2)
lines(l) # wie points mit als Standard type="l"
plot(k, type="o") # "b": Punkte und Linien. "o" zieht die Linien durch
par(new=T) # Ein danach versandter plot wird über das alte geplottet.
plot(l, type="o", col=6, col.axis=6) # Achsen vermurkst!
# (Zweiseitige Achsen werden im Kapitelteil 4.12 vorgestellt.)
plot(k, type="b", lty=2) # Plot von gerade wieder (new=T gilt nur für 1 Graphik)
par(new=T)
plot(l, type="o", ylim=c(-10,8), pch="@")
# Die y-Achsengrenze des ersten Datensatzes reicht nicht, daher ylim: Kap 4.3.
# 4.3 Graphikausschnitte, Ränder, Barplots -------------------------------------
# Plotausschnitt einstellen: xlim = c( untere Grenze, obere Grenze )
plot(Daten$StackLoss,xlim=c(-3,50), ylim=c(-4,80))
# Graphik ohne 4%-Erweiterung des Achsenbereiches: yaxs="i" (Y-AXis Style)
plot(c(3,2,5,4), type="o")
plot(c(3,2,5,4), type="o", yaxs="i", xaxs="i")
# Ränder (Margins) einstellen: mar = c(bottom, left, top, right)
# Die Nutzung von par ist im nächsten Abschnitt erklärt.
par(mar = c(5, 4, 4, 2)+0.1 ) # Default Settings (Voreinstellung)
plot(Daten$Water) # mar muss immer mit par für das Device gesetzt werden
par(mar=c(2.2, 2.0, 0.2, 2.1) )
plot(Daten$Water)
# Barplots:
tabelle <- read.table(header=TRUE, text="
LandN. Verdunstung
Wald 200
Weide 400
Soja 359")
tabelle
barplot(tabelle$Verd, names=tabelle$LandN, main="Verdunstung",
col=c("forestgreen","green2","wheat"))
# 4.4 par ----------------------------------------------------------------------
# viele der Graphik-Einstellungen können auch mit par eingestellt werden.
# Die danach im gleichen Fenster erstellten Plots haben dann die Eigenschaften.
# Weitere Fenster können mit windows() eröffnet werden.
windows() # geht nur mit dem Betriebssystem Windows, sorry Linuxer.
# (Macser: selbst Schuld 😉 Kauft was billigeres oder kostenloses...)
# verschiebe mal das neue Fenster (Device 3)
plot(Daten$Air, Daten$Water)
dev.off() # schließt das aktuelle (active) Fenster
graphics.off() # macht alle zu
# Aktuelle Settings anzeigen (Default, da alle Graphikfenster jetzt zu sind)
par()
?par
# Alte Graphikeinstellungen (Parameter) im Graphics:Device speichern
altepar <- par(mar=rep(6,4)) # speichert das alte (!) mar
altepar # in Funktionen oft old.par oder op genannt
plot(1:8)
par(altepar) # vorherigen Rand wieder einstellen
par(new=T, col="red", col.axis="red", col.lab="red") # Farben Siehe 4.5
plot(1:8)
# Alternativ möglich:
old.par <- par(no.readonly=T)
# Nur nötig, wenn Änderungen im Device mit par(___) vorgenommen werden
# Nicht nötig, wenn Einstellungen in einzelnen Graphiken erfolgen. plot(x,y,___)
dev.off()
# Randplatz minimieren, Graphikraum maximal nutzen:
par(mar=c(3,3,1,1), mgp=c(1.8, 0.6, 0), las=1)
plot(Daten$Air, Daten$Water)
# Einstellungen mit par gelten für Graphen, die im aktuellen, aktiven, Fenster
# erstellt werden. Wird dieses geschlossen, müssen sie erneut gesetzt werden.
# Graphik History (zeichnet einzelne Ausgaben im Graphics:Device auf)
# Bild wechseln mit "Bild up" und "Bild dn" Tasten oder Mausklick
windows(record=T) # oder im Device selbst unter History anklicken
# praktisch: Rechtsklick in der Graphik - "bleibe im Vordergrund" anschalten
graphics.off()
# 4.5 Farben, Legende ----------------------------------------------------------
# Farben setzen: col = ___
plot(Daten$Water, Daten$StackLoss, col = "red") # oder col=2
plot(Daten$Water, Daten$StackLoss, col = rgb(0,1,0), pch=16)
# RGB steht für RedGreenBlue. Ermöglicht exakte Farben zB nach Corporate Design
rgb(0,1,0) # Hat Standardwerte zwischen 0 und 1
col2rgb("green")/255 # Hat Standardwerte zwischen 0 und 255
plot(1:3, pch=16, cex=3, col= rgb(c(0,.3,.8), c(.3,.5,.5), c(.1,.9,.2) ), las=1)
# Siehe auch ?hcv und ?hsl
plot(Daten$StackLoss,Daten$Air, col.axis=6, font=2)
# font=2 ist fette Schrift (bold). font.lab ist für die Achsentitel
colors() # Siehe Anhang und http://fany.checkts.net/Farben.html
# Die acht Schnellfarben, die standardmäßig in palette() sind:
plot(1:8, col=1:8, pch=16, cex=4, xlab="plot ( 1:8, col=1:8 )")
# Hintergrundfarbe einstellen
par(bg=3) # BackGround
plot(1,1, pch=16)
dev.off()
# Achsenfarben setzen: par(col=__)
old.par <- par(col=2, cex.axis=2, bg=8)
plot(1,1)
plot(Daten$Water, Daten$Acid)
# alte Grafikeistellungen wiederherstellen
par(old.par)
plot(Daten$Water, Daten$Acid)
# Alternativ das Graphikfenster schließen (geht auch normal übers "Kreuzigen"):
dev.off() # Beim nächsten Plotten sind die Settings wieder im Ursprungsformat
# Farben der Achsenwerte (Einheiten) : col.axis
# und der Beschriftungen (Labels): col.lab
plot(1, col.axis = "sky blue", col.lab = "brown")
# Symbole mit Farben füllen
plot(1,1, pch= 1, cex=4) # (Point CHaracter, Character EXpansion)
plot(1,1, pch=21, cex=4) # pch 21 bis 25 können eingefärbt werden.
plot(1,1, pch=21, cex=4, bg=3) # --> Übersicht der Symbole im Anhang
# bg hat hier also eine andere Bedeutung, als wenn es in par() gesetzt wird.
# bedingt einfärben: col = ifelse(wenn, dann, sonst)
plot(Daten$Water, Daten$StackLoss,
col = ifelse(Daten$StackLoss > 18, "blue", "red"), pch=15)
plot(Daten$Water, Daten$StackLoss,
col = ifelse(Daten$Acid > 87, "blue", "red"), pch=15)
# Legende hinzufügen: legend(x,y-position der linksoberen Ecke, Einträge, ___)
legend(18, 35, c("Acid > 87", "Acid <= 87"), col=c(4,2), pch=15)
plot( c(1,3,2,5), type="l") # Eine Übersicht der TYPEs ist im Anhang.
legend("topleft", c("Hallo", "R-Lerner"), lty=1) # Line TYpe siehe Anhang
legend(3, 4, c("CheckThisOut",1,2), pch=1:3, bg="white")
legend(1.5, 3, c("Party","People"), title="Legende", bty="n", fill=3:4)
legend("bottom", c("Und", "tschüss"), inset=0.01, lty=c(NA, 2), pch=c(3,NA))
# Farbenrad aus der demo(graphics)
pie(rep(1,24), col=rainbow(24), radius=0.9)
k <- c(3,1,10,8,6,3)
pie(k) # PieCharts sind fast nie geeignet, Daten adäquat zu visualisieren:
browseURL("www.storytellingwithdata.com/2011/07/death-to-pie-charts.html")
browseURL("http://www.perceptualedge.com/articles/08-21-07.pdf")
browseURL("http://www.stevefenton.co.uk/Content/Pie-Charts-Are-Bad/")
browseURL("http://blog.jgc.org/2009/08/please-dont-use-pie-charts.html")
browseURL("http://www.juiceanalytics.com/writing/the-problem-with-pie-charts/")
browseURL("http://eagereyes.org/techniques/pie-charts")
# 4.6 Linien, Text, Gitter -----------------------------------------------------
# Linie(n) einzeichnen: abline(a=Y-Schnittpunkt, b=Steigung, h=Y-wert, v=X-wert)
plot(Daten$Water, Daten$StackLoss)
abline(h=20) # horizontale, waagerechte Linie
abline(v=24) # vertikale, senkrechte, Linie
abline(a=-38, b=3) # schräge Gerade: a = Y-Schnittpunkt, b = Steigung (!)
abline(5, 0.7, col=6, lty=2, lwd=2) # Gerade ist also y = a + bx
# Graphikfenster öffnen, ohne zu zeichnen: type="n" für none/nothing/nichts/nada
# Verhältnis X zu Y (ASPect ratio): asp = ___ (je gleich groß: asp=1)
plot(c(-2,3), c(-1,5), type="n", xlab="x", ylab="y", asp=1)
abline(h=0, v=0, col="gray60")
# Text hinzufügen: text(x,y-Koordinaten, "text", Ausrichtung vom Punkt aus
text(1, 0, "abline( h = 0 )", col = "gray60", adj = c(0, -.1)) # adj: s. Anhang
text(1, 2, "fette Schrift", col = "gray60", font=2) # font: Übersicht im Anhang
# Schriftart verwenden:
windowsFonts(MyFont = windowsFont("Georgia"))
windowsFonts() # Für alle Bilder: am Skriptanfang par(family="MyFont")
plot(1:100, type="l") # einzeln: mit family in plot, text, mtext, etc
plot(1:100, type="l", family="sans") # Arial, bei mir Default
plot(1:100, type="l", family="MyFont") # Georgia
plot(1:100, type="l", family="mono") # Courier New
plot(1:100, type="l", family="serif") # Times New Roman
windowsFonts(Bookman = windowsFont("Bookman Old Style"),
Comic = windowsFont("Comic Sans MS"), Symbol = windowsFont("Symbol"))
plot(1:100, type="l", family="Bookman")
plot(1:100, type="l", family="Comic")
plot(1:100, type="l", family="Symbol")
# Linientyp einstellen: lty = ___
abline(h = -1:5, v = -2:3, col = "lightgray", lty=3)
abline(a=1, b=2, col = 2)
text(1, 3, "abline( 1, 2 )", col=2, adj=c(-.1,-.1))
plot(8:5, pch=22, cex=3, bg=5)
grid(col=2) # überzeichnet die Ergebnisse, diese können neu gelegt werden
points(8:5, pch=22, cex=3, bg=5)
# LinienEnd-typen einstellen (geht heutzutage auch im Plot-Befehl)
plot(c(5,5), lwd=8, type="l")
par(lend="butt") # 0="round" (Default), 1="butt", 2="square"
points(c(4,4), lwd=8, type="l")
# 4.7 Achsenbeschriftungen 2.0 -------------------------------------------------
# Achsenbeschriftungen hinzufügen: axis(welche, wo, was)
# welche: 1=unten, 2=links, 3=oben, 4=rechts
plot(67:88, yaxt="n", xaxt="n")
axis(2, 69 ) # einzelne Werte
axis(2, 76, "Grenze", las=1) # mit Text (oder anderen Zahlen)
axis(2, 83, "Gr2", las=1, col.axis=2, col=8) # schön farbig
axis(2, 86, tick=F, las=1) # ohne Striche
axis(3, 5, labels=F) # ohne Werte
axis(1, 5:7, col=5, col.axis=2) # Intervallachse bei mehreren Positionen
axis(1, pretty(10:30) ) # pretty wählt schöne Zahlen
bm <- prettyNum(c(10000, 20000, 30000), big.mark="'"); bm # Tausend-Trennpunkte
axis(4, c(70,75,80), bm)
axis(3, 13:16, col="orange", lwd=0, lwd.ticks=3) # no line
axis(3, pos=par("usr")[3], col.axis=4, font.axis=2) # Striche nach Innen
axis(2, 78, tcl=0.5, mgp=c(0,-1.8,0), las=1) # Striche nach Innen (besser)
# die letzte Zahl von mgp ist hier für die Ticks (= Striche) und nicht
# für die Achsen selbst, wie im Abschnitt Basics
axis(1,15.5, "", tcl=-0.4, col.ticks="purple")
plot(5:9, xlab=expression(Das~underline(sind)~italic(meine)~bold(Daten)~hier) )
?expression
?plotmath
?axTicks()
plot(c(5,-2,7,8), main=c("per Hand geteilte",""), col.main=2)
title(main=c("","Überschrift"), col.main=4, cex.main=.7)
# 4.8 Logarithmierte Achsen mit Linien im Plot ---------------------------------
# Beschriftungen:
10^(-3:4)
options("scipen"=4) # scientific penalty, siehe ?options
10^(-3:4)
options("scipen"=4)
format(10^(-3:4), big.mark="'", trim=TRUE )
options("scipen"=0)
format(10^(-3:4), big.mark="'" )
install.packages("berryFunctions")
require(berryFunctions)
plot(1, ylim=c(1e-3, 1e4), log="y", yaxt="n", yaxs="i")
logAxis(2)
lines(seq(0.5, 1.5, len=50), 10^runif(50, -3, 4), col=2)
str(logVals()) # vals: values, labs: formatted labels, all: all 10^x values for lines
# 4.9 Anordnung mehrerer Graphiken in einem Bild -------------------------------
# mehrere Graphiken pro Fenster (Layout-Aufteilung)
par(mfrow=c(3,2)) # 3 Zeilen, 2 Spalten - 6 Graphiken in einem Fenster
plot(stackloss[,1:2])
plot(stackloss[,2:3])
# 6 (verschiedene) plot-Befehle an R schicken
# Zurücksetzen auf eine Grafik pro Fenster
par(mfrow=c(1,1)) # oder schließen und mit plot(x, y, ___) neu öffnen
# flexiblere Anordnung mehrerer Graphiken im Graphics:Device
layout(matrix(1:6, 3, 2), widths=c(7, 5), heights=c(1,1,1))
plot(stackloss[,1:2]) # widths und heights werden jeweils anteilig gesetzt
plot(stackloss[,2:3])
?layout
# mehrere Graphiken auf einem Blatt: absolute Größen (in cm)
# zB. für quadratische plots. Grafikfenster vorher ausreichend groß ziehen!
layout(matrix(1:3, 1 ,3), widths=lcm(c(5, 5, 5)), heights= lcm(5))
# Layout anzeigen
m <- layout(matrix(1:6, 3, 2), widths=c(8,3), heights=c(2,4,1))
layout.show(m)
lay <- layout(matrix(c(1,1,2,2, 5,5,2,2, 3,3,4,4, 3,3,4,4), ncol=4, byrow=TRUE))
layout.show(lay)
lay <- layout(matrix(c(2,2, 2,2, 1,5, 3,5, 4,4, 4,6), ncol=2, byrow=TRUE) )
layout.show(lay)
dev.off()
# margin-text (Text in den Rändern)
layout(matrix(1:2, ncol=2), width=c(2.7,1))
plot(1)
mtext("rechts", side=4, col=4)
mtext("plot region", adj=1)
mtext("outer margin", adj=1, outer=TRUE, line=-1)
box("figure", col=2)
mtext("Figure", at=par("fig")[2], adj=1, outer=TRUE, line=-2)
# Graphiken übereinander zeichnen:
browseURL("http://www.statmethods.net/advgraphs/layout.html")
graphics.off() # etwaige alte Einstellungen löschen - neu beginnen
plot(rnorm(10))
par(fig=c(0,0.6, 0.6, 0.95), mar=rep(0,4), new=T)
plot(1, ann=F, axes=F)
pu <- par("usr") ; rect(pu[1],pu[3],pu[2],pu[4], col=3); rm(pu)
par(mar=c(3,4,1,1), new=T)
plot(rnorm(80))
# Aspect Ratio (potentiell wichtig)
x <- (1:50)*0.92 ; y <- sin((1:50)*0.92) # http://www.amazon.com/dp/0521762936
plot(x,y)
plot(x,y, asp=3)
# 4.10 Graphiken mit Umrandung -------------------------------------------------
par(oma=c(1,2,2,4)) # outer margins
lay <- layout(matrix(c(1,1,2,3, 4,4,4,4, 5,5,5, 6), ncol=4, byrow=TRUE) )
layout.show(lay)
# Wenn man mit Layout arbeitet, sollte man vorm ersten Plot par(cex=1) setzen!
plot(1:9)
par(cex=1)
plot(2:-5)
box(lwd=2, col=7)
box("figure", col=2)
box("inner", col=3)
box("outer", col=4, lwd=4)
graphics.off()
hist(rnorm(20)) # 20 Zufallszahlen aus der Normalverteilung. Mehr in Kap. 5
# Histogramm siehe Kapitel 6.
box(col = "blueviolet", lwd=3, lty=2)
# zB wichtig nach abline, da das den Rand der Graphik selbst überplottet.
win.graph(width=3, height=2) # Öffnet Graphikfenster bestimmter Größe (inches)
# 4.11 Speichern von Bildern: bmp("dateiname.bmp") -----------------------------
# automatisch im WorkingDirectory, ansonsten im Dateinamen den Pfad einfügen
png("MeinBild.png") # Defaultbildgröße: 480 x 480 Pixel
# danach die Graphikbefehle abschicken
plot(1:10)
dev.off() # Device auch wieder schließen!
getwd() # Schau dir das Bild hier im Explorer deines Computers an...
# praktisch unter Windows:
# Im Graphikfenster - Datei - Speichern als - png
plot(c(7,5,11,2,4), type="b") # aber nicht reproduzierbar!
# *.wmf file that can be copy-pasted into a word file:
win.metafile("mygraph.wmf") # funktioniert nicht immer gut!
plot(cumsum(rnorm(30)), type="l") # Und nur unter Windoof.
dev.off()
# Zum Publizieren geeignete Bildformate:
png("MeinBild.png", width=600, height=400, units="px")
# png ist ein gutes, kleines, genaues Format
# 600 x 400 ist ganz gut, um direkt in DIN A4 (Word-)Dokumente einzufügen.
plot(1:10)
dev.off()
png("MeinBild.png", width=4.5, height=5, units="in", pointsize=14, res=200)
# Für Publikation Schriftgröße sinnvoll setzen! http://civilstat.com/?p=557
plot(1:10)
dev.off()
# eps-äquivalente (Encapsulated Postscript, kann nachbearbeitet werden)
postscript("versuch.eps", onepage=TRUE, horizontal=FALSE)
plot(rnorm(10)) # wenn dir eps nichts sagt, einfach überspringen.
dev.off()
# Fortlaufende Nummerierung bei mehreren Bildern: %03d im Dateinamen
png("Bild%d.png") # %03d hängt führende Nullen (leading zeros) vor
par(bg=2)
plot( c(6,9,0,2,8), type="h", lwd=7 )
plot( 1:8 )
# danach:
dev.off() # sonst ist das letzte Bild noch leer!
# Verschiedene Speichermöglichkeiten:
dir.create("PNG_Test_Bilder")
oldwd <- setwd("PNG_Test_Bilder")
png("Test_A_%d.png") # ohne führende Nullen von 1:n durchnummeriert
for(i in c(1:8, 837, 2087, 2333)) plot(i) # i ist nicht im Dateinamen!
dev.off()
png("Test_B_%4d.png") # mit Leerzeichen
for(i in 1:13) plot(i)
dev.off()
png("Test_C_%04d.png") # mit Nullen
for(i in 56:68) plot(i)
dev.off()
for(i in c(1:8, 837, 2087, 2333)) # mit i im Dateinamen
{
png(paste("Test_D_",i,".png", sep=""))
plot(i)
dev.off()
}
setwd(oldwd)
unlink("PNG_Test_Bilder", recursive=T) # Seid vorsichtig mit diesem Befehl, der
# kann auch das Skript, mit dem ihr es geschrieben habt, löschen!
# Und alles mögliche andere! Wildcards nur nutzen, wenn ihr genau wisst, was die
# auswählen! Und die gelöschten Dateien landen nicht im Papierkorb!
# 4.12 Beidseitige y-Achsenbeschriftung ----------------------------------------
# Primär und Sekundärachse (Zweite Y-Achse mitbenutzen)
# Lesetipp: http://rwiki.sciviews.org/doku.php?id=tips:graphics-base:2yaxes
# Graphiken mit zwei unkorrelierten Achsen sind problematisch!
# perceptualedge.com/articles/visual_business_intelligence/dual-scaled_axes.pdf
# http://junkcharts.typepad.com/junk_charts/2006/06/illusion_of_suc.html
# http://junkcharts.typepad.com/junk_charts/2006/05/the_crossover_l.html
# Möglichkeit 1, günstig für gleiche y-Bereiche
par(mar=c(2,4,2,4)) # wichtig: rechts mehr Platz machen
a <- c(6,2,9,0,4) ; b <- c(8,7,1,4,6) ; d <-c(43,65,12,23,27)
plot (a, type="l", las=1, ylab="Variable 1")
lines(b, col=2)
axis(4, 2:9, letters[1:8], col.axis=2, las=1)
mtext("Variable 2", 4, 2.5, col=2) # mtext("Text", Achse, Abstand)
points(4, 7, cex=2) # Y-Koordinaten immer von links ("f" geht nicht)
# Bei ungleichen Bereichen muss auf linkes ylim normiert werden
plot (a, type="l", las=1, ylab="Variable 1")
lines(d * 9/65, col=3)
axis(4, pretty(d)*9/65, pretty(d), col.axis=3, las=1)
mtext("Variable 2", 4, 2.5, col=3)
# Möglichkeit 2, da spart man sich die ganze Normierung:
dev.off() # Schließt altes Graphikfenster
werte <- seq(-6,6,.01)
par(mar=c(2,4,3,4))
plot( werte, dnorm(werte), type="l", las=1, ylab="Dichtefunktion",
main="Normalverteilung")
# pnorm und dnorm werden im Kapitel 5 vorgestellt.
par(new=T, xaxt="n")
plot(werte, pnorm(werte), type="l", col=2, ylab="", xlab="", yaxt="n")
axis(4, pretty(0:1), col.axis=2, las=1)
mtext("kumulierte Wahrscheinlichkeit", 4, 2.5, col=2)
legend("topleft", c("dnorm", "pnorm"), col=1:2, lty=1)
abline(h=0.3, lwd=2)
par()$usr # Graphikbereich ist jetzt der vom zweiten Plot!
# 4.13 Achsen unterbrechen, Pakete installieren und laden ----------------------
# Cran Mirror auf Berlin setzen (Siehe Kap 10 für dauerhaftes Setting)
options(repos='http://mirrors.softliste.de/cran')
# Package plotrix installieren (einmalig), kann man anschließend auskommentieren
# R > Pakete > Installiere Paket(e)... > CRAN > plotrix ODER:
install.packages("plotrix")
# offline Tipp: Zips der Pakete in R/Library speichern.
# dann kann man auf einem anderen Rechner per USB-Stick ohne Internet arbeiten.
# Jetzt das Paket laden (Nach jedem Start von R)
library(plotrix) # Mit "Anführungsstrichen" ist auch OK
# Breaks werden zum Plot hinzugefügt, bei Überlappung evt. Daten neuzeichnen
plot(4:15, ylim=c(2,15))
axis.break(2, 2.9, style="slash")
# Plot mit größeren Lücken
gap.plot(c(rnorm(10)+4,rnorm(10)+20), gap=c(8,16),xlab="Index",
ylab="Zufallszahlen",col=c(rep(2,10),rep(3,10)), ytics=c(3,4,5,19,20,21))
legend(12,6,c("High group","Low group"), pch=1, col=3:2)
dev.off()
# 4.14 mehrdimensionale Info in Graphiken --------------------------------------
# Farbe zeigt dritte Dimension:
# Matrix:
volcano[1:10, 1:10]; dim(volcano)
image(volcano)
# Punktdaten:
install.packages("berryFunctions")
library(berryFunctions)
i <- c( 22, 40, 48, 60, 80, 70, 70, 63, 55, 48, 45, 40, 30, 32)
j <- c( 5, 10, 15, 20, 12, 30, 45, 40, 30, 36, 56, 33, 45, 23)
k <- c(175, 168, 163, 132, 120, 117, 110, 130, 131, 160, 105, 174, 190, 183)
colPoints(i,j,k, cex=1.5, pch="+", add=FALSE)
# Dichte der Daten zeigen: sunflowerplot
par(mar=c(3,0,3,0))
sunflowerplot(1,1, xlim=c(1,20))
for (i in 2:20) sunflowerplot(rep(i,i), rep(1,i), xlim=c(1,20), add=T)
title(main="sunflowerplot")
# for - Schleifen wiederholen Operationen. Details in Kapitel 11.
# Es geht hier erstmal um die Darstellung der Sonnenblumen.
# Transparenzplots
x <- rnorm(2000) # 2000 Zufallszahlen aus der Normalverteilung
y <- rnorm(2000) # Siehe Kapitel 5 für Details.
plot(x, y, col = rgb(0, 0, 0, alpha = 0.2), pch = 16)
# kontinuierliche Dichte zeigen: hdr.boxplot.2d
install.packages("hdrcde") # nur einmal
require(hdrcde)
x <- c(rnorm(200,0,1),rnorm(200,4,1))
y <- c(rnorm(200,0,1),rnorm(200,4,1))
par(mfrow=c(1,2))
plot(x,y, pch="+", cex=.5)
hdr.boxplot.2d(x,y)
# Kontinuierliches Histogramm install.packages("UsingR")
library(UsingR)
simple.violinplot(decrease ~ treatment,
data=OrchardSprays, col="bisque", border="black")
# 4.15 Dreidimensionales Plotten -----------------------------------------------
# einige Basic-Befehle
x <- 1:10 ; y <- 1:10 ; z <- matrix(rnorm(100), ncol = 10)
round(z, 2)
filled.contour(x, y, z)
contour(x, y, z)
persp(x,y,z)
# 3d-plot von allen Seiten betrachten
# install.packages("rgl") # nur einmal zum downloaden und installieren.
library(rgl) # in jeder Sitzung.
t <- seq(0,30,0.01) ; x <- 1/pi*t ; y <- sin(t) ; z <- cos(t)
open3d() # Fenster in Sichtweite ziehen
plot3d(x,y,z, type = "l", xlab = "X", ylab = "Y", zlab = "Z")
# mit der Maus klicken und ziehen. Cool, was?
r3dDefaults
x <- seq(-5,5,len=41) ; y <- x
z <- outer(x,y,function(x,y){ifelse(x*x+y*y==0,NA,
sin(3*sqrt(x*x+y*y))/sqrt(x*x+y*y))})
open3d()
plot3d(x,y,z, type="l", xlab="X-Achse", ylab="Y-Achse", zlab="Z-Achse")
# Persp benutzen
?volcano ; v <- volcano
persp(1:nrow(v), 1:ncol(v), v, theta=135, phi=30, shade=0.75, border=NA)
x <- y <- seq(-5,5,len=41)
z <- outer(x,y,function(x,y){ifelse(x+y==0,NA,(x^2-2*x+3*y)/(x+y))})
for (i in 1:360)
persp(x,y,z, theta=i, shade=0.75, border=NA, ticktype="detailed", expand=0.5)
z2 <- outer(x,y,
function(x,y) ifelse( x*x+y*y==0, NA, sin(3*sqrt(x^2+y^2))/sqrt(x^2+y^2) ) )
for (i in 1:360)
persp(x,y,z2, theta=i, phi=i, shade=0.75, border=NA, ticktype="detailed",expand=0.5)
persp(x, y, z, theta=30, phi=30, expand=0.5, col="lightblue")
persp(x, y, z, theta=30, phi=30, expand=0.5, col="lightblue", ltheta=120,
shade=0.75, ticktype="detailed", xlab="X", ylab="Y", zlab="Z")
graphics.off() # schließt alle Graphikfenster
x <- seq(-10, 10, length = 50) ; y <- x
rotsinc <- function(x,y)
{ sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y }
10 * sinc( sqrt(x^2+y^2) ) }
z <- outer(x, y, rotsinc) ; dim(z)
par(bg = "white",mfrow=c(1,2),mar=rep(1.5,4))
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")
title(sub=".")## work around persp+plotmath bug
title( main = expression( z == Sinc(sqrt(x^2 + y^2)) ) )
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed",
xlab = "X", ylab = "Y", zlab = "Z")
dev.off()
# Baumgrößen in 2D-darstellen --------------------------------
graphics.off() # bestehende Devices mit ihren Einstellungen schließen
Orange
str(Orange) # komischer Objekttyp
Orange$Tree <- as.numeric(Orange$Tree) # factor mit anderen Zahlen...
# Siehe auch Beispiele in ?Orange
plot(Orange[Orange$Tree==1, 2:3], main="Orange Trees",
xlim=range(Orange$age), ylim=range(Orange$circumference), type="o")
for(i in 2:5 ) lines(Orange[Orange$Tree==i, 2:3], col=i, type="o")
text(x=c(rep(1600,4), 1250), y=c(130,150,170,195,195), paste("Tree", 1:5),
adj=1, col=1:5)
# Falls gewünscht: für jeden Baum ein eigenes Plot, je mit gleicher Skalierung
par(mfrow=c(2,3), mar=c(3.5,3,2,1), mgp=c(2,1,0))
for(i in 1:5 )
{
plot(Orange[Orange$Tree==i, 2:3], main=paste("Tree Nr.", i), col=i,
xlim=range(Orange$age), ylim=range(Orange$circumference), type="o")
box("figure")
}
# Baumgrößen in 3D ------------------------------------------
open3d()
plot3d(Orange$age, Orange$Tree, Orange$circumference, type="s", col=Orange$Tree)
# Übersichtlicher ist es nicht unbedingt... Aber schön explorativ mit der Maus.
# 4.16 Fehlerbalken (Error bars) -----------------------------------------------
# Funktionen für Konfidenzintervall definieren (Siehe Kap 12)...
cil <- function(dat) return( t.test(dat, conf.level=.95)$conf.int[1] )
ciu <- function(dat) return( t.test(dat, conf.level=.95)$conf.int[2] )
# ... und anwenden:
cil( rnorm(30, mean=5) ) # untere Grenze des 95%igen Konfidenzintervals
# Integrierter Datensatz Iris (siehe Kapitel 7)
head(iris)
str(iris)
# Werte zuordnen ( Vektoren mit Positionen der Fehlerbalken definieren )
mitte <- tapply(iris$Sepal.Length, iris$Species, mean)# Funktion mean anwenden
unten <- tapply(iris$Sepal.Length, iris$Species, cil) # auf die Sepal.Length,
oben <- tapply(iris$Sepal.Length, iris$Species, ciu) # getrennt nach Spezies
# Barplot und x-Positionen der Balken
x.values <- barplot(mitte, ylim=c(0,7), col="cyan")
# Werte anschauen mit 3 Stellen (2 nach Komma)
print( data.frame(x.values, mitte, unten, oben) , digits=3)
# Fehlerbalken einfügen
arrows(x.values, unten, x.values, oben, angle=90, code=3, length=0.3)
# mit length lässt sich die Strichbreite einstellen.
# Oder mit errbar im Paket Hmisc: install.packages("Hmisc")
library(Hmisc)
barplot(mitte, space=0, ylim=c(0,7)) ; axis(1, 0:3)
errbar(1:3-.5, mitte, oben, unten, add=T)
barplot(mitte, space=0, ylim=c(0,7)) ; axis(1, 0:3)
errbar(1:3-.5, rep(NA,3), oben, unten, add=T) ; box()
# Im Paket "gplots" gibt es die Funktion "plotCI" install.packages("gplots")
library(gplots)
?plotCI
# Example with confidence intervals and grid
( hh <- t(VADeaths)[, 5:1] ) ; ci.l <- hh * 0.85 ; ci.u <- hh * 1.15
barplot2(hh, beside=TRUE, legend=colnames(VADeaths), ylim=c(0, 100),
col=c("lightblue", "mistyrose", "lightcyan", "lavender"),
main="Death Rates in Virginia", font.main = 4, plot.grid=TRUE,
sub="Faked 95 percent error bars", col.sub = "gray20",
cex.names=1.5, plot.ci=TRUE, ci.l=ci.l, ci.u=ci.u,)
browseURL("http://cran.r-project.org/doc/contrib/Lemon-kickstart/kr_ptline.R")
# Auch plotrix hat eine Funktion: install.packages("plotrix")
library(plotrix)
plotrix::plotCI(1:3, mitte, ui=oben, li=unten, pch=1, slty="solid", scol=4)
y <- runif(10) ; err.x <- runif(10) ; err.y <- runif(10)
plotCI(1:10,y, err.y, pt.bg=2, pch=21, xlim=c(0,11))
plotCI(1:10,y, err.x, pt.bg=2, pch=21, err="x", add=TRUE)
# 4.17 Fläche unter Graphen füllen ---------------------------------------------
bspx <- c(1,3,4,5,7,9,10,13)
bspy <- c(7,4,9,2,6,5,9,11)
plot(bspx,bspy, ylim=c(0,15))
letzte <- length(bspx)
polygon( c(bspx, bspx[letzte], bspx[1]) , c(bspy, 1, 1), col=2 )
# Einheitsganglinien mit Mittelwert und Konfidenzband
v <- c(0,1,3,5,6,7,8,7,6,5,4,4,3,3,2,2,2,1,1,1,0,0,0,0)/30 ; lv <- length(v)
set.seed(17) # gibt dem Zufallszahlengenerator einen Startpunkt.
( m <- matrix(v+runif(lv*10, 0, 0.13), nrow=lv) )
plot(m[,5], type="n", las=1, ylab="Ordinaten",
xlab="Stunden nach Ereignisbeginn", main="Einheitsganglinie")
for(i in 1:10) lines(m[,i])
cia=0 ; cib=0
for (i in 1:lv) { cia[i]=t.test(unlist(m[i,]), na.rm=T)$conf.int[1]
cib[i]=t.test(unlist(m[i,]), na.rm=T)$conf.int[2] }
cia ; cib
polygon(c(1:lv,lv:1), c(cia, cib[lv:1]), col="orange", border="transparent")
cic=0 ; cid=0
for (i in 1:lv)
{ cic[i]=t.test(unlist(m[i,]), conf.level = 0.7, na.rm=T)$conf.int[1]
cid[i]=t.test(unlist(m[i,]), conf.level = 0.7, na.rm=T)$conf.int[2] }
polygon(c(1:lv,lv:1), c(cic, cid[lv:1]), col=2, border="transparent")
for(i in 1:10) lines(m[,i])
lines(apply(m, 1, mean, na.rm=T), col=3, lwd=3)
par(lend="butt")
legend("topright",
c("Einzelne UHs","95% Konf. Int.", "70% Konf. Int.", "Mittelwert"),
lwd=c(1,7,7,3),col=c(1,"orange",2,3))
mtext("10 Ereignisbasierte Unit Hydrographs sind zuwenig!",1,-2, adj=.4)
# Einfügen von Rechtecken
plot(1:10)
rect(1,5,3,7, col="green")
# polygon, segments, symbols, arrows, curve, abline, points, lines
# http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=154
# SEHR EMPFEHLENSWERTE SEITE!!! (Siehe Source Code für die Befehle)
# 4.18 Funktionen zeichnen ----------------------------------------------------
curve(x^3-3*x, -2, 2)
curve(x^2-2, add = TRUE, col = "violet")
plot(cos, xlim = c(-pi,3*pi), n = 1001, col = "blue")
chippy <- function(x) sin(cos(x)*exp(-x/2))
curve(chippy, -8, 7, n=2001)
curve(chippy, -8, -5)
# Interferenz
x <- seq(-pi,pi,len=1000)
windows(record=T)
par(mar=rep(0,4))
cool <- function (a) plot(a*sin(x), a*cos(x), asp=1, type="l")
# Die for-Schleife zeichnet für jeddes i ein Bild. Mehr in Kap 11.
for (i in 1:500/500) cool(c(i,1)) # Vorischt: höhere Rechen- und plot-zeit
for (i in -600:300/500) cool(c(i,i-.43))
# Jetzt kann man mit Pg Up und Pg Dn die 901 Bilder durchklicken...
# Oder die speichern und als GIF animieren...
# 4.19 zooming in graphics -----------------------------------------------------
install.packages("berryFunctions") # mit meinem Paket ---
require(berryFunctions)
a <- rnorm(90); b <- rexp(90)
plot(a,b, las=1)
pointZoom(a,b,,1,40)
d <- data.frame(a,b)
class(d)
plot(d)
pointZoom(d)
install.packages("playwith") # mit einer tollen interaktiven Geschichte ---
require(playwith)
?playwith
playwith( {plot( 1:10, rnorm(10) )
lines( (1:10) + 0.5, rnorm(10,1,3), col='red' ) })
install.packages("TeachingDemos") # mit einem weiteren Paket ---
require(TeachingDemos)
?zoomplot
# 4.20 und ein würdiger Abschluss des Kapitels ---------------------------------
install.packages("onion")
require(onion)
data(bunny)
p3d(bunny,theta=3,phi=104,box=FALSE)
# 4.21 Übersicht der wichtigsten Graphikbefehle --------------------------------
# Noch to be ergänzt. Für Hinweise wäre ich sehr dankbar!
# high-level plots (fangen neue Graphik an)
plot(x, y)
pairs(data.frame)
coplot(a ~ b | c + d)
qqnorm(x)
qqline(x)
qqplot(x, y)
hist(x)
dotchart(x, ...)
image(x, y, z, ...)
contour(x, y, z, ...)
persp(x, y, z, ...)
barplot(y)
boxplot(y)
curve()
pie()
# low level plots (fügen zur bestehenden Graphik hinzu)
points()
lines() # wie points, mit type="l" als default (spart Tipparbeit)
text()
abline()
polygon()
rect()
legend()
title()
axis()
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 5: statistische Verteilungen #######################################
# Wahrscheinlichkeitsverteilungen (Probability distributions)
# Verteilungen Parameter:
# 5.1 diskret: Binomialverteilung Stichprobenumfang (n), Wahrsch. (p)
# 5.2 diskret: Poissonverteilung Mittelwert (m)
# 5.3 stetig: Normal / Gauss Mittelwert (m), Standardabweichung (sd)
# 5.4 stetig: Exponentialverteilung Rate (r)
# 5.5 stetig: Beta-Verteilung shape1 (a/alpha), shape 2 (b/beta)
# 5.6 weitere Verteilungen
# 5.7 Simulationen
# Vorsilben (Prefixes): Beispiel:
# d probability density function Wahrscheinlichkeitsfunktion dexp(...)
# p cumulative probability function Verteilungsfunktion pbinom(...)
# q quantiles of the distribution qnorm(...)
# r random numbers generated from the distribution rpois(...)
# 5.1 Diskret: Binomialverteilung ----------------------------------------------
# Beispiele: Torwandschießen; Erfolg = Treffer, Misserfolg = kein Treffer
# Ziehen von Kugeln aus Urne; Erfolg = weiße Kugel, Misserfolg = schwarze Kugel
# Relevant ist die Anzahl der Erfolge in n Versuchen
# Wenn es um das Ergebnis eines einzelnen Versuches geht: Bernoulli-Verteilung
# Werte definieren
werte <- seq (from=0, to=30, by=1) # kurz: seq(0,30,1) # kürzer 0:30
werte
# Binomial-Wahscheinlichkeits-Dichte-Funktion:
y.dbinom <- dbinom(x=werte, size=10, prob=0.5) # 10 mal ziehen,
plot(werte, y.dbinom) # Erfolgswahrscheinlichkeit gleich 0,5
?dbinom # (Trefferquote ist 50 %)
# öfter ziehen
plot(werte, dbinom(werte, size=10, prob=0.5), ylim=c(0,0.25), pch=16,
main="Binomialverteilung, Erfolgswahrs. = 0.5", xlab="Anzahl Erfolge")
points(werte, dbinom(werte, 20, 0.5), ylim=c(0,0.25), pch=16, col=2)
points(werte, dbinom(werte, 50, 0.5), ylim=c(0,0.25), pch=16, col=3)
legend("topright", paste("n =", c(10,20,50)), title="Anzahl Versuche", pch=16,
col=1:3, inset=0.02)
# Erfolgsquote verändern
plot(werte, dbinom(werte, 20, 0.85), ylim=c(0,0.25), pch=16,
main="Binomialverteilung, Stichprobenumfang n = 20", xlab="Anzahl Erfolge")
points(werte, dbinom(werte, 20, 0.5), ylim=c(0,0.25), pch=16, col=2)
points(werte, dbinom(werte, 20, 0.2), ylim=c(0,0.25), pch=16, col=3)
legend("topright", paste("prob =",c(.85,.5,.2)), title="Erfolgswahrsch.",
pch=16, col=1:3, inset=0.02)
n <- 8 # Anzahl Münzwürfe: wie wahrscheinlich x mal Kopf?
plot(0:n, dbinom(0:n, n, 0.5), pch=16, xlab="x = Anzahl 'Kopf'", las=1,
main=paste(n, "-facher Münzwurf", sep=""), ylab="Wahrsch. für x-fach Kopf")
abline(h=0.025, col=6)
# Bei 8-maligem Werfen kommt es in 3% der Versuche vor, dass nur 1 mal Kopf
# auftritt. Wenn jemand 8 mal wirft, und nur 1 mal Kopf (bzw. Zahl) auftritt,
# kann ich dennoch nicht annehmen, dass die Münze gezinkt ist (bei 2.5%-Grenze).
# Mord in Palermo / Werwölfe
werte <- 0:10
windows(record=T)
for (i in 0:40) { plot(werte, dbinom(werte, i, 1/8), ylim=c(0,1), xlab="x",
main=paste("Wahrscheinlichkeit, Mörderkarte x mal zu bekommen\nbei",i,"Runden"))
legend("topright", "1 von 8 Karten Mörder") }
# Die Zufallsvariable X sei 'Anzahl Sechser beim viermaligen Würfeln'.
# Die Einzelwahrscheinlichkeiten der fünf möglichen Realisationen sind
# P(X=0) = (5/6)^4 = 0.48225
# P(X=1) = 4(1/6) (5/6)^3 = 0.38580
# P(X=2) = 6(1/6)^2(5/6)^2 = 0.11524
# P(X=3) = 4(1/6)^3(5/6) = 0.01543
# P(X=4) = (1/6)^4 = 0.00077
dbinom(0:4, 4, 1/6)
plot(0:4, dbinom(0:4, 4, 1/6), pch=16)
# Zufallszahlen ziehen
rbinom(n=40, size=10, prob=0.5)
# 5.2 Diskret: Poissonverteilung -----------------------------------------------
# Zähldaten, z.B. Anzahl Telefonanrufe pro Stunde im Callcenter
# oft auch bei geringen Wahrscheinlichkeiten für große Werte (p klein, n groß)
# Beispiel: radioaktiver Zerfall; Aus einer sehr großen Anzahl von Atomen
# zerfällt in einer Zeiteinheit nur ein sehr kleiner Anteil der Atome.
# Dieser Zerfall ist zufällig und unabhängig von den schon zerfallenen Atomen.
# Binomial-Wahscheinlichkeits-Dichte-Funktion
# lambda ist mittleres Auftreten von Ereignissen
plot(0:20, dpois(0:20, lambda = 0.5), pch=16)
points(0:20, dpois(0:20, lambda = 1), pch=16, col=2)
points(0:20, dpois(0:20, lambda = 2), pch=16, col=3)
points(0:20, dpois(0:20, lambda = 5), pch=16, col=4)
points(0:20, dpois(0:20, lambda =10), pch=16, col=5)
# Zufallszahlen ziehen
rpois(n=5, lambda=0.8)
# 5.3 Stetig: Gauss-(Normal)-verteilung ----------------------------------------
# analog zum Vorgehen bei den stetigen Verteilungen
werte <- seq (from = -10, to = 10, by = 0.01)
head(werte)
plot(werte)
plot(werte, dnorm(werte, mean = 5, sd = 1.0))
windows(record=T) # Speichere Bilder, sodass mit Pg Up und Pg Dn durchklickbar
for(i in 1:10) # führe für jedes i folgende Zeilen aus:
{
plot(-20:20, dnorm(-20:20, 0, sd=i ), main=paste("St. Abw. =", i),
xlab="x", ylab="p(x)", type="o")
abline(v=c(0+i, 0-i))
} # Mehr zu For-Schleifen in Kapitel 10
dev.off() # schließe Graphikfenster
# Zufallszahlen als Spalte
zuf.zahl <- rnorm(10, mean = 5.0, sd = 0.5)
data.frame(zuf.zahl)
rnorm(10)
# Plotten von Zufallszahlen der Standardnormalverteilung (mean=0, sd=1)
plot(rnorm(10))
plot(rnorm(1000), pch=16, cex=0.5)
plot(rnorm(500,mean=0,sd=1),rnorm(500,mean=0,sd=1),
pch=16,cex=0.5,abline(v=0, h=0),xlim=c(-4,4), ylim=c(-4,4))
# Variationskoeffizient
x <- -1:130
plot(x, dnorm(x,5,2), type="l", main="Unsinn des Variationskoeffizienten",
ylab="dnorm(x,5,2) | dnorm(x,100,2)")
lines(x, dnorm(x,100,2))
text(10, 0.15, "CV = sd/mean = 0.4", adj=0)
text(105, 0.15, "CV = 0.02",adj=0)
# VK also nur bei Vergleichen zweier Gruppen mit ähnlichem Mittelwert sinnvoll
set.seed(123) # Gibt dem "Zufallszahlen"-Generator den Startpunkt, sodass
# jeder Lauf die gleiche Zahlen ergibt.
d <- rnorm(29); hist(d)
# Homer's Favorite Food - normalverteilt angenähert...
par(mar=c(0,0,0,0), pch="." )
plot(10,10, xlim=c(0,20), ylim=c(0,20))
points( rnorm(100000, 10, 2) + 5*sin(seq(0, 2*pi, len=100000)),
rnorm(100000, 10, 2) + 5*cos(seq(0, 2*pi, len=100000)) )
# Points-Befehl häufiger wiederholen...
# Anwendung kumulierter Wahrscheinlichkeit (pnorm)
# Umfrage zur Körpergröße von Personen
# Mittelwert 170 cm, Standardabweichung 8 cm in Normalverteilung
graphics.off()
plot(150:190, dnorm(150:190, mean=170, sd=8), type="l") ; #abline(h=0:5/100)
par(new=T); plot(150:190,pnorm(150:190,170,8), type="l", axes=F, col=2, ylab="")
axis(4,pretty(0:1), col.axis=2) ; abline(v=c(160,185, 165.8), col=c(4,4,3))
for(h in 1:9*.1) lines( c(qnorm(h,mean=170, sd=8),200), rep(h,2), col=2)
# Wahrscheinlichkeit, dass eine zufällig gewählte Person kleiner als 160 ist
pnorm(160, mean=170, sd=8)
# Wahrscheinlichkeit, dass eine zufällig gewählte Person größer als 185 ist
1 - pnorm(185, mean=170, sd=8)
# Wahrscheinlichkeit, dass eine Person zwischen 160 und 185 groß ist
pnorm(185, mean=170, sd=8) - pnorm(160, mean=170, sd=8)
# Mit 30%iger Wahrscheinlichkeit ist eine Person kleiner als ...
qnorm(.3, mean=170, sd=8)
# 5.4 Stetig: Exponentialverteilung --------------------------------------------
# Beispiel: Zeitabstand zwischen zwei Erdbeben / Naturkatastrophen
# Beispiel: Lebensdauer von elektronischen Bauelementen;
# Die Wahrscheinlichkeit, dass ein älteres Bauteil noch t Tage hält,
# muss genauso groß sein, wie die, dass ein neues Bauteil t Tage hält.
# Die "Ausfallrate" ist konstant (nur für positive Werte).
# analog zum obigen Vorgehen
werte <- seq (from = 0, to = 1, by = 0.001)
plot(werte, dexp(werte, rate = 3), type="l", main="Exponentialverteilung")
lines(werte, dexp(werte, rate = 5), col=2)
text(0.3, 2.5, "rate=5", col=2) ; text(0.1, 1.5, "rate=3")
# Zufallszahlen
rexp(8, rate = 3)
# 5.5 Stetig: Beta-Verteilung --------------------------------------------------
# Zwischen 0 und 1; wird viel in Bayes-Statistik verwendet
werte <- seq (from = 0, to = 1, by = 0.001)
plot(werte, dbeta(werte, shape1=0.5, shape2=5))
# Zufallszahlen
rbeta(8, shape1=3, shape2=1)
# grafisch dargestellt
browseURL("http://www.uni-konstanz.de/FuF/wiwi/heiler/os/vt-index.html")
# Siehe auch library("berryFunctions"); ?betaPlot
# Normalverteilungen beschreiben Werte nicht immer richtig! (meistens nicht):
set.seed(2) ; d <- rbeta(200, 2,4)*350 + 80
hist(d, col=3, breaks=20)
text(350,25, paste("mean =", round(mean(d),2),"\nsd =",round(sd(d),2)), adj=1)
par(new=T)
plot(80:400, dnorm(80:400, mean(d),sd(d)), type="l", lwd=3, axes=F, ann=F)
# wo die Normalverteilung die meisten Werte erwartet, sind weniger!
# 5.6 weitere Verteilungen -----------------------------------------------------
werte <- seq(-3, 3, .01)
plot(werte, dunif(werte, -1, 1), type="l") # density uniform distribution
lines(werte, dnorm(werte), col=2) # "-" normal "-"
plot(werte, punif(werte, -1, 1), type="l") # probability distribution: uniform
lines(werte, pnorm(werte), col=2) # pd: normal (cumulative probability)
# Plotten von Zufallszahlen einer Gleichverteilung: runif(n, von, bis)
# alle Werte (zwischen 0 und 1) sind gleich wahrscheinlich (uniform)
plot(runif(10)) # Random numbers from UNIForm distribution
plot(runif(1000), pch=16, cex=0.5)
# Hier eine Auswahl weiterer Verteilungen, mit Dank an Boris Schröder:
rnorm(n, mean=0, sd=1) # Gaussian (normal)
rexp(n, rate=1) # exponential
rgamma(n, shape, scale=1) # gamma
rpois(n, lambda) # Poisson
rweibull(n, shape, scale=1) # Weibull
rcauchy(n, location=0, scale=1) # Cauchy
rbeta(n, shape1, shape2) # beta
rt(n, df) # 'Student' (t)
rf(n ,df1, df2) # Fisher-Snedecor (F)
rchisq(n, df) # Pearson (Chi^2)
rbinom(n, size, prob) # binomial
rgeom(n, prob) # geometric
rhyper(nn, m, n, k) # hypergeometric
rlogis(n, location=0, scale=1) # logistic
rlnorm(n, meanlog=0, sdlog=1) # lognormal
rnbinom(n, size, prob) # negative binomial
runif(n, min=0, max=1) # uniform
rwilcox(nn, m, n) # Wilcoxon's statistics
rsignrank(nn, n) # Wilcoxon's statistics
browseURL("http://cran.r-project.org/doc/manuals/R-intro.pdf")
# Kapitel 8 Probability distributions (2013-03-01 Seite 41 PDF / 35 print)
browseURL("www.bb-sbl.de/tutorial/verteilungen")
# Eine ganze Reihe (Extremwert-)verteilungen an Daten fitten und analysieren:
# package extremeStat (on github), siehe Kapitel 15.7
# 5.7 Simulationen -------------------------------------------------------------
# Wenn man nicht berechnen kann oder will, kann man oft auch simulieren.
# z.B. die Verteilung der Summenzahl dreier Würfeln:
f <- function() sample(1:6, 1e7, replace=TRUE)
erg <- f() + f()+ f() # Vektorisieren, wann immer es geht
plot(table(erg)/1e7, las=1, main="empirische Wahrscheinlichkeit
für erg Anzahl Augen mit 3 fairen Würfeln")
# z.B. die Wahrscheinlichkeit, drei Mädchen hintereinander zu bekommen:
gender <- t(sapply(1:1e5, function(x) sample(0:1, 3, TRUE) ))
head(gender)
plot(table(rowSums(gender))/1e5, las=1, main="Anzahl Mädchen bei 3 Kindern")
# gilt nur, wenn
# - Geschlecht von dem der Geschwister unabhängig ist
# - die Wahrscheinlichkeit für ein Geschlecht 50% beträgt
# - kein Survivorship Bias besteht (Mädchen öfter abgetrieben als Jungs zB)
# - ...
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 6: Lagemaße, Boxplot, Histogramm, qq-Plot ##########################
# (exploratory data analysis, EDA)
# 6.1 Daten
# 6.2 Histogramme
# 6.3 Lagemaße
# 6.4 Boxplots (und warum man sie nicht so oft verwenden sollte)
# 6.5 QQ-Plots
# 6.1 Daten --------------------------------------------------------------------
# hydraulische Leitfähigkeit (Ksat) in Tiefe A und B
set.seed(28) # set.seed sorgt dafür, dass jeder Lauf die gleichen Zahlen ergibt.
KsA <- round(rexp(40, .02), digits=2)
KsB <- round(rexp(40, .04), digits=2)
KsA # Hier habe ich 40 Zufallszahlen der Exponentialverteilung gewählt.
KsB
Ks <- data.frame(KsA, KsB) ; Ks
# 6.2 Histogramme --------------------------------------------------------------
hist(KsA, col=3, breaks=10)
# zwischen 0 und 20 liegen 15 Werte, usw.
histdata <- hist(KsB, col=2)
histdata
rm(histdata) # wieder aus dem Workspace löschen
par(mfrow=c(1,2))
hist(KsA, cex.axis=0.7, breaks=20, col=3)
hist(KsB, cex.axis=0.7, col=2)
par(mfrow=c(1,1))
# Oder in einem Plot:
hist(KsA, border=3, col=rgb(0,1,0, alpha=.2)) # Alpha für Teiltransparenz
hist(KsB, border=2, col=rgb(1,0,0, alpha=.2), add=T); box()
legend("topright", c("KsA","KsB"), fill=rgb(0:1,1:0,0, alpha=.2), border=3:2 )
# Auswirkung der Anzahl Breaks:
?volcano
hist(volcano, breaks=5, main="5", ylim=c(0,2000), xlim=c(50,200))
hist(volcano, breaks=10, main="10", ylim=c(0,2000), xlim=c(50,200))
hist(volcano, breaks=15, main="15", ylim=c(0,2000), xlim=c(50,200))
hist(volcano, breaks=40, main="40", ylim=c(0,2000), xlim=c(50,200))
hist(volcano, breaks=80, main="80", ylim=c(0,2000), xlim=c(50,200))
hist(volcano, breaks=150, main="150", ylim=c(0,2000), xlim=c(50,200))
# Rand nur außenrum:
h <- hist(volcano, col=3, border=NA)
polygon(x=c(rep(h$breaks, each=2)), y=c(0, rep(h$counts, each=2), 0), lwd=2 )
# 6.3 Lagemaße -----------------------------------------------------------------
# Minimalwert, Maximalwert, Wertespanne
min(KsA)
max(KsA)
range(KsA)
# arithmetisches Mittel vs Median
mean(KsA) # 54 mm/h
median(KsA) # 50% der Werte liegt unter 32 mm/h, 50% der Messwerte sind größer
# Der Mittelwert wird von den wenigen sehr hohen Werten hochgezogen.
# Der Median ist ausreißerunabhängig. (Wobei hohe Ksatwerte keine Messfehler
# sind, also nicht wirklich als Ausreißer zu bezeichnen sind.)
# Beide Maßzahlen für die mittlere Größe einer Variable haben auch
# ein Maß für die Streuung: Standardabweichung (sd) und median Abweichung (mad)
sd(KsA) # 60
# Mit mean und sd kann jeder die Verteilung nachvollziehen:
install.packages("berryFunctions"); library("berryFunctions")
normPlot(mean(KsA), sd(KsA), ylim=c(0, 0.02), cum=F)
hist(KsA, breaks=10, freq=F, add=T, col=2)
# Die Verteilung wird mit den beiden Parametern völlig falsch beschrieben.
# Mittelwert und Stabw sind NUR für normalverteilte Daten sinnvoll (siehe Kap 9)
# MAD: Median der absoluten Abweichungen vom Median (median absolute deviation)
# wird normalerweise normiert; damit das nicht passiert: constant = 1
mad(KsA, constant=1) # durchnittlich liegen die Werte 24 mm vom median entfernt.
# Quantile = Percentile, x% der Werte unter Quantil
quantile(KsA)
# 75% der Werte liegt unter 70 mm/h, dennoch gibt es Werte bis zu 296!
# 25, 50 und 75% werden Quartile (Viertelbereiche) genannt, Q1, Q2 und Q3
# InterQuartilsDistanz: Abstand zwischen erstem (25%) und drittem (75%) Quartil
IQR(KsA) # IQR = Inter Quartil Range, auch Fourth-spread
# 50 % der Werte liegen zwischen Q1 und Q3 mit einer Breite von 54 mm/h
# IQR entspricht der Streungsbreite ohne Extremwerte. Im Boxplot (siehe 6.4)
# wird ein Wert als Ausreißer gesehen, wenn er 1.5*IQR über Q3 / unter Q1 liegt.
# Fünf-ZahlenMaß: Min, 1.Quartil, Median, 3.Quartil, Max
fivenum(KsA)
?boxplot.stats # Section Details erklärt, warum quantile und fivenum Q1 und Q3
# unterschiedlich herausgeben
# Spalten vergleichen: Übersicht einiger Kennzahlen
summary(KsA)
# 9-Zahlenmaß: Median, Quartile (Fourth), Oktile (Eighth), Hexadezile, Range
quantile(KsA, prob=c(0, 1/16, 1/8, 1/4, 1/2, 3/4, 7/8, 15/16, 1))
# Quantile kann man mit mehreren Methoden abschätzen:
quantile(KsA, prob=c(0.90, 0.95, 0.99, 1), type=1)
quantile(KsA, prob=c(0.90, 0.95, 0.99, 1), type=2)
# in Kap 11.5 gibts eine systematische Analyse
# EMPIRISCHE kumulative Verteilungsfunktion
plot(ecdf(KsA), verticals=T, do.points=F)
# 6.4 Boxplots -----------------------------------------------------------------
# Auch Box-Whisker-Plot
# Whiskers: Maximalwert oder, wenn Ausreisser vorhanden, 1.5*IQR (ungefähr 2*sd)
boxplot(KsA) # sind aber nicht wirklich Ausreißer, also lieber:
boxplot(KsA, range=0)
axis(1, c(.5,1,1.5)) # um zu entscheiden, wo Text, Pfeile, etc hin soll.
arrows(.7, 100, .98, 160, col=2, lwd=2)
plot(sort(KsA), pch=16); points(sort(KsB), col=2, pch=22)
Ks <- data.frame(KsA, KsB) ; head(Ks)
boxplot(Ks, range=0)
# Notch (5%-Vertrauensintervall des Median, bei Statistikern umstritten!)
boxplot(KsA, notch=T, xlab="KsA", ylab="Ksat [mm/h]", col=3, range=0)
# Warnmeldung: einige notches liegen außerhalb der hinges ('box'):
# evtl. notch=FALSE setzen -> zeigt, dass Daten breit gestreut sind.
# Mittelwert in Boxplot einzeichnen
boxplot(mean(KsA), add=T, border="red") # add=T: add boxplot to current plot
# Zeichnet einen weiteren Boxplot, bestehend aus nur einem Wert, hinein.
# Lagemaße einzeichnen
abline(h=quantile(KsA, prob=c(1/16,1/8,1/4,1/2,3/4,7/8,15/16)), col ="blue")
# Mehrere Datensätze in ein Fenster boxplotten
boxplot(KsA, KsB, .8*KsA^1.2, notch=T, col=2:4)
boxplot(Ks, notch=T, col=2:3) # automatische Beschriftung an der x-Achse
# Mittelwerte in mehreren Boxen einzeichnen (mit Sternchen)
a <- as.data.frame(matrix(rnorm(8*5), nrow=8, ncol=5))
bp <- boxplot(a)
m <- apply(a, 2, mean) # wende die funktion mean spaltenweise (2) an. --> Kap.11
# Falls nötig: na.rm=T hinter mean: remove missing data (NA)
points(1:5, m, pch=8, col="red", cex=1.5)
# Um auch das mit Boxplots hinzuzufügen, muss m ein data.frame sein.
m ; class(m); is.vector(m) # Sieht aus wie ein dataframe, ist aber ein Vektor!
ma <- as.data.frame( t(m) ) # t = transpose: Zeile <-> Spalte
ma ; class(ma)
boxplot(ma, add=TRUE, border="red")
# Spalten vergleichen: Boxplots
boxplot(KsA, KsB, names=c("Ks A", "Ks B"), notch = T,
ylab = "Ks [mm/h]", col = "grey", main = c("Ks der beiden Tiefen",""))
# Mittelwert addiert
boxplot(data.frame(t(apply(Ks, 2, mean))), add=TRUE, border=2, names=c("",""),
main=c("","arithmetisches Mittel"), col.main="red", cex.main=0.8, notch=T)
# Allerdings verstecken Boxplots die zugrundeliegende Verteilung der Daten.
# Man sieht nur den Median und die Quartile, graphisch sehr ansprechend, aber
# zur Untersuchung von Datenverteilungen ungeeignet. Siehe dieses Beispiel:
set.seed(007) ; dontuseboxplots <- rbeta(100, shape1=2, shape2=4)
X11() # Ordentliches Graphikfenster, falls du in Rstudio arbeitest
op <- par(mfrow=c(3,1))
boxplot(dontuseboxplots, horizontal=T)
title(main="eine leicht rechtsschiefe Normalverteilung, richtig?")
hist(dontuseboxplots, breaks=20, col=8, main="")
title(main="Nichts da: innerhalb des Boxes sind die Werte fast gleichverteilt!")
plot(density(dontuseboxplots), xlim=range(dontuseboxplots), main="")
title(main="Spannend ist auch eine KernDichteSchätzung")
# oder dieses:
par(mfrow=c(4,1), mar=c(2,2,1,1))
set.seed(42) ; dub1 <- c(rnorm(800,5,1), rnorm(200,9,0.5))
dub2 <- rbeta(1000, 6, 50)*39+1.2 # dub: Don't Use Boxplots (unwisely)
boxplot(dub2, col=4, notch=T, ylim=c(1,12), horizontal=T)
boxplot(dub1, col=3, notch=T, ylim=c(1,12), horizontal=T)
title(main="Die sind sehr änhlich, nicht wahr?", line=-1, adj=0.9)
hist(dub2, breaks=30, xlim=c(1,12), col=4, main=""); box()
hist(dub1, breaks=30, xlim=c(1,12), col=3, main=""); box()
title(main="Oh, doch nicht.", line=-1, adj=0.9)
par(op)
# Boxplots machen nur Sinn, wenn man mehr als zwei Datensätze mit gleicher
# zugrundeliegender Verteilung in einem Plot vergleichen möchte.
# 6.5 QQ-Plots -----------------------------------------------------------------
# qqnorm: Vergleich Datensatz zur Normalverteilung
qqnorm(KsA, main = "Ks_A")
qqline(KsA) # Punkte weit weg von der Linie --> Nicht normalverteilt
# Gängig für Ksat-Werte, liegt hier natürlich auch an der Wahl der Daten
plot(sort(KsA)) # das ist letztendlich, was qqnorm macht...
qqnorm(KsB, main = "Ks_B")
qqline(KsB)
qqplot(KsA, KsB) # qqplot: Vergleich zweier (Datensätze) gegeneinander
# Ich verwende eher Histogramme mit möglichst vielen breaks (in den einzelnen
# Klassen sollten trotzdem viele Werte drin sein), um zu sehen, ob Daten
# normalverteilt sind.
# Siehe auch Tests auf Normalverteilungen in Kap 9.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 7: Große Datensätze: iris ##########################################
# 7.1 Datensatz
# 7.2 Subsets (Teildatensätze, indizieren)
# 7.3 Graphiken
# 7.4 EDA - Maßzahlen, Boxplots
# 7.5 Histogramme, QQ-Plots
# 7.1 Datensatz ----------------------------------------------------------------
# Integrierter Datensatz mit Messungen von Kronblättern (Petalen;oben; oft groß)
# und Kelchblättern (Sepalen; unten; meist kleiner) verschiedener Iris Spezies
iris # Sepalen bei Schwertlilien größer als Petalen
class(iris)
str(iris)
summary(iris)
pairs(iris, col=iris$Species)
library(psych)
pairs.panels(iris[,1:4])
# 7.2 Subsets (Teildatensätze, indizieren) -------------------------------------
# Anzeige + Verarbeitung ausgewählter Daten (Siehe Kapitel 3)
iris$Petal.Length # Spalte "Petal.Length" als Vektor
iris["Petal.Length"] # Spalte "Petal.Length" als Spalte
iris[1:50,"Petal.Length"] # Zeilen 1 bis 50 der Spalte "Petal.Length"
iris[iris$Species == "setosa",] # Datensatz für die Species Iris setosa
iris$Petal.Length[iris$Species=="virginica"] # Petal.Length von Iris virginica
# 7.3 Graphiken ----------------------------------------------------------------
# Petal.Length (y) gegen Petal.Width (x)
plot(iris$Petal.Width, iris$Petal.Length, main = "Iris Petal-Dimensionen",
xlab = "Breite der Petalen [cm]", ylab = "Länge der Petalen [cm]",
col = ifelse(iris$Species=="setosa", "red", "blue" ))
# allen drei Arten verschiedene Farben zuordnen
plot(iris$Petal.Width, iris$Petal.Length, col=iris$Species)
plot(sort(iris$Petal.Width), sort(iris$Petal.Length), col=iris$Species)
# allen drei Arten verschiedene ausgewählte Farben zuordnen
plot(iris$Petal.Width, iris$Petal.Length, col=ifelse(iris$Species == "setosa",
"red", ifelse(iris$Species == "versicolor", "blue", "green" )))
# Zuordnungen zu Kürzeln für schnelles Plotten
xl <- "Breite der Petalen [cm]" ; yl <- "Länge der Petalen [cm]"
sl.seto <- iris$Sepal.Length[iris$Species == "setosa"]
sw.seto <- iris$Sepal.Width [iris$Species == "setosa"]
pl.seto <- iris$Petal.Length[iris$Species == "setosa"]
pw.seto <- iris$Petal.Width [iris$Species == "setosa"]
mn.seto <- "Iris Petalen\nI. setosa"
sl.vers <- iris$Sepal.Length[iris$Species == "versicolor"]
sw.vers <- iris$Sepal.Width [iris$Species == "versicolor"]
pl.vers <- iris$Petal.Length[iris$Species == "versicolor"]
pw.vers <- iris$Petal.Width [iris$Species == "versicolor"]
mn.vers <- "Iris Petalen\nI. versicolor"
sl.virg <- iris$Sepal.Length[iris$Species == "virginica"]
sw.virg <- iris$Sepal.Width [iris$Species == "virginica"]
pl.virg <- iris$Petal.Length[iris$Species == "virginica"]
pw.virg <- iris$Petal.Width [iris$Species == "virginica"]
mn.virg <- "Iris Petalen\nI. virginica"
# Anmerkung: mit der Funktion by() ließen sich einige Zuordnungen umgehen.
# mehrere Graphiken auf einem Blatt: relative Größen
layout(matrix(1:3, 1, 3), widths=c(5, 5, 5))
plot(pw.seto, pl.seto , main=mn.seto, xlab=xl, ylab=yl, ylim=c(1,7))
plot(pw.vers, pl.vers , main=mn.vers, xlab=xl, ylab=yl, ylim=c(1,7))
plot(pw.virg, pl.virg , main=mn.virg, xlab=xl, ylab=yl, ylim=c(1,7))
# mehrere Graphiken auf einem Blatt: absolute Größen (in cm)
# zB. für quadratische plots. Grafikfenster ausreichend groß ziehen!
layout(matrix(1:3, 1 ,3), widths=lcm(c(5, 5, 5)), heights= lcm(5))
plot(pw.seto, pl.seto , main=mn.seto, xlab=xl, ylab=yl, ylim=c(1,7))
plot(pw.vers, pl.vers , main=mn.vers, xlab=xl, ylab=yl, ylim=c(1,7))
plot(pw.virg, pl.virg , main=mn.virg, xlab=xl, ylab=yl, ylim=c(1,7))
# Zurücksetzen der Graphik
dev.off()
# 7.4 EDA - Maßzahlen, Boxplots ------------------------------------------------
# 5-Zahlenmaß + mean für jeden Parameter, nach Iris-Arten getrennt
summary(iris[iris$Species == "setosa",])
summary(iris[iris$Species == "versicolor",])
summary(iris[iris$Species == "virginica",])
# Boxplot-Vergleiche (zwischen den 3 Iris-Arten) pro Parameter
par(mfrow = c(2,2))
boxplot(sl.seto, sl.vers, sl.virg, notch = T, ylab = "Länge der Sepalen [cm]",
names = c("I.setosa", "I.versicolor", "I.virginica"), ylim=c(0,8))
means <- data.frame(t(c(mean(sl.seto),mean(sl.vers),mean(sl.virg))))
boxplot(means, add=TRUE, border=2, names=c("", "", ""))
boxplot(sw.seto, sw.vers, sw.virg, notch = T, ylab = "Breite der Sepalen [cm]",
names = c("I.setosa", "I.versicolor", "I.virginica"), ylim=c(0,8))
means <- data.frame(t(c(mean(sw.seto),mean(sw.vers),mean(sw.virg))))
boxplot(means, add=TRUE, border=2, names=c("", "", ""))
boxplot(pl.seto, pl.vers, pl.virg, notch = T, ylab = "Länge der Petalen [cm]",
names = c("I.setosa", "I.versicolor", "I.virginica"), ylim=c(0,8))
means <- data.frame(t(c(mean(pl.seto),mean(pl.vers),mean(pl.virg))))
boxplot(means, add=TRUE, border=2, names=c("", "", ""))
boxplot(pw.seto, pw.vers, pw.virg, notch = T, ylab = "Breite der Petalen [cm]",
names = c("I.setosa", "I.versicolor", "I.virginica"), ylim=c(0,8))
means <- data.frame(t(c(mean(pw.seto),mean(pw.vers),mean(pw.virg))))
boxplot(means, add=TRUE, border=2, names=c("", "", ""))
par(mfrow = c(1,3))
boxplot(iris[iris$Species=="setosa",1:4], ylim=c(0,8), main="Iris setosa")
boxplot(iris[iris$Species=="versicolor",1:4],ylim=c(0,8),main="Iris versicolor")
boxplot(iris[iris$Species=="virginica",1:4], ylim=c(0,8), main="Iris virginica")
par(mfrow = c(1,1))
# Zusammenhänge zwischen Elemente in großen Datensätzen
pairs(iris, pch=16, cex=0.8, col=ifelse(iris$Species=="setosa",
"red", ifelse(iris$Species=="versicolor", "blue", "green")))
# 7.5 Histogramme, QQ-Plots ----------------------------------------------------
# Histogramme des Parameters Sepal.Length
par(mfrow=c(1,3))
sl <- "Länge der Sepalen [cm]"
hist(sl.seto, xlab=sl, ylab = "Häufigkeit", main = "I. setosa",
xlim=c(4,8), ylim=c(0,18), breaks=seq(4, 8, 0.25))
hist(sl.vers, xlab=sl, ylab = "Häufigkeit", main = "I. versicolor",
xlim=c(4,8), ylim=c(0,18), breaks=seq(4, 8, 0.25))
hist(sl.virg, xlab=sl, ylab = "Häufigkeit", main = "I. virginica",
xlim=c(4,8), ylim=c(0,18), breaks=seq(4, 8, 0.25))
# QQ-Plots des Parameters Sepal.Length mit eingezeichneten QQ-Linien
qqnorm(sl.seto, main = "I.setosa")
qqline(sl.seto)
qqnorm(sl.vers, main = "I.versicolor")
qqline(sl.vers)
qqnorm(sl.virg, main = "I.virginica")
qqline(sl.virg)
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 8: Bayes-Statistik und -Hypothesentests ############################
# 8.1 Paket: Bolstad
# 8.2 Urnenmodell
# 8.3
# 8.4 Wasserqualität
# 8.5 Hypothesentests
# 8.1 Paket: Bolstad -----------------------------------------------------------
# R > Pakete > Installiere Paket(e)... > CRAN > Bolstad
# oder kurz:
install.packages("Bolstad") # Paket Downloaden, entpacken, installieren
# Danach auskommentieren, muss nur einmal ausgeführt werden!
library(Bolstad) # lädt das Paket, muss in jeder Sitzung geschehen
# Für Binomial sampling with a discrete prior steht folgende Funktion bereit:
?binodp(x,n, pi=NULL, pi.prior=NULL, n.pi=10)
# 8.2 Urnenmodell --------------------------------------------------------------
# Binomialverteilung, diskreter Prior
# n = 5 Kugeln rot/grün; 6 gleichwahrsch. Hypothesen über Anzahl roter Kugeln
# 1 rote gezogen (Erfolg) -> Funktion binodp (Seite 5 Bolstad Manual)
binodp(1, 1, pi = rep(1/6, 6), pi.prior=NULL, n.pi=6 )
# 8.3 Mittelwertschätzung ------------------------------------------------------
# Inferenz: Schluß von einer Stichprobe auf die Grundgesamtheit
# Schätzung des Mittelwertes unter Berücksichtigung von Vorwissen / Schätzung
# Berechnung der Likelihood und des Posteriors eines jeden Wertes
# 5 Hypothesen des Mittelwertes des Gewichts von Studentinnen in Potsdam
mu <- c(55, 60, 64, 66, 73) # mu-Hypothesen
prior <- rep(0.2, 5) # jede ist gleichwahrscheinlich richtig
stichprobe <- c(72, 65, 60, 71, 70, 60, 59)
Bolst <- normdp(x=stichprobe, sigma.x=sd(stichprobe), mu=mu, mu.prior=prior)
# Plot nicht vergessen anzuschauen 😉
Bolst
round(Bolst$posterior*100,2) # 0.00 2.71 45.62 51.58 0.09
# Gegeben diese Stichprobe, dass sind die Wahrscheinlichkeiten, dass die
# jeweilige Hypothese richtig ist.
# In diesem Fall würde man in der Realität kontinuierlich rechnen lassen:
normgcp(x=stichprobe, sigma.x=sd(stichprobe), density="uniform",
params=c(45,80), n.mu=80)
# Der Posterior wird nur für n.mu diskrete Werte ausgegeben.
# 8.4 Wasserqualität -----------------------------------------------------------
# Water Quality: n=116 water samples; y=17 samples contaminated.
# Binomialverteilung (n=116, pi) (siehe Kap. 5)
# pi = true probability that a water sample is contaminated
# posterior distribution of pi given y with a beta(l,4) prior for pi
binobp(17, 116, a=1, b=4)
binobp # zeigt die zugrundeliegende Funktion (Siehe Kap 12)
# Frequentist estimator (Erfolgswahrscheinlichkeit) = Zahl Erfolge / Samples
( pi.f <- 17/116 ) # = 0.09
# normal approximation to the posterior distribution g(pi|y).
werte <- seq(0,1, len=10000)
plot(werte ,dnorm(werte, mean=0.1487603, sd=0.0322173), type="l", ylim=c(0,14),
main="n=116; y=17; beta(1,4)-prior", ylab="Density", xlab="pi")
# 95% credible interval for pi with the normal approximation:
abline(v=qnorm(c(0.025, 0.975), mean=0.1487603, sd=0.0322173), col=2)
text(0.5, 10, "95% Interval with normal approx", col=2)
# and with the binobp-approach:
abline(v=c(0.0913896, 0.2171069), col=3)
text(0.45, 8, "95% Interval with binobp", col=3)
# Ergebnistabelle abspeichern (Datei wird im Working Directory angelegt),
getwd() # zeigt dieses an
sink(file="bayes-output.txt", append=T , split=T)
# append=T: outputs werden nicht überschrieben, sondern angehängt
# split=T : output wird auch in der console noch angezeigt
# Benutzung eines anderen Priors: Betrachte Mittelwert und Varianz.
binobp(17, 116, a = 1, b = 4)
binobp(17, 116, a = 0.5, b = 1)
binobp(17, 116, a = 1, b = 0.5)
binobp(17, 116, a = 1, b = 1)
binobp(17, 116, a = 2, b = 1)
binobp(17, 116, a = 3, b = 2)
# Beispiel für weitere Daten
binobp(10, 174, a=1, b=4)
# Schreiben in Datei beenden
sink()
# Posterior distribution g(pi|y) mit beta(1,lO) als prior für pi
# Ziehungen = 116, Erfolge = 17, Misserfolge = 116-17 = 99
# a' = a.neu = a.alt + Zahl.Erfolge b'= b.neu = b.alt + Zahl.Misserfolge
# hier: a.neu = 1 + 17 b.neu = 10 + 116 - 17
binobp(17, 116, a=1, b=10)
binobp(17, 116, a=18, b=109) # posterior-beta-Funktion beta.neu = beta(18, 109)
# 95% credible interval for pi: R-Ausgabe Quantile 0.025, 0.975, oder
qbeta(c(0.025, 0.975), 18, 109)
# Vergleich mit Normalverteilung mit gleichem Mittelwert / Standardabweichung
qnorm(c(0.025, 0.975), mean= 0.1440329, sd= 0.0224784)
lines(0:100/100, dnorm(0:100/100, mean=0.1440329, sd=0.0224784), col=3, lty=2)
# 8.5 Hypothesentests ----------------------------------------------------------
# H0: pi = .l0 H1: pi != .10 Level of significance: 5%
# Nullhypothesen-Wert pi= .10 liegt im Wahrscheinlichkeitsintervall
# Nullhypothese wird nicht verworfen
# H0: pi>= .15 H1: pi < .15 Level of significance: 5%
# kumulierte Wahrscheinlichkeitsverteilung (pbeta) der posterior-beta-Verteilung
# Gesucht ist der Wert, der in der Wahrscheinlichkeitsdichte-Verteilung der
# Fläche unter der Kurve rechts des Wertes 0.15 entspricht. In der kumulierten
# Wahrscheinlichkeitsverteilung bedeutet dies 1-(Wert bei 0.15)
1-pbeta(.15, shape1= 18, shape2=109)
# Dieser Wert 0.37 ist größer als das Signifikanzlevel 0.05
# Die Null-Hypothese kann nicht verworfen werden
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 9: klassische Statistik ############################################
# 9.1 Hypothesentests Allgemein
# 9.2.1 Tests auf Normalverteilung
# 9.2.2 Transformation hin zur Normalverteilung
# 9.3 Tests auf Varianzenhomogenität
# 9.4 Unterschied zweier Datensätze
# 9.5 Anova
# 9.1 Hypothesentests Allgemein ------------------------------------------------
# P-Value: Wahrscheinlichkeit, ein Ergebnis (Datenverteilung zB) wie im Test
# zu bekommen, wenn die Nullhypothese zutrifft.
# Wenn p < alpha: H0 ablehnen alpha= 1-Sign.Niveau, oft 0.05 bei = 95%
# alpha = Wahrsch., Fehler 1. Art zu machen (H0 ablehnen, obwohl zutreffend)
browseURL("http://de.wikipedia.org/wiki/Statistische_Signifikanz")
# Hypothesentest sind schwarz-weiss-Malerei, wo die Welt doch farblich (und
# ultraviolett) ist, machen also nur in 0-1-Entscheidungen Sinn.
# Sind immer noch kritisch zu betrachten, siehe zB:
browseURL("http://www.eecs.qmul.ac.uk/~norman/blog_articles/p_values.pdf")
A <- rnorm(30, mean=154) ; B <- rnorm(30, mean=392)
t.test(A, B)
wilcox.test(A, B)
cor.test(A, B, method="p")
# Überprüfung der Signifikanz der Übereinstimmung zweier Verteilungen
# pairwise.wilcox.test(A, B, p.adjust.method="bonf", paired=FALSE)
fakedaten <- rnorm(1000, mean=97, sd=8.9)
t.test(fakedaten, mu=100, alternative="less") # p-value < 2.2e-16 < alpha
# also H0 ablehnen, und davon ausgehen, dass der wahre Mittelwert aller
# potentiellen Werte kleiner als 100 ist.
hist(fakedaten, col=8, freq=F)
lines(density(fakedaten)) # kerndichteschätzung
lines(70:130, dnorm(70:130,mean=mean(fakedaten),sd=sd(fakedaten)), col=2, xpd=T)
lines(70:130, dnorm(70:130,mean=97,sd=8.9), col=3, xpd=T)
# 9.2.1 Tests auf Normalverteilung ---------------------------------------------
browseURL("www.bb-sbl.de/tutorial/verteilungen/ueberpruefungnormalverteilung")
?binom.test
?chisq.test
shapiro.test(fakedaten)
ks.test(fakedaten, "pnorm", mean(fakedaten), sd(fakedaten))
# Wenn p > 0.05: Nullhypothese annehmen, daten sind normalverteilt.
if(!"nortest" %in% installed.packages()[,1]) install.packages("nortest")
library(nortest) # mit Lilliefors-Korrektur!
help(package="nortest")
?lillie.test
lillie.test(fakedaten)
ad.test(fakedaten) # Anderson-Darling test for normality
cvm.test(fakedaten) # Cramer-von Mises test for normality
lillie.test(fakedaten) # Lilliefors (Kolmogorov-Smirnov) test for normality
pearson.test(fakedaten) # Pearson chi-square test for normality
sf.test(fakedaten) # Shapiro-Francia test for normality
# verschiedene Möglichkeiten zur Transformation
# 1/x, 1/sqrt(x), log(x), x^(1/4), x^(1/3), sqrt(x), x^2, x^3, boxcox (unten)
# 9.2.2 Normalisierung mittels BoxCox-transformation ---------------------------
browseURL("http://de.wikipedia.org/wiki/Box-Cox-Transformation")
browseURL(
"http://en.wikipedia.org/wiki/Power_transform#Box.E2.80.93Cox_transformation")
# Daten generieren
set.seed(42) ; Ksat <- rbeta(150, 0.3417, 1.8320)*220
# Verteilung der Daten ist nicht gaussisch:
hist(Ksat, breaks=50, col=8, freq=F)
x <- seq(-5, 200, len=200)
lines(x, dnorm(x, mean(Ksat), sd(Ksat)), lwd=2, col=5) # doesn't fit!
if(!require(car)) {install.packages("car") ; library(car) }
lambda <- powerTransform(Ksat)$lambda ; lambda # 0.2254471
Ksatbc <- bcPower(Ksat, lambda) #bc: boxcox
hist(Ksatbc, col=8, breaks=20, freq=F)
lines(x, dnorm(x, mean(Ksatbc), sd(Ksatbc)), lwd=2, col=5) # fits much better
# Alternative, hier aber gleiches Histogramm als Ergebnis
library(geoR)
lambdageoR <- boxcoxfit(Ksat)$lambda ; lambdageoR # 0.2254406
hist(bcPower(Ksat, lambdageoR), col=8, breaks=20, freq=F)
# für Lineare Modelle siehe
library(MASS)
?boxcox() # Box-Cox Transformations for Linear Models
# 9.3 Tests auf Varianzenhomogenität -------------------------------------------
fligner.test( A ~ B )
# p < 0.05 impliziert: Varianzen sind unterschiedlich
var.test() # anfällig gegenüber Ausreißern
# install.packages("car")
library(car)
leveneTest()
summary(aov( A ~ B))
TukeyHSD(aov(A ~ B))
kruskal.test(A, B)
# 9.4 Unterschied zweier Datensätze --------------------------------------------
?t.test
# Wann welcher Signifikanztest sinnvoll ist, ist hier dargestellt:
browseURL("http://de.wikibooks.org/wiki/GNU_R:_Signifikanztests")
# There are two groups of values. The confidence interval overlaps, which means
# that the two groups do not significantly differ.
# However, the t-test returns a p-value smaller than the divide of 0.05,
# indicating a significant difference of the means.
# Why are the resulting interpretations contradictory?
# If two confidence intervals overlap, the difference between the two means
# still may be significantly different.
browseURL("http://www.jstor.org/stable/view/2983411")
browseURL("http://www.cmaj.ca/content/166/1/65.long")
# Group 1 # Echte Daten aus einer Psychologie-Bachelorarbeit (2013)
G1 <- c(92,98,102,108,99,99,88,97,114,102,114,104,101,88,102,103,102,86,85,90,
92,93,102,100,100,96,99,97,97,103,101,105,92,93,92,91,102,105,87,92,100,104,99,
93,102,94,88,94,91,100,93,83,85,90,88,93,99,100,104,102,87,92,104,106,88,98,92,
90,92,94,93,101,96,92,94,92,90,99,92,97,93,97,98,103,99,95,94,102,98,99,99,103,
107,111,106,108,104,108,97,94,92,90)
# Group2
G2 <- c(115,107,111,108,102,83,94,90,102,109,106,90,113,85,109,103,102,109,94,
93,82,96,95,100,107,93,102,108,100,118,88,87,86,87,86,85,93,95,81,99,106,97,90,
93,84,83,98,119,98,100,87,80,94,90,98,99,114,103,117,123,99,106,107,126,94,91,
86,90,95,91,90,105,93,89,97,107,102,96,102,104,88,94,109,103,102,100,98,102,95,
108,92,94,95,99,102,108,97,96,107,123,83,81,89,88,114,82,96,108,109,104,93,121,
116,117,113,110,87,95,112,97,101,102,105,102,102,97,102,105,99,92,90,86,91,116,
92,93,109,107,114,108,111,105,116,114,111,96,106,104,106,101,94,100,102,97,81,
87,85,99,89,82,90,93,109,111,100,96,101,102,86,93)
# the graphic with confidence intervals
CI_1 <- t.test(G1)$conf.int
CI_2 <- t.test(G2)$conf.int
plot(1:2, 1:2, ylim=c(95,101), type="n", xaxt="n", ann=F, las=1)
arrows(x0=1.4, y0=CI_1[1], y1=CI_1[2], angle=90, code=3, length=0.2)
points(1.4, mean(G1), pch=16)
arrows(x0=1.6, y0=CI_2[1], y1=CI_2[2], angle=90, code=3, length=0.2)
points(1.6, mean(G2), pch=16)
# Confidence Intervals around the mean overlap. WRONG interpretation: so the
# two groups are not significantly different from each other
t.test(G1, G2) # p-value = 0.0349 < 0.05 -> reject H0, accept Alternative
# true difference in means is not equal to 0 -> signif. diff. of mean of groups
# at conf.lev=0.95, you will WRONGLY reject H0 in 5% of (endless many) samples!
# Control interpretation of p-value
t.test(rnorm(100, mean=1), rnorm(100, mean=15)) # p-value < 2.2e-16
# -> difference is significant
# control settings of t.test-argument var.equal
var.test(G1, G2) # p-value = 4.167e-06 -> true ratio of var.s is not equal to 1
var(G1) ; var(G2) # 42 and 100 -> var.equal = FALSE is correct
# Konfidenzintervälle der klasisschen (frequentistischen) Statistik sind böse:
browseURL("http://www.entsophy.net/blog/?p=53") # und viele mehr!
# 9.5 Anova --------------------------------------------------------------------
# ANOVA: ANalysis Of VAriances (t.test for >2 groups)
# Compare the ratio of between group variance to within group variance.
# If the variance caused by the interaction between the samples is much larger
# when compared to the variance that appears within each group,
# then it is because the means aren't the same.
browseURL("http://sociology.about.com/od/Statistics/a/Analysis-of-variance.htm")
browseURL("http://people.richland.edu/james/lecture/m170/ch13-1wy.html")
browseURL("http://en.wikipedia.org/wiki/F_test#One-way_ANOVA_example")
# AOV Beispiel nach:
browseURL("http://users.minet.uni-jena.de/~jschum/biostat/ANOVA.pdf")
peas <- read.table(header=T, text=
"Wdh Kontr. Gluk.2 Frukt.2 Gluk.1_Sacch.2 Frukt.1
1 71 57 58 58 62
2 68 58 61 59 66
3 70 60 56 58 65
4 74 59 58 61 63
5 68 62 57 57 64
6 71 60 56 56 62
7 70 60 61 58 65
8 67 57 60 57 65
9 73 59 57 57 62
10 69 61 58 59 67")
peas
peas <- stack(peas[,-1])
names(peas) <- c("length", "group")
head(peas, 15)
summary(peas)
tapply(peas$length, peas$group, mean)
par(mar=c(3,6,2,1), las=1)
boxplot(length~group, data=peas, col=5, horizontal=T)
AOV <- aov(length ~ group, data=peas)
AOV
summary(AOV)
par(mfrow = c(2, 2))
plot(AOV)
# ANOVA Beispiel nach:
browseURL("www.stat.wisc.edu/~yandell/st571/R/anova.pdf")
# Daten erstellen:
y1 <- c(18.2, 20.1, 17.6, 16.8, 18.8, 19.7, 19.1)
y2 <- c(17.4, 18.7, 19.1, 16.4, 15.9, 18.4, 17.7)
y3 <- c(15.2, 18.8, 17.7, 16.5, 15.9, 17.1, 16.7)
y = c(y1, y2, y3)
group <- rep(1:3, each=7)
# exploratory data analysis (EDA):
plot(group, y, pch=16, col=rgb(0,0,0, 0.2)) # keine überlappenden Punkte
boxplot(y~group)
info <- function(x) c(sum=sum(x), mean=mean(x), var=var(x), n=length(x))
info(y)
tapply(X=y, INDEX=group, FUN=info) # Wende die Funktion "info" für jede "group"
# getrennt auf "y" an (anwenden = apply). Siehe auch Kapitel 11.3
# Linear model and ANOVA
mod <- lm(y~group) # Mehr dazu im Kapitel 14 (Regression)
mod
plot(y~group); abline(mod)
summary(mod)
str(mod)
mod_aov <- anova(mod)
mod_aov
methods("print")
?print.anova
summary(mod_aov)
str(mod_aov)
# tabled F values:
Df <- mod_aov$Df # Degrees of freedom of treatment and error
Df
qf(p=c(0.05, 0.01), df1=Df[1], df2=Df[2], lower.tail=FALSE)
# confidence interval on the pooled variance:
SSTrt <- mod_aov["Residuals", "Sum Sq"] # SSTrt: residual sum of squares
SSTrt
SSTrt/qchisq(p=c(0.025, 0.975), df=Df[2], lower.tail=FALSE)
# Versuch, das per Hand zu machen. Klappt noch nicht:
# One-way between groups ANOVA variation != variance
m1 <- mean(G1); m1 # group mean
m2 <- mean(G2); m2
n1 <- length(G1); n1
n2 <- length(G2); n2
mg <- mean(c(G1,G2)); mg # overall mean (Grand mean)
ssw <- var(G1)*(n1-1) + var(G2)*(n2-1); ssw # within group variation SS(W)
sum((G1-m1)^2, (G2-m2)^2) # the same
ssb <- n1*(m1-mg)^2 + n2*(m2-mg)^2; ssb # between group variation SS(B)
(ssb/1) / (ssw/(n1+n2-2)) # stimmt so noch nicht!
var.test(G1,G2)$stat
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 10: if, for (Programmieren) ########################################
# 10.1 Bedingungen (conditions): if, ifelse
# Schleifen (loops): for, repeat, while
# 10.2 Allgemeines zu for-Schleifen
# 10.3 repeat- und while-Schleifen
# 10.4 Anwendungsbeispiel
# 10.5 Geschwindigkeit for-Schleifen
# 10.1 Bedingungen (conditions) ------------------------------------------------
# Syntax: if(Bedingung) {Ausdruck1} else {Ausdruck2} bei 1 Wahrheitswert
# Syntax: ifelse(Bedingung, Ausdruck1, Ausdruck2) für Vektor mit mehreren T/F
# Wenn Bedingung wahr, dann Auswerten von Ausdruck 1,
# wenn Bedingung falsch, dann Auswerten von Ausdruck 2.
7-3 > 2 # TRUE
class(7-3 > 2 ) # logical, Wahrheitswert
if(7-3 > 2) 18 # Bedingung ist TRUE, also gebe 18 zurück
if(7-3 > 5) 18 # Bedingung ist FALSE, also passiert nichts
if(7-3 > 5) 18 else 17 # Bed. F, also Rückgabe von 17
par(mfrow=c(2,1), cex=1)
if(TRUE) plot(1:10) else plot(10:1) # Bedingung True, also aufsteigende Punkte
# Mehrere Befehle müssen mit {} zusammengehalten werden:
if(FALSE)
{
plot(1:10)
box("figure", col=4, lwd=3)
} else
{
plot(10:1)
box("figure", col=2, lwd=3, lty=2)
}
# Code einrücken macht ihn menschenlesbar!
# Führe den Block auch mal mit einer erfüllten Bedingung aus
# Tinn R: R Send: contiguous (echo=T) # Siehe Kap. 10.1 falls .trPaths fehlt.
# Unterschied zwischen "if(..) .. else .." und "ifelse(...)":
# if wenn die Bedingung nur einen Wert hat
# ifelse wenn ein ganzer Vektor mit Bedingungen vorliegt.
bart <- c(13, 14, 15, 16, 17)
bart>14
bart
ifelse(bart>14, bart+1, bart-10)
if(bart>14) bart+1 else bart-10 # Nimmt nur 1. Element, und warnt dabei.
simpson <- -3:8
sqrt(simpson) # Warnt vor Wurzel aus negativen Zahlen
sqrt(ifelse(simpson >= 0, simpson, NA)) # läuft
sqrt(if(simpson >= 0) {simpson} else NA) # geht nicht, denn simpson ist Vektor.
?"if"
# Zeilenumbrüche zwischen '}' und 'else' sind zu vermeiden, wenn via Editor
# die Befehle Zeilenweise an R geschickt werden!
ergebnis <- 1/4 # Ergebnis vorheriger Berechnungen
if( ergebnis > 0.2 ) message("Ergebnis ist > 0.2 und das ist relevant.\n")
# Sowas kann man gut in Funktionen einbauen!
# Eine ganze Reihe von Bedingungen ist daher in Kapitel 12 verwendet.
# Zur Abkürzung T und F siehe Kapitel 19.5 Logische Operatoren
# Wenn Rechenzeit für dich wichtig wird:
# "ifelse() is for vectorized conditions. if(){} else{} is more efficient
# for length 1 conditions than ifelse()" Uwe Ligges, R-help 2004
# Schleifen (loops) ------------------------------------------------------------
# 10.2 Allgemeines zu for-Schleifen --------------------------------------------
# mehrfaches Ausführen eines Codeblocks, zB Erstellen mehrerer Bilder
# Syntax: for ( EinLaufbuchstabe in EineSequenz) { machwas }
# Zwischen den {} können mehrere Befehle stehen, die i-fach / für i verschiedene
# Werte durchgeführt werden, auch über mehrere Zeilen getrennt.
# Bei einem einzelnen Befehl können die {} weggelassen werden.
help("for")
for(i in c(2,5,9) ) print(1:i) # man verwendet oft i für index
# im ersten Durchlauf wird an den Befehl i mit dem Wert 2 übergeben.
# Der Befehl wird ausgeführt, dann fängt die Schleife von vorne an und führt
# den Code mit dem nächsten Wert aus.
# mit for einen Vektor füllen
v <- vector("numeric", 20)
v
for(i in 3:17) v[i]= (i-2)^2 # für jedes i wird der Code einmal ausgeführt
v
# Ist in R sehr langsam - vektorisieren ist 100 fach schneller. -> 10.2.4
# Erstellen mehrerer Graphiken. Besonders sinnvoll in Kombination mit anderen
# Devices, zB png, bmp oder pdf (Siehe Kap 4.10)
windows(record=T) # Aufzeichnung der Einzelbilder starten
plot(1, ylim=c(0,7))
plot(3, ylim=c(0,7))
for (i in 1:7) plot(i, ylim=c(0,7)) # mit Bild Up, Bild Dn durchklicken
# Addieren mehrerer Linien zum Plot
z3 <- seq(0,25,0.1)
plot(z3,z3, type="n")
lines(z3, 12.5*5/factorial(4)*(z3/8*5)^(4)*exp(-z3/8*5), col=2 )
lines(z3, 12.5*6/factorial(5)*(z3/8*6)^(5)*exp(-z3/8*6), col=2 )
for (n in 7:25)
lines(z3, 12.5*n/factorial(n-1)*(z3/8*n)^(n-1)*exp(-z3/8*n), col=2 )
# Laufende Dateinamen oder Objektnamen erstellen:
getwd() # mit setwd auf einen Ordner mit Schreibberechtigung setzen
# Ja, beschränkte Rechte auf C:\ sind für viele (unerfahrene) Nutzer sinnvoll!
meinevar <- c(12,14,15,28,76)
for(k in 1:4) { out <- rexp(k)
write.table(out, paste("MeinName", meinevar[k], ".txt", sep="") , quote=F)}
cat("Hier sind deine Textdateien:\n", getwd(), "\n") # \n newline, Kap 17.1
# Jetzt Dateien wieder löschen
unlink(paste("MeinName", meinevar[1:4], ".txt", sep="") ) # riskanter Befehl!
# Das Skript kann sich zB selbst löschen, wenn man den ganzen Ordner löscht!
# 10.3 repeat- und while-Schleifen ---------------------------------------------
# Wiederholung von Befehlsaufrufen, Iterationen, Simulationen
# Syntax: repeat {Ausdruck} # Wiederholung des Ausdrucks
# next # Sprung zum nächsten Iterationsschritt
# break # Verlassen der Schleife
# Syntax: while(Bedingung){Ausdruck} # Wiederholung, solange Bed. erfüllt
i <- 0 # Lege Objekt i mit dem Wert 0 an
repeat # wiederhole nachfolgenden Codeblock immer wieder
{ # Codeblock beginnt
i <- i + 1 # überschreibe Objekt i mit neuer Zahl
if (i < 3) print(i) # wenn i kleiner ist als 3, dann Ausgabe in der Konsole
if (i == 3) break # wenn i gleich 3 ist, wird die Schleife abgebrochen.
} # Codeblock endet
i # beim letzten Durchlauf wurde 1 zu 2 hinzugefügt (=3), und das als i angelegt
i <- 0
while(i < 3) i <- i + 1 # Während (while) 1<3 ist, wird Code ausgeführt
i
# Werden beides in der Umwelt-Datenanalyse sehr selten angewendet, repeat habe
# ich sogar noch nie in der Realität gebraucht.
# 10.4 Anwendungsbeispiel ------------------------------------------------------
# Punkte zufällig festlegen, aber mit Mindestdistanz zueinander.
set.seed(42) # Startpunkt für den Zufallszahlengenerator festlegen. So hast du
# die gleichen Ergebnisse wie ich.
# Zielvektor für 20 zufällig verteilte Punkte erstellen:
x <- rep(NA, 20)
# Ersten Wert reinschreiben:
x[1] <- runif(1, 10, 80) # Zufallszahlen gleichmäßig zwischen 10 und 80 verteilt
x
count <- 1 # Anzahl Versuche, die benötigt werden
# 19 Punkte hinzufügen
for(i in 2:20) # Code zur Übersichtlichkeit einrücken!
{
x[i] <- runif(1, 10, 80) # random aus uniform-verteilung
count <- count+1
# Wenn eine minimale Distanz von 1.7 zu den anderen Punkten nicht gegeben
# ist, wird der Punkt ersetzt, also nochmal zufällig gezogen
while( min( abs(x[i] - x[-i]), na.rm=T ) < 1.7 )
{
x[i] <- runif(1, 10, 80)
count <- count+1
} # Ende while-Schleife | schreibe ich bei mehreren } oft dazu, damit
} # Ende for-Schleife | ich schnell sehen kann, wozu es gehört.
x
count # 45 versuche, 20 Punkte zwischen 10 und 80 mit dist > 1.7 zu verteilen.
plot(x, rep(1,20))
# Jetzt sind die Punkte zufällig gesetzt, aber immer mit einer minimalen Distanz
# zueinander. Dies ist für Versuchsplanungen (sampling design) interessant.
# Um Methoden der statistischen Inferenz (Schluss von einer Stichprobe auf die
# Grundgesamtheit) anwenden zu können, müssen die Stichpr. zufällig gezogen sein
points(runif(80, 10, 80), rep(c(0.6, 0.8, 1.2, 1.4),each=20), col=2, pch=17)
# Liegen mitunter sehr geklumpt und decken weniger gut den ganzen Bereich ab.
# 10.5 Geschwindigkeit for-Schleifen -------------------------------------------
# Generell sind Schleifen in R in der Berechnung um ein Vielfaches langsamer als
# vektorisierte Rechnungen. Bei größeren Datensätzen (> 500 Zeilen zB) solltest
# du sie nach Möglichkeit vermeiden und lieber apply (Abschnitt 10.3) verwenden.
# Sind sie trotzdem nötig, macht es immer Sinn, erst den gesamten Vektor mit
# NA oder 0 Werten zu erstellen, und diese dann einzeln zu füllen. Denn sonst,
# wenn ein leerer Vektor mit sovielen Elementen gefüllt wird, muss er nach jedem
# Schleifendurchlauf neu gespeichert werden --> "cannot allocate"-Fehler.
# Beweis: Vektorisieren ist besser!
x <- rnorm(5e6) # 5 Mio Zufallszahlen aus der Standardnormalverteilung
system.time( for (i in seq(along=x)) x[i] <- x[i]^2 ) # 11.42 Sekunden bei mir
system.time( x <- x^2 ) # vektorisierte Alternative # 0.03 Sekunden!
# Beweis: Vektor mit richtiger Länge vorbereiten ist besser!
x <- rnorm(1e5)
y <- numeric() # Leerer Vektor der Länge Null
system.time( for (i in seq(along=x)) y[i] <- x[i]^2 ) # 12.58 Sek
y <- numeric(1e5) # Leerer Vektor mit 0.1 Mio Elementen
system.time( for (i in seq(along=x)) y[i] <- x[i]^2 ) # 0.22 Sek
system.time( y <- x^2 ) # 0 Sek
# Merke: Schleifen sind langsam in R! (in z.B. FORTRAN oder C nicht)
# - Vermeide lange Schleifen durch Vektorisierung
# - für Informatiker und Nerds: Code-Auslagerung "Calling C and Fortran from R":
browseURL("http://cran.r-project.org/package=Rcpp")
browseURL("http://users.stat.umn.edu/~geyer/rc/")
# - Oder lerne, apply Funktionen anzuwenden. Die sind schnell und oft ganz toll,
# aber leider für Anfänger nicht immer leicht zu lernen.
# Da sich das aber wirklicht lohnt, soll das nächste Kapitel (11) dabei helfen.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 11: apply-Funktionen ###############################################
# Motiviert durch niedrige Geschwindigkeit for-Schleifen (Kapitel 10.5)
# apply und Co sind außerdem eleganter!
# 11.1 die Anwende(apply)-Funktionen in Übersicht
# 11.2 apply: "anwenden"
# 11.3 tapply: tagged-apply (pr Kategorie)
# 11.4 lapply: list-apply und sapply: simplify-apply
# 11.5 mapply: multiple-apply
# 11.6 Rechenzeit for, lapply, apply
# 11.7 Parallel Computing unter Windows
# 11.8 Übungsaufgabe
# 11.1 die Anwende(apply)-Funktionen in Übersicht ---Argumente----------------
?apply # Spalten- oder zeilenweise Funktion ausführen X, MARGIN, FUN
?tapply; ?by # table/tag-apply in Gruppen (Tags) X, INDEX, FUN, simplify
?aggregate # tapply für data.frames statt nur Vektoren x, by, FUN, simplify
?lapply # list-apply X, FUN
?sapply # simplify-apply, sonst wie lapply X, FUN, simplify
?mapply # multiple-apply FUN, MoreArgs, SIMPLIFY
# Leider sind die Argumente recht unterschiedlich benannt und geschrieben...
# Daher lasse ich sie oft unbenannt. In der Default-reihenfolge geht das ja.
# Eine sehr stark verkürzte Variante dieses Kapitels (auf englisch):
browseURL("http://stackoverflow.com/a/7141669/181638")
# Auch diese Seite finde ich ziemlich gut:
browseURL(
"http://nsaunders.wordpress.com/2010/08/20/a-brief-introduction-to-apply-in-r")
# 11.2 apply: "anwenden" -----------------------------------------------------
m <- matrix(sample(1:100, 24), ncol=6)
rownames(m) <- letters[1:4] ; colnames(m) <- LETTERS[21:26]
m
apply(X=m, MARGIN=2, FUN=mean) # mean spaltenweise (MARGIN=2) anwenden
apply(X=m, MARGIN=1, FUN=max) # max zeilenweise (MARGIN=1) anwenden --> vector
apply(X=m, MARGIN=1, FUN=range) # --> matrix
colMeans(m) # für einfache Funktionen gibt's die auch direkt
rowSums(m)
# Übung: Erstelle eine Tabelle mit mu und sd von m, mit Zeilen- und Spaltennamen
# Anwendungsbeispiel apply: Array
# Wirtschaftsdaten mehrerer Länder, Jahre, und Sektoren
load("Daten/11_econ.Rdata") # fiktive Daten aus Dominik Reusser's R-Kurs 2011
# falls bisher nicht runtergeladen:
load(url("http://dl.dropbox.com/u/4836866/Rclick/Daten/11_econ.Rdata"))
# Siehe Kap. 10.3
ls()
econ
str(econ)
class(econ) # array, also matrix mit (hier) 3 Dimensionen - gut für regelmäßige
mode(econ) # numeric, also alles mit Zahlen Zeitreihen geeignet
# Indizierung von Arrays
dimnames(econ)
econ["DE", ,] # Daten von Deutschland anzeigen
econ[,"Industrie",] # Industrie-Daten
econ[,,"2003"] # Daten aus 2003, Analog:
econ[,,2] # zweites Element der dritten Dimension (hier: "2002")
# Mittelwerte über Länder, Jahre, und Sektoren
# Das MARGIN-Argument gibt an, welche Dimension behalten wird:
apply(econ, MARGIN=c(1, 2), FUN=mean)
apply(econ, 2:3, mean) # bei Angabe in richtiger Reihenfolge können
apply(econ, 1, mean) # Argumentnamen weggelassen werden
apply(econ, 2, mean)
apply(econ, 3, mean)
apply(econ, 3, quantile)
# Siehe auch:
browseURL("http://de.wikibooks.org/wiki/GNU_R:_apply")
# 11.3 tapply: tagged-apply ---------------------------------------------------
data() # Zur Verfügung stehende Datensätze
head(chickwts, 15) # Hühnerwachstum abhängig von Fütterungstyp
tapply(X=chickwts$weight, INDEX=chickwts$feed, FUN=mean)
# Der Datensatz X wird anhand des INDEX gruppiert, bevor auf jede Gruppe die
# Funktion FUN angewendet wird. Die Ausgabe ist vom Typ array, mit Namen.
# t-apply steht ggf für tagged apply (the tags identify the subsets)
sort(tapply(chickwts$weight, chickwts$feed, mean))
# Für histogramme pro Gruppe siehe # 11.5 mapply
# Übung: schaue via tapply nicht nur mean, sondern auch min und max an
# Eine Funktion, die je zwei Werte zurück gibt, wird als array ausgegeben:
tapply(chickwts$weight, chickwts$feed, range)
futtergew <- tapply(chickwts$weight, chickwts$feed, range)
View(futtergew) # Objekts wird dabei intern in ein data.frame umgewandelt
sapply(futtergew, I) # macht ein echtes Dataframe draus
# sapply ist in 11.4 und 11.5 erklärt und wendet hier die Funktion I auf
# jedes Element der Liste an. I gibt im Wesentlichen den jeweiligen Input ohne
# Attribute zurück.
class(tapply(chickwts$weight, chickwts$feed, mean))
class(tapply(chickwts$weight, chickwts$feed, range))
# Anwendungsbeispiel tapply: iris-Daten
head(iris) # Größen von Kelch- und Kronenblätter dreier Schwertlilienarten
?iris
str(iris)
tapply(X=iris$Petal.Length, INDEX=iris$Species, FUN=mean)
tapply(iris$Petal.Length, iris$Species, max)
summary(iris)
tapply(iris$Petal.Length, iris$Species, summary)
# Für Tabellen verwendet man "by", annstatt INDEX heißt das Argument da INDICES
by(iris[,1:4], INDICES=iris$Species, FUN=sum)
# mit mean funktioniert das nicht so wie bei Neil Saunders.
# ToDo: rausfinden, warum nicht.
# Für ganze Datensätze aggregate verwenden (tapply nimmt als X nur einen Vektor)
state.x77 # Info über USA-Staaten
state.region # Regionen der Staaten der USA
cbind(state.x77, state.region) # schön sichtbar, dass die Regionen factors sind,
# die intern mit ganzen Zahlen beschrieben werden
aggregate(state.x77, by=list(Region = state.region), FUN=mean)
# Datensätze nach Kriterium auftrennen:
tapply(chickwts$weight, chickwts$feed, function(x) x) # Funktion gibt input
# geht eleganter mit # unverändert zurück
tapply(chickwts$weight, chickwts$feed, I)
# und noch eleganter mit
split(chickwts$weight, chickwts$feed)
# Anwendungsbeispiel tapply: Gruppierung nach mehreren Faktoren
head(warpbreaks)
str(warpbreaks)
# mittlere Anzahl Garn-Brüche (X) beim Weben nach Fadenspannung (INDEX):
tapply(X=warpbreaks$breaks, INDEX=warpbreaks$tension, mean)
# mittlere Anzahl Brüche nach Wolle-typ:
tapply(warpbreaks$breaks, warpbreaks$wool, mean)
# mittlere Anzahl Brüche nach Spannung und Wolle-typ:
tapply(warpbreaks$breaks, warpbreaks[,2:3], mean)
# Anwendungsbeispiel tapply: "gleitendes Mittel" (moving average)
# 10'000 äquidistante Daten. Mittelwert alle 100 Punkte.
x <- cumsum(rnorm(10000)) # hypothetische Daten
plot(x, type="l", col=8)
mx <- tapply(x, rep(1:100, each=100), mean)
lines(seq(0,10000, leng=100), mx, col=2, lwd=2)
wid <- 100 # Breite des gleitenden Fensters fexibler
tapply(x, rep(1:(length(x)/wid), each=wid), mean)
# Andere Lösung:
rowMeans(matrix(x, ncol=100, byrow=TRUE))
# Ich bin zoo-Gegner, ich sehe den Nutzen nicht. Überzeuge mich doch 🙂
install.packages("zoo")
require(zoo)
?rollapply # from package zoo
z2 <- zoo(round(rnorm(15),1)); z2
rollapply(z2, 3, mean, by = 3) # means of nonoverlapping groups of 3
# PS: ein echtes gleitendes Fenster hätte Überlappungen
# Hab deswegen eine rollapply Äquivalente für normale Daten in Kap 12 (movAv)
# 11.4 lapply: list-apply und sapply: simplify-apply ---------------------
bart <- list(1:10, 2:30, 3:40)
bart
lapply(bart, mean) # list
sapply(bart, mean) # vector (s für simplify list to vector)
# lapply statt for-Schleifen: Spalten gegeneinander plotten
dat.a <- as.data.frame(matrix(runif(20),ncol=4))
dat.a
dat.b <- as.data.frame(matrix(runif(20),ncol=4))
par(mfrow=c(2,2))
for(i in 1:4) plot(dat.a[,i], dat.b[,i], pch=16, main=i)
# lange for-Schleifen ersetzen mit lapply (rechnet um ein vielfaches schneller):
dummy <- lapply(1:4, function(i){plot(dat.a[,i], dat.b[,i], pch=16, main=i)} )
dummy # Hier nicht benötigter Kram, daher auch dummy genannt
graphics.off()
# Rechenzeit:
dat <- matrix(rnorm(5e6), ncol=50)
dim(dat) # 100'000 Zeilen, 50 Spalten
Erg <- NA # Vektor mit Ergebnis je Zeile: Wert genau zwischen min und max
system.time( for(i in 1:100000) Erg[i] <- mean(range(dat[i,])) )
rm(Erg) # 53 Sekunden (Es lebe SSD! Bei alten PCs mehr)
system.time( Erg <- sapply(1:1e5, function(i) mean(range(dat[i,])) ) )
# 3.7 Sekunden = 7% der Rechenzeit der Schleife !!
hist(Erg)
# Man kann also (nichtrekursive) Schleifen ganz einfach schneller machen:
# Einfach "sapply(1:n, function(i)" und ")" drumherum schreiben.
# Und man muss vorher keine leeren Objekte erstellen.
Erg <- rep(NA, 100000) # auch 3.7 Sekunden
system.time( for(i in 1:100000) Erg[i] <- mean(range(dat[i,])) )
# weiteres Beispiel - lapply und sapply im Vergleich:
x <- list(a = 1:10, beta = exp(-3:3), logic = c(TRUE,FALSE,FALSE,TRUE))
x
dim(x); class(x) # list mit 3 elementen
length(x) # jedes dieser Bestandteile hat selbst wieder eigene Elemente
lapply(x, mean)
sapply(x, mean) # simplify list to vector
# Anwendungsbeispiel sapply: Quantile
set.seed(2) # legt einen Startpunkt für den Zufallszahlengenerator fest, dadurch
# sind die reproduzierbar: du kriegst die gleichen "Zufalls"zahlen wie ich.
Ksat <- round(rexp(n=45, rate=0.02),1) # Random EXPonentially distributed
quantile(Ksat)
# 25% der Daten ist kleiner als 16.3
# 100%-75% = 25% der Daten ist größer als 73.7 (75% Quantil)
# 9-Zahlenmaß: Median, Quartile (Fourth), Oktile (Eighth), Hexadezile, Range
# Quantile kann verschiedene Algorithmen für die Berechnung verwenden (type)
# Um das schneller zu vergleichen, schreiben wir eine Funktion:
myquant <- function(tp) quantile(Ksat,prob=c(0,1,2,4,8,12,14,15,16)/16, type=tp)
# Die führt dann jeweils abhängig von tp quantile von Ksat aus
myquant(tp=1) # Mehr zu Fuktionen in Kap 12. Wenn du nachvollziehen kannst,
myquant(3) # was hier passiert, reicht das aber erstmal.
# insbesondere bei den hohen Quantilen ist das Ergebnis abhängig vom type!
# (weil dort in diesem Datensatz weniger Werte sind)
# Vergleich in einer Tabelle:
quants <- sapply(X=1:9, FUN=myquant); quants # simplify-apply
# deutlich eleganter (und schneller) als mit einer for-Schleife.
boxplot(t(quants), at=c(0:2,4,8,12,14:16)/16, boxwex=0.1, xlim=0:1,
las=1, horizontal=TRUE)
# Methoden unterscheiden sich stark wenn wenig Werte da sind...
lines(ecdf(Ksat), verticals=T, do.points=F, col=2, lwd=3)
# Anderes Beispiel für Ökonomen:
# Fake Data: Einkommen von Leuten in verschiedenen Einkommensklassen
set.seed(42) # Startpunkt für Zufallszahlengenerator (zur Reproduzierbarkeit)
Einkommen <- round(c(rexp(20)*1e4, rexp(20)*1e5, rexp(20)*1e6))
Klasse <- rep(c("1.low", "2.mid", "3.high"), each=20)
# Funktion definieren, die die benötigten Kennwerte rausgibt:
dezile <- function(input) quantile(input, probs=c(0.1, 0.9) )
dezile(Einkommen) # funktioniert schonmal
# Jetzt nach Klassen getrennt:
tapply(Einkommen, INDEX=Klasse, FUN=dezile) # list Ausgabe unübersichtlich
sapply(X=tapply(Einkommen, INDEX=Klasse, FUN=dezile), FUN=I)
# lapply kann eine "list" zu einem "data.frame" umwandeln, siehe l2df im
# Paket BerryFunctions
install.packages("berryFunctions"); library("berryFunctions")
?l2df
# Ähnlich wie lapply funktioniert auch do.call:
eg3 <- list(KA=data.frame(SW=1:6, SB=2:-3), LS=data.frame(SW=23:24, SB=c(-3,2)))
eg3
do.call(rbind, eg3)
# 11.5 mapply: multiple-apply ------------------------------------------------
check.math <- function(a,b,c)
cat(paste(a,"+",b,"=",c, ifelse(a+b==c,"ist richtig","ist falsch"),"\n"))
check.math(3,4,7)
check.math(4,5,7)
alang <- c(6,2,4,5)
blang <- c(5,3,5,9)
clang <- c(3,5,8,7)
dummy <- mapply(check.math, alang, blang, clang)
# Histogramme einzelner Gruppen mit Überschrift
par(mfrow=c(2,3))
mapply(hist, split(chickwts$weight, chickwts$feed), main=levels(chickwts$feed),
MoreArgs=list(xlim=range(chickwts$weight), col=2, breaks=1:10*50))
# Für ProgressBars in *apply-Funktionen siehe:
#browseURL(paste("http://ryouready.wordpress.com/2010/01/11/",
# "progress-bars-in-r-part-ii-a-wrapper-for-apply-functions", sep=""))
install.packages("pbapply")
library("pbapply")
dummy <- pbsapply(rep(0.02, 100), Sys.sleep)
# 11.6 Rechenzeit for, lapply, apply -------------------------------------------
n <- 1e3
dat <- matrix(rnorm(n*50), ncol=50)
library("microbenchmark")
mr_bad_for <- function() {Result <- NA
for(i in 1:n) Result[i] <- mean(range(dat[i,])) ; Result}
mr_good_for <- function() {Result <- rep(NA, n)
for(i in 1:n) Result[i] <- mean(range(dat[i,])) ; Result}
mr_sapply <- function() {Result <- sapply(1:n, function(i) mean(range(dat[i,])) ); Result}
mr_apply <- function() {Result <- apply(dat, MARGIN=1, function(x) mean(range(x))); Result}
system.time(mr_good_for())
system.time(mr_good_for())
# will be slightly different each time
# Folgende Berechnung dauert relativ lange! - n ggf. runtersetzen
speed <- microbenchmark(mr_bad_for(), mr_good_for(), mr_sapply(), mr_apply(), times=20)
speed
plot(speed) # if there is one extreme outlier, run it again
par(mar=c(3,7,2,1))
plot(speed, horizontal=TRUE, las=1)
library("ggplot2")
autoplot(speed)
# microbenchmark for other tasks:
library("microbenchmark")
a <- sample(0:1, size=3e6, replace=TRUE)
speed <- microbenchmark(a != 0, ! a == 0, times=50)
boxplot(speed, notch=TRUE, unit="ms", log=F)
# 11.7 Parallel Computing unter Windows ----------------------------------------
# Motivation: kürzere Rechenzeit ein Literaturbeispiel von Tausenden:
browseURL("www.fernuni-hagen.de/FACHSCHINF/1727/Zusammenfassung.pdf")
# Für UNIX-Systeme gibt es die Pakete parallel und multicore - Für Windows snow
if(!require(snow)) {install.packages("snow") ; require(snow) }
cl <- makeCluster(4) # 4 Kerne auf dem Rechner
# Ggf. vorher notwendig: Objekte an Cluster übergeben:
clusterExport(cl, c("centroide","patchsize","distance"))
browseURL("http://stackoverflow.com/questions/10095956")
# Nimm eine kleinere Zahl, falls der PC für dieses Beispiel zu langsam ist:
sinnlos <- function(i) quantile(rnorm(1e8))
system.time( sinnlos(1) ) # 12 sek
system.time( sapply(1:10, sinnlos) ) # 131 sek
system.time(parSapply(cl, 1:10, sinnlos) ) # 42 sek
# Nie vergessen:
stopCluster(cl)
# Noch mehr RAM usage leeren:
gc()
# Das sind die Speicherfresser
sort(sapply(ls(), function(x) object.size(get(x)))/1e6, dec=T) # in MB
# 11.8 Übungsaufgabe -----------------------------------------------------------
# Gegeben sind folgende (fiktive) Daten aus unregelmäßigen Messkampagnen:
UebDat <- read.table("Daten/11_applytest.txt", header=T)
# oder
File <- "http://dl.dropbox.com/u/4836866/Rclick/Daten/11_applytest.txt"
UebDat <- read.table(url(File), header=T)
UebDat
# Diese Zeitreihen sollen so als Boxplots dargestellt werden, dass das Skript
# auch mit weiteren Daten laufen wird. Nutzungstypen sind zu berücksichtigen.
# Zudem will man für jede Nutzung und jedes Jahr den Mittelwert mit einzeichen.
# Eine mögliche Lösung findet sich weiter unten, aber probiere doch erst selbst.
# Bessere oder elegantere Lösungen sind - wie immer - willkommen!
View(UebDat)
# Parameter fürs Plotten ermitteln
Jahre <- min(UebDat$Jahr):max(UebDat$Jahr)
nJ <- length(Jahre) # Anzahl Jahre
Range <- range(UebDat$Wert) # Wertebereich
Range[2] <- Range[2] + 0.05*diff(Range) # zweiter Wert höher für Anzahlangabe
# Daten nach Landnutzung auftrennen
Dat <- split(UebDat, UebDat$Nutzung)
class(Dat)
str(Dat)
nN <- length(Dat); nN # 2, number of Nutzungen
# Nutzungen je in einem eigenem Fenster darstellen
par(mfrow=c(nN,1), mar=c(2,3,1,1), mgp=c(3,0.5,0))
for(i in 1:nN)
{
# evtl NAs in fehlenden Jahren hinzufügen -> Bereich auf allen Achsen gleich
fehl <- Jahre[!Jahre %in% Dat[[i]][,"Jahr"]]
toplot <- Dat[[i]][2:3]
if(length(fehl)>0) toplot<-rbind(Dat[[i]][2:3], data.frame(Jahr=fehl,Wert=NA))
# Das eigentliche Plotten
posit <- boxplot(tapply(toplot[,2], toplot[,1], function(input) input),
ylim=Range, main=names(Dat)[i], las=1, col=sample(colors(), nJ))
points(1:nJ, tapply(toplot[,"Wert"], toplot[,"Jahr"], mean), pch=8)
text(1:nJ, posit$stats[5,], paste("n =", posit$n), adj=c(0.5,-0.2))
rm(fehl, toplot, posit, i)
}# Ende for-Schleife
# Das alles ist jetzt unproblematisch nochmal für folgende Daten möglich:
set.seed(8)
UebDat <- data.frame(Nutzung=sample(c("Wald","Wiese","Acker"), 90, rep=T),
Jahr=sample(2005:2013, 90, rep=T), Wert=sample(100:300, 90, rep=T) )
# entsprechende Zeilen oben nochmal abschicken...
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 12: Funktionen #####################################################
# 12.1 Funktionseinführung
# 12.2 Einfache Beispiele
# 12.3 Verweis auf Paket berryFunctions
# 12.4 Weitere mögliche Argumente
# 12.5 Lexical scoping - wo kommen die Variablen her?
# 12.6 Input prüfen ist wichtig (Exception handling)
# 12.7 Beispiele für praktische Funktionen
# 12.7.1 Rechteck manuell in Graphik zeichnen
# 12.7.2 Gleitendes Mittel
# 12.7.3 Fonttester
# 12.7.4 SCS Verfahren zur Bestimmung des Abfussbeiwertes
# 12.1 Funktionseinführung -----------------------------------------------------
# Funktionen schreibt man oft, wenn man Code mehrfach ausführen muss/will, oder
# wenn man Routinen für andere Leute schreiben will.
# Generelle Struktur:
Funktionsobjektname <- function(argument1, argument2, ...) {"Rechenanweisungen"}
?"function"
# 12.2 Einfache Beispiele ------------------------------------------------------
# Daten sollen geplottet und mit 7 multipliziert zurückgegeben werden:
myfunct <- function(grappig)
{
plot(grappig, type="l")
return(grappig*7)
}
# grappig ist Niederländisch und symbolisiert die freie Benennung der Argumente.
myfunct( c(5,11,3,7) )
# return() # gibt das was drin steht zurück, und danach wird die Funktion
# nicht weiter ausgeführt! Sollte also erst am Ende stehen.
# Es muss nicht gesetzt werden: die letze Anweisung ("expression") wird auch
# ohne return zurück gegeben:
myfunct <- function(grappig) {plot(grappig, type="l"); grappig*7 }
myfunct( c(5,11,3,7) )
# Nullstellen (P-Q-Formel) linearer Funktionen y = x^2 + px + q
pq <- function(p,q)
{
w <- sqrt( p^2 / 4 - q )
c(-p/2-w, -p/2+w)
} # Ende der Funktion
# Code einrücken erhöht die Menschenlesbarkeit enorm!
# Funktion verwenden
pq(3, -12) # -5.274917 2.274917
# Gleich mit Plotten im angegebenen Bereich
pq <- function(p,q,a=-5, b=5) # Standardwerte mit "=" vorgeben (defaults)
{x <- seq(a, b, len=100)
plot(x, x^2+p*x+q, type="l"); abline(h=0, v=0)
w <- sqrt( p^2 / 4 - q )
c(-p/2-w, -p/2+w)
}
pq(3, -12) # a und b haben bereits Standardwerte, müssen nicht angegeben sein.
pq(3, -8, a=-7) # können aber (auch einzeln) geändert werden
pq(3, -8, , 3)
# Nicht namentlich genannte Argumente müssen in der richtigen Reihenfolge sein,
# bei Änderung der Reihenfolge muss man sie aber angeben:
pq( a=-2, q=-12, p=3)
# Funktionen werden zudem oft in Solvern benötigt, zB um ein Minimum finden:
optimize( f=function(x) x^2 + 3*x - 12 , interval=c(-5, 5) )
# Siehe Kapitel 18: Optimierung
# Wurzel aus Skalarprodukt zweier Vektoren
?"%*%"
wsk <- function(x, y) drop( sqrt(x%*%y) ) # drop macht aus matrix ein vector
wsk(1:4, 4:1)
wsk(c(5.2,pi,-3), c(7,-4,7.5))
# Barcodes erstellen
barcode <- function() { # einfache Verwendung von function: Code zusammenfassen
par(mar=c(0,0,0,0)); layout(matrix(1,1,1), widths=lcm(15), heights=lcm(9))
t4 <- data.frame(a=1:500, b=ifelse(rnorm(500)>0.9, 1, 0))
plot(t4, type="n", xaxt="n", yaxt="n" )
abline(v=t4[t4$b==1,1], lwd=3); box(col="white") }
windows(record=T)
for (i in 1:20) barcode()
# mit PG up/dn (Bild) kann man die Bilder durchklicken
# Erste und letze Zeilen einer Tabelle
flast <- function(dat, n=3) return(dat[c(1:n,(nrow(dat)-n+1):nrow(dat)),])
flast( data.frame(letters, rt(26, 5)) )
# eleganter mit: head, tail
flast <- function(dat, n=3) rbind( head(dat, n), tail( dat, n) )
flast( data.frame(letters, rt(26, 5)) )
# Interferenz durch Pixelabstände
x <- seq(-pi,pi,len=1000)
par(mar=rep(0,4))
coola <- function (a) plot(a*sin(x), a*cos(x), asp=1, type="l")
coolb <- function (a) lines(a*sin(x), a*cos(x))
coola(10)
for(i in 5:99/10) coolb(i) # Zieh mal die Graphik kleiner und größer...
# In einer Funktion auf Objekte zugreifen: get()
obj_1 <- cumsum(rnorm(1000))
obj_2 <- cumsum(1:10*c(-2,3))
testf <- function(k) plot( get(paste0("obj_", k)), type="l" )
testf(1)
testf(2)
# Integrale
f <- function(x) {-x*x +7}
plot(-10:40/10, f(-10:40/10), type="l")
abline(v=0, h=0)
f(1)
integrate(f, 0, 1.5)
lines(c(1.5,1.5,0,0:15/10), c(f(1.5),0,0,f(0:15/10)), col=2, lwd=2)
text(.75, 3, 9.4, col=2)
# 12.3 Verweis auf Paket berryFunctions ----------------------------------------
browseURL("http://cran.r-project.org/web/packages/berryFunctions/index.html")
# Dieses Paket hat unter anderem auch eine Menge meiner Helper-Functions
install.packages("berryFunctions")
library(berryFunctions)
?berryFunctions # Feedback ist herzlich willkommen!
# 12.4 Weitere mögliche Argumente --------------------------------------------
testdots <- function(x,...) plot(x, ...) # dots expansion, three dots, ellipsis
testdots(3:8)
testdots(3:8, main="zwo") # oder was auch immer.
# to pass ... to a different function:
myfunc <- function(g, ...) {plot(g); list(...)}
myfunc(8, buh=6)
# Für jede Funktion eine eigene Argumentliste:
tp <- function(...) {plot(7:8, ...) ; legend("top", "Text hier", ...)}
tp()
tp(type="o") # legend hat kein type! die Lösung: berryFunctions::owa
# different Ellipsis passed to different functions, each with their own defaults
?owa # argument overwriting and combining. Read the examples.
# 12.5 Lexical scoping - wo kommen die Variablen her? --------------------------
x <- 2
ls() # x ist ein Objekt (Variable) im Global Environment, dem user workspace
f <- function(y) c(x,y)
f(7) # x wird aus dem globalenv() gezogen, es gibt kein x innerhalb der Funktion
f <- function(y) {x <- 3; print(environment()); c(x,y)}
f(7) # hier kommt x aus dem funktionseigenen temporären environment
x # im globalenv hat x nach wie vor die Zahl 2
# Das ist bei geschachtelten Funktionen komplizierter. Für Fortgeschrittene:
browseURL("http://andrewgelman.com/2014/01/29/stupid-r-tricks-random-scope/")
# Für Anfänger eher geeignet, mit Übungsaufgaben:
browseURL("http://adv-r.had.co.nz/Functions.html")
tf <- function(x) {if(x)buh=8; list(exists("buh", inherits=FALSE), buh)}
tf(T)
tf(F) # geht nicht, buh existiert nicht
buh=7
tf(T)
tf(F) # buh existiert nicht im Function environment, aber in globalenv()
# Nur einmal pro Sitzung warnen/informieren
# r.789695.n4.nabble.com/Display-warning-only-once-in-session-td4696005
myFunc <- local({ # (c): Bill Dunlap, TIBCO Software
notWarnedYet <- TRUE
function(x) {
if (notWarnedYet) {
warning("myFunc is funky")
notWarnedYet <<- FALSE # note use of <<-
}
sqrt(log2(x))
}
})
myFunc(1:5)
myFunc(1:5)
myFunc(1:5)
# 12.6 Input prüfen ist wichtig (Exception handling) ---------------------------
# z.B. sollen in folgender Berechnung der Fläche eines planaren Dreiecks
# 3-stellige Vektoren mit je den x und den y-Koordinaten der Ecken eingehen.
# Bei anderen Eingaben soll mit Fehlermeldung abgebrochen werden.
dreieck <- function(x, y, digits=3)
{
if( !is.vector(x) | !is.vector(y) ) stop("Input müssen Vektoren sein!")
if(length(x) != 3 | length(y) !=3 ) stop("Vektoren sind nicht 3-stellig")
A <- 0.5*(x[1] * (y[2]-y[3]) + x[2] * (y[3]-y[1]) + x[3] * (y[1]-y[2]))
return(round(abs(A),digits))
}
a <- c(1,5.387965,9); b <- c(1,1,5)
plot(a[c(1:3,1)], b[c(1:3,1)], type="l", asp=1)#; grid()
dreieck(a,b)
dreieck(a,b[1:2])
# Always try to make your functions foolproof. That way, if someone else uses it
# wrongly, they'll know what they did wrong. And by the way, this includes you
# yourself, three years later. (or 3 weeks, actually)
# cycle a vector, eg change 1 2 3 4 5 6 7 to 4 5 6 7 1 2 3
Cycle <- function(vector, shift)
{ # make nice indents to keep things readable
# make sure only vectors are supplied as input
if(!is.vector(vector)) stop("input has to be a vector")
# don't mess with shift: It has to be ONE integer, less than vectorlength
if(length(shift) > 1) stop("shift should be ONE number")
if(is.na(as.integer(shift))) stop("shift has to be an integer")
n <- length(vector)
if(shift >= n) stop("shift has to be smaller than length(vector)")
# Well, here's the actual thing happening:
vector[c((shift+1):n, 1:shift)]
} # end of function
# Test the function
Cycle(1:8, 5) # 6 7 8 1 2 3 4 5
Cycle(letters, 8)
is.integer(5) # FALSE # cannot be used as criterion in function
is.integer(5:6) # TRUE # see ?":"
typeof(5) # double
# Wat use is this function, you may wonder...
# read on 😉
# Draw polygon with n corners
n_gon <- function(n)
{ # function start
op <- par(mar=rep(0,4)) # don't mess with other people's graphics:
# assign the current par and reset it at the end of your function.
# Keep track of what's going on:
# the next line creates n vertices spaced equally along the unit circle
v <- data.frame(x=cos(1:n*2*pi/n), y=sin(1:n*2*pi/n))
# Set the plot, but don't plot anything
plot(v, asp=1, type="n")
# Here we go: the lines. see comment below about how I figured this out
segments(x0=v$x, y0=v$y, x1=Cycle(v$x, floor((n-1)/2)),
y1=Cycle(v$y, floor((n-1)/2)) )
# Set the old parameters again:
par(op)
} # end of function
n_gon(8)
windows(record=T)
for (i in 3:30) n_gon(i)
# repeatedly clicking 'page up' on this one is not recommended for epileptics:
for (i in 20:300) n_gon(i)
n_gon(50) # nice interference.
n_gon(100) # Look at the white bands which are not really existent...
# I guess this is nothing really new, but I was just wondering whether there is
# any regularity behind drawing concentric polygons. I found it is necessary to
# draw a line once from every vertex to the one that is (n-1)/2 vertices later.
# Then I realized, I need to recycle the vector to do that.
# And then I thought it nice to take this as an example of foolprofing functions
# Which is really called exception handling, but foolprofing sounds cooler.
# Useless as this was contentwise, I still hope you learned something...
# 12.7 Beispiele für praktische Funktionen -------------------------------------
# 12.7.1 Rechteck manuell in Graphik zeichnen ----------------------------------
rechteck <- function(...){
cat("Bitte wähle den Bereich im Graphikfenster aus.",
"Klicke zuerst linksoben, dann rechtsunten.\n", sep="\n")
flush.console() # Damit die Aufforderung in die Console geschrieben wird
# Ausgaben werden sonst gebuffert und erst am Ende der Funktion geschrieben
werte <- locator(2)
polygon(c(werte$x, rev(werte$x)), rep(werte$y,each=2), ...)
cat("Rechteck gezeichnet.\n",
"Danke Berry für diese tolle Funktion: berry-b@gmx.de.\n", sep="")
}
plot(rexp(300), rt(300,5))
rechteck()
rechteck(col=2)
# eine ähnliches Prinzip steckt hinter berryFunctions:::bzoom, siehe 12.7
rechteck(border=3, lwd=4)
rechteck(border=NA, col="white") # damit kann man auch ganz gut Text löschen
?polygon # für die Argumente, die gesetzt werden können.
loeschen <- function() {cat("linksobere, dann rechtsuntere Ecke klicken.\n")
flush.console() ; werte <- locator(2)
polygon(c(werte$x, rev(werte$x)), rep(werte$y,each=2), border=NA, col="white")
}
plot(rexp(300), rt(300,5))
loeschen()
# 12.7.2 Gleitendes Mittel -----------------------------------------------------
# Siehe auch BerryFunctions::movAv sowie decompose sowie zoo::rollaply
movav <- function(dat, steps=7)
{ # Funktionsanfang
# halbe Fensterbreite (in eine Richtung), dabei gerade Zahlen auffangen
# (soll immer ungerade sein, damit das Fenster eine definierte Mitte hat):
st <- round(steps/2-.1,0)
# Länge des Eingangsvektors:
ld <- length(dat)
# leeren Ergebnisvektor erstellen:
v <- vector("numeric", ld)
# Bereich für NA ermitteln (Hälfte des gleitenden Fensters am Anfang + Ende)
# bei Standard 7 steps erste und letzte drei Werte jeweils := NA
v[c(1:st, (ld-st+1):ld)] <- NA
# für den Rest jeweils das Mittel aus dem Fenster um den einzelnen Wert herum
for(i in (st+1):(ld-st)) v[i]= mean(dat[(i-st):(i+st)],na.rm=TRUE)
return(v) # Ergebnisausgabe
} # Funktionsende
set.seed(29); a <- runif(60, 5,50)
data.frame(a[1:20], movav(a[1:20]))
plot(a, type="o", pch=16, las=1) # sieht wirr aus
lines(movav(a), col=2, lwd=4) # bringt Übersicht rein, zeigt Trends
lines(movav(a,3), col=4, lwd=4)
lines(movav(a,15), col=3, lwd=4) # unterscheidlich starke Glättung
graphics.off()
# Weitere Beispiele in Dokumentation library(berryFunctions); ?movAv
# Visualisierung des Funktionsprinzips von movav
# dir.create("MovingAverage_Funktionsprinzip")
# png("MovingAverage_Funktionsprinzip/movav_%02d.png", 1024, 600)
set.seed(1); a <- runif(30, 5,50)
b <- movAv(a)
par(mar=c(3,3,1,1), cex=2)
plot(a, type="o", pch=16, las=1, ylim=c(5,55))
#blines(1:30,a)
windows(record=T)
for (i in 4:27) {
plot(a, type="o", pch=16, las=1, ylim=c(5,55))
#blines(1:30,a)
polygon(x=c(i-4,i-4,i+3,i+3)+0.5, y=c(5,55,55,5), col=rgb(1,0,0, alpha=0.2),
border=2)
#blines(1:(i-1), b[1:(i-1)], type="o", pch=16, col=2, lwd=3)
points(i, b[i], pch=16, col=2)
}
# mit PgUp und PgDn (Bild hoch und runter)-Tasten Bildersequenz ablaufen lassen
windows()
plot(a, type="o", pch=16, las=1, ylim=c(5,55))
lines(1:30, b, type="o", pch=16, col=2, lwd=3)
# dev.off()
# 12.7.3 Fonttester ------------------------------------------------------------
fonttester <- function(f){ par(mar=c(0,0,0,0))
plot(0:100, rep(0,101), type="n", ylim=c(-0.3,2.8)) ; cx <- 1.5
text(0,2.5, "A B C D E F G H I J K L", adj=0 ,cex=cx)
text(33,2.5, "M N O P Q R S T U V W X", adj=0 ,cex=cx)
text(66,2.5, "Y Z ! ' § $ % & / ( ) =", adj=0 ,cex=cx)
text(0,2, "a b c d e f g h i j k l", adj=0 ,cex=cx)
text(33,2 , "m n o p q r s t u v w x", adj=0 ,cex=cx)
text(66,2 , "y z 1 2 3 4 5 6 7 8 9 0", adj=0 ,cex=cx)
text(0,1.5, "ß ´ } ] [ { \U00B3 \U00B2 @ \U20AC ' ~", adj=0 ,cex=cx)
text(33,1.5, "* ü ä ö Ü Ä Ö ? # ` ^ °", adj=0 ,cex=cx)
text(66,1.5, "< > | ; : . , - _ + @ \U20AC", adj=0 ,cex=cx)
abline(h=1.25) ; abline(v=c(31,64))
text(0,1 , "A B C D E F G H I J K L", adj=0,font=f ,cex=cx)
text(33,1 , "M N O P Q R S T U V W X", adj=0,font=f ,cex=cx)
text(66,1 , "Y Z ! ' § $ % & / ( ) =", adj=0,font=f ,cex=cx)
text(0,0.5, "a b c d e f g h i j k l", adj=0,font=f ,cex=cx)
text(33,0.5, "m n o p q r s t u v w x", adj=0,font=f ,cex=cx)
text(66,0.5, "y z 1 2 3 4 5 6 7 8 9 0", adj=0,font=f ,cex=cx)
text(0,0 , "ß ´ } ] [ { \U00B3 \U00B2 @ \U20AC ' ~", adj=0,font=f ,cex=cx)
text(33,0 , "* ü ä ö Ü Ä Ö ? # ` ^ °", adj=0,font=f ,cex=cx)
text(66,0 , "< > | ; : . , - _ + @ \U20AC", adj=0,font=f ,cex=cx)
text(102,2.5, paste("font =",f), adj=1, font=2)
text(98,1.8,"ß", font=5, cex=3) }
fonttester(5)
png("fonts%02d.png", width=1000, height=200, units="px")
for(i in 1:40) fonttester(i)
dev.off()
# 12.7.4 SCS Verfahren zur Bestimmung des Abfussbeiwertes ----------------------
Neff <- function(N) { CN <- 0:100; neff <- (N-5080/CN+50.8)^2/(N+20320/CN-203.2)
par(mar=c(3.2,3,3,3.5))
plot(CN, neff, xlim=c(CN[which(neff<=N)[1]],100),ylim=c(0,N),
main="EffektivNiederschlag in Abhängigkeit des CN-Faktors",
type="l", las=1, mgp=c(2,1,0), yaxs="i", xaxs="i")
legend("topleft", paste("Niederschlag =", N, "mm"), bg="white")
axis(4,seq(0,N,len=11)[c(1,3,5,7,9,11)],c(0,2,4,6,8,10)*10, col.axis=2, las=1)
mtext("Abflussbeiwert Psi",4,2, col=2)}
Neff(N=15)
windows(record=T)
plot(1, type="n", xlab="CN", ylab="Neff [mm]", xlim=c(0,110), ylim=c(3,100),
main="EffektivNiederschlag in Abhängigkeit des CN-Faktors", las=1)
CN <- 1:100
for (N in 1:10*10) { neff <- (N-5080/CN+50.8)^2/(N+20320/CN-203.2)
lines(CN[which(neff<=N)],neff[which(neff<=N)], col=N/10) }
for (i in 1:10*10) { text(101,i, paste(i,"mm"), col=i/10, adj=0) }
text(95,100, adj=1, "bei Niederschlag =") ; text(40,50,"Neff = ", adj=1)
text(40,50,"(N-5080/CN+50.8)^2 /\n(N+20320/CN-203.2)", adj=0)
#
plot(1, type="n", xlab="CN", ylab="Neff [mm]", xlim=c(0,110), ylim=c(3,100),
main="EffektivNiederschlag in Abhängigkeit des CN-Faktors", las=1)
CN <- 1:100
for (N in 1:10*10) { neff <- (N-1270/CN+12.7)^2/(N+24130/CN-241.3)
lines(CN[which(neff<=N)],neff[which(neff<=N)], col=N/10) }
for (i in 1:10*10) { text(101,i, paste(i,"mm"), col=i/10, adj=0) }
text(95,100, adj=1, "bei Niederschlag =") ; text(40,50,"Neff = ", adj=1)
text(40,50,"(N-1270/CN+12.7)^2 /\n(N+24130/CN-241.3)", adj=0)
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 13: Zeit und Datum #################################################
# 13.1 Datumsformat einlesen und formatieren
# 13.2 Zeitformat einlesen und formatieren
# 13.3 Plot und Achsenbeschriftung
# 13.4 Zeitreihen aggregieren
# 13.5 Hydrologisches Jahr
# 13.6 nicht-equidistante Zeitreihen
# 13.7 Übungsaufgaben
# 13.8 Lösungen Übungsaufgaben
# Beispieldaten erzeugen:
Daten <- read.table(as.is=T, header=T, text="
Datum Zeit Wert
21.01.2009 13:30:00 19
21.01.2009 13:45:00 28
21.01.2009 14:00:00 17
21.01.2009 14:15:00 29
21.01.2009 14:30:00 24
21.01.2009 14:45:00 34
21.01.2009 15:00:00 28
21.01.2009 15:15:00 12")
# 13.1 Datumsformat einlesen und formatieren -----------------------------------
zeit1 <- as.Date(Daten[,1], format="%d.%m.%Y")
# alternative Formate "%m-%d-%Y" "%Y/%m/%d" usw.
class(zeit1) ; as.numeric(zeit1) # intern ganze Zahlen für Tage seit 1970-01-01
as.Date("20090921","%Y%m%d") # geht auch ohne Leerzeichen
xyz <- c("2014-02","2014-04","2015-07","2015-11")
# Datumsreihen im Format mm.yyyy einlesen geht nicht mit
as.Date( xyz , "%Y-%m") # sondern nur mit
as.Date(paste0("01.",xyz), "%d.%Y-%m") # Datumsobjekte brauchen immer einen Tag
# locale-specific sind Sachen wie %A, die auf jedem Sys.getlocale() anders sind:
format(Sys.Date(), "%a %d.%b.%Y (= ein %A im %B)")
# Mehr Formate unter Details in der Hilfe von
?strptime
format(Sys.time(), " %H Uhr %M")
format(Sys.Date(), "%A, der %d. %b im Jahre %Y")
tenweeks <- seq(Sys.Date(), length.out=10, by="1 week") # next ten weeks
tenweeks
weekdays(Sys.Date())
months(tenweeks)
# Exakt 1 Jahr zum Datum dazu addieren:
k1 <- as.Date("15.3.2007", "%d.%m.%Y")
k1
k1 + 1
k1 + 365 # geht nicht bei Schaltjahren. workaround:
k2 <- paste(as.numeric(format(k, "%Y"))+1 , format(k, "%m-%d") , sep="-" )
k2 # es gibt aber mit Sicherheit eine einfachere Lösung...
k2 + 7 # geht noch nicht, denn:
mode(k2) # k2 ist kein Datum, sondern character
k3 <- as.Date(k2, "%Y-%m-%d")
k3
mode(k3) # zeigt numeric, denn intern ist ein Datum, wie bei Excel, eine Zahl
class(k3)
k3 + 7
# 13.2 Zeitformat einlesen und formatieren -------------------------------------
# Bei Uhrzeit (strptime bedeutet sowas wie string parse time)
zeit2 <- strptime(Daten[,2], format="%H:%M:%S") # format ist für input
# automatisch wird das aktuelle Datum hinzugefügt
class(zeit2) # Portable Operating System Interface = POSIX: internat. Standard
as.numeric(zeit2) # Ganze Zahlen für Sekunden seit 1970-01-01 01:00:00 CET
# Bei Datum und Uhrzeit in 2 Spalten: Spalten mit paste zusammenfügen
Zeit <- strptime( paste(Daten[,1],Daten[,2]) , "%d.%m.%Y %H:%M:%S")
data.frame(Zeit)
str(Zeit) # ein POSIX long time Objekt
# Abkürzungen: "%Y-%m-%d %H:%M" = "%F %R"
format(Sys.time(), "%F %T") # "2014-12-04 12:49:37"
# Ausgabe formatieren (z.B. für Achsenbeschriftungen)
format(Zeit, "%H:%M")
# alternativ lässt sich hier auch strftime (string format time) benutzen.
strftime(Zeit, "%d.%m., %H:%M:%S") #format ist für output
?strftime # für eine ganze Reihe von Formatvorlagen inkl Wochentag etc.
# Einzelne Komponenten aus Daten herausnehmen:
format(Zeit, "%M")
which( format(Zeit, "%M")==15 )
Daten[ format(Zeit, "%M")==15 , ]
# Zeitreihen
anfa <- as.Date("2007/07/18")
ende <- as.Date("2007/08/05")
anfa; ende
seq(anfa, ende, by="day") # berücksichtigt Schaltjahre
strftime(anfa, "%b")
strftime(anfa, "%B")
# Zeitobjekt in Spalte einfügen ging früher nicht bei POSIXlt
Daten$Zeit2 <- Zeit
Daten$Zeit2 <- as.POSIXct(Zeit) # hätte man früher gebraucht...
# 13.3 Plot und Achsenbeschriftung ---------------------------------------------
plot(Zeit, Daten[,3], type="l")
# Selber mit axis() Achsen beschriften ist problematisch.
# bei reinem Datum geht plotten und Beschriften noch recht einfach
y <- c(54, 63, 48)
x <- as.Date(c("2010-08-05", "2010-08-12", "2010-08-14")) # Default format
plot(x, y)
plot(x, y, xaxt="n")
axis(1, as.Date("7.8.2010","%d.%m.%Y"), expression( bold(hier) ) )
axis(1, x, format(x,"%d. %b\n%Y"))
plot(x, y, xaxt="n")
axis(1, pretty(x), format(pretty(x),"%d.%b %y"))
# bei Uhrzeit geht das besser mit axis.POSIXct
y2 <- c(54,63,48)
d2 <- c("05.08.2010 12:13:07", "05.08.2010 12:18:54", "05.08.2010 12:29:14")
x2 <- strptime(d2, "%d.%m.%Y %H:%M:%S")
plot(x2, y2) # X-Achsenbeschriftung in Stunden, Minuten und Sekunden
plot(x2, y2, xaxt="n")
k <- "2010-08-05 12:17:00 MEZ"
# automatisiert: (braucht CET: CentralEuropeanTime statt MEZ)
k2 <- as.POSIXct("05.08.2010 12:26:07", "%d.%m.%Y %H:%M:%S", tz="CET")
axis.POSIXct(1, at=k, format="%H:%M:%S")
axis.POSIXct(1, at=k2, format="%H:%M:%S")
axis.POSIXct(1, at="2010-08-05 12:22:00 MEZ", labels="hier")
# oder normal mit axis, wenn mit POSIXct-Objekten gearbeitet wird:
plot(x2, y2, xaxt="n")
axis(1, at=k2, labels=format(k2, "%H:%M:%S") )
abline(v=as.POSIXct("2010-08-05 12:22:00 MEZ"), lty=2) # http://vmware.ie/supp
# ort/developer/vc-sdk/visdk400pubs/ReferenceGuide/timezone.html
# Eigene Funktion mit labels an regelmäßigen Tagen des Monats:
source(url("http://dl.dropbox.com/u/4836866/BerryFunctions/R/monthLabs.R"))
monthLabs()
monthLabs(2014,2014, 3) # 3 days per month
monthLabs(2014,2014, npy=3) # 3 months per year, equally spaced
monthLabs(2014,2014, npy=4) # 4 months per year
plot(as.Date(c("2013-01-05", "2015-06-9")), 1:2)
plot(as.Date(c("2013-01-05", "2015-06-9")), 1:2, xaxt="n")
labs <- monthLabs(2013, 2015, npy=2) # npy: number of labels per year
axis(1, labs, format(labs, "%d.%m.%y") )
# 13.4 Zeitreihen aggregieren --------------------------------------------------
zr <- strptime(c("05.08.2010 12:15", "09.02.2013 14:45"), "%d.%m.%Y %H:%M")
zr <- seq(zr[1], zr[2], by=15*60) # Sequenz in 15-Minuten-Schritten
head(zr)
tail(zr)
fd <- cumsum(rnorm(length(zr))) # fake data aus Random Walk
plot(zr, fd, type="l")
plot(zr, fd, type="l", xaxt="n", main=c("täglicher Random Walk",""))
axis.POSIXct(1, as.Date(paste0(2010:2015,"-01-01")), labels=NA)
yearlabs <- as.Date(paste0(2010:2015,"-07-01"))
axis.POSIXct(1, at=yearlabs, labels=format(yearlabs, "%Y"), tick=F)
# Monatliche Mittelwerte ausrechnen (aggregieren):
monmean <- tapply(fd, format(zr, "%Y-%m"), mean) # apply the function mean on
monmean # the values fd that are first categorized into groups by month
# Mittelpunkt des Monats als Datumsobjekt:
monmid <- as.POSIXct(paste0(names(monmean), "-15"), format="%Y-%m-%d")
lines(monmid, monmean, col=2, lwd=2)
title(main=c(""," mit Monatsmittel"), col.main=2, cex.main=0.9)
# Summarize daily records by quarter of the year
set.seed(1492)
logins <- data.frame(date=as.Date("2014-01-01")+0:364,
failures=round(rlnorm(365)*20) )
# Extract month and quarter:
logins$month <- strftime(logins$date, "%m")
logins$qtr <- ceiling(as.numeric(logins$month)/3)
# Summarize by quarter:
tapply(X=logins$failures, INDEX=logins$qtr, FUN=sum)
tapply(X=logins$failures, INDEX=logins$qtr, FUN=median)
# unregelmäßige Messungen zu viertelstündlichen Werten aggregieren:
ourdata <- data.frame(Zeit=as.character(sort(Sys.time()+rnorm(15, sd=60*30))),
Wert=rexp(15))
str(ourdata)
times <- strptime(ourdata[,1], format="%Y-%m-%d %H:%M", tz="Asia/Manila")
# Achtung: sekunden werden hier einfach ignoriert, nicht gerundet:
times$min <- ceiling(times$min/15)*15
tapply(ourdata$Wert, as.POSIXct(times), FUN=sum, na.rm=T)
# für Niederschlagswerte sum, für Temperaturen mean, für Abflussmaxima max, etc.
# 13.5 Hydrologisches Jahr -----------------------------------------------------
# Effizientes Programmieren ist bei langen Zeitreihen wichtig !
# Nur einmal die StringParsedTimes erstellen, und dann davon jeweils Formate
# verwenden, zB beim Ermitteln der Jahresmaxima. (zr und fd aus 13.4)
# Anfang hydrologisches Jahr in Deutchlands Klima:
hyja <- which( format(zr,"%d.%m %H:%M") =="01.11 01:00")
Max <- sapply(1:(length(hyja)-1),
function(i) max(fd[ hyja[i]:(hyja[i+1]-1) ], na.rm=T))
# auch hier auf vorhandene zr zurückgreifen, und nur die gewollten auswählen
Erg <- data.frame(HydroJahr=format(zr[hyja[1:length(Max)]],"%Y"), Max)
Erg
# 13.6 nicht-equidistante Zeitreihen -------------------------------------------
# kumulierte Niederschlagswerte
set.seed(1234) # Startwert für den Zufallszahlengenerator (Reproduzierbarkeit)
raindata <- data.frame(time=strptime("2011-01-28 02:08", "%Y-%m-%d %H:%M") +
sort(sample(1:(400*24*3600), 150)),
rain=sample(0:10, 150, replace=TRUE))
head(raindata)
tail(raindata)
par(mfrow=c(3,1), mar=c(4,2.5,1,0.5))
plot(raindata[ 1: 50,], xlab="2011", type="h", col=4, las=1)
plot(raindata[ 51:100,], xlab="2011", type="h", col=4, las=1)
plot(raindata[101:150,], xlab="2011/2012", xaxt="n", type="h", col=4, las=1)
axis.POSIXct(side=1, at=monthLabs(2011,2012,npm=1), format="%b" )
# "zeitlich gebündelte" Regenfälle erkennbar
monthly_rain <- tapply(raindata$rain, format(raindata$time, "%Y-%m"), sum)
monthly_rain
monthly_rain <- monthly_rain[-1] # Januar entfernen, weil nur 1 Wert
par(mfrow=c(1,1), mar=c(5.5,4,2,0.5), lend="butt")
plot(monthly_rain, type="n", las=1, xlab="", ylab="Niederschlagssumme",
xaxt="n", main="Niederschlag", mgp=c(2.5, 1, 0),
ylim=c(0, max(monthly_rain)*1.04), yaxs="i")
abline(h=0:6*20, col=8) ; box()
points(monthly_rain, type="h", col=4, lwd=11) # kann man auch mit barplot machen
axis(1, 1:13, names(monthly_rain), las=2, tick=F, mgp=c(3,.5,0))
title(xlab="Monat", mgp=c(4,1,0)) # mgp[1] ist für Abstand des Achsentitels
# 13.7 Übungsaufgaben ----------------------------------------------------------
# 1. In welchen Jahren seit der Wende war der 15.10. ein Mittwoch?
# (Lässt sich mit 75 Zeichen in zwei Kommandozeilen herausfinden)
# 2. Ignorieren Sequenzen mit 20 minuten Abstand die Schaltsekunden?
# 3. Lese ein Datum ein, das wie folgt formatiert ist: "110809" (9.Aug.2011)
# Ausgehend davon, erstelle eine Sequenz von einer Stunde im Sekundenabstand
# 4. Schreibe eine Funktion die abhängig vom Geburtszeitpunkt dein Alter in
# Sekunden, Minuten, Stunden, Tage, Wochen, Monate, Jahre ausgibt.
# Optional: formatiere alle zahlen schick rechtsbündig untereinander.
# Stelle fest, wann du runde Jubiläen in den jeweiligen Einheiten feierst.
# Wann überholte dein Alter in Minuten die Weltbevölkerung?
# 5. Oftmals sind lange Zeitreihen nicht äquidistant:
# an welchen Stellen liegen Unterbrechungen vor und wie groß sind diese?
# Gib den Unterschied auch in Stunden an.
anfang <- strptime("05.08.2010 18:15:00", "%d.%m.%Y %H:%M:%S")
set.seed(42) ; zeitreihe1 <- anfang + c(0:5*15, sort(runif(25, 60, 600000)))
# Zeige die Zeitpunkte nach denen Messungen fehlen, mit Dauer der Lücke dahinter
zeitreihe2 <- ( anfang + 0:200*5*60 ) [-runif(30,1,201)]
# 13.8 Lösungen Übungsaufgaben -------------------------------------------------
# Aufg1: Wochentage
d <- as.Date(paste(1990:2014, "10-15", sep="-"))
d[ format(d, "%a") == "Mi"] # 1997, 2003, 2008, 2014
# Aufg2: Sequenzen, Schaltsekunden
tail(.leap.seconds)
kk <- strptime(c("2012-06-30 23:20", "2012-07-01 03:20"), "%Y-%m-%d %H:%M")
seq(from=kk[1], to=kk[2], by=20*60 ) # ignoriert sie (Minuten sind noch Rund)
kk[1] + 1:12*20*60 # ditto
data.frame(k=strptime("2012-07-01 01:59:50", "%Y-%m-%d %H:%M:%S") + 1:20)
# Mein System ignoriert die, wie erklärt unter Details der Hilfe zu
?.leap.seconds
# Aufg3: Formatkonvertierungen
myDate <- as.POSIXct("110809", format = "%y%m%d")
myDate
head(myDate + 0:3600)
# Aufg4: Alter
alter <- function(geburt, abfrage=Sys.time(), units=c("s","m","h","d","w"))
{
weite <- function(k)
{ifelse(k>=1e9, 14, ifelse(k>=1e6, 15, ifelse(k>=1e3, 16, 17)))}
oldopt <- options(scipen=10)
for (i in 1:length(units))
{
k <- difftime(abfrage, strptime(geburt, "%Y-%m-%d %H:%M"), units=units[i])
cat(format(k, big.mark=" ", nsmall=2, digits=2, width=weite(k)), "\n")
}
k <- difftime(abfrage, strptime(geburt, "%Y-%m-%d %H:%M"), units="d")
cat(format(as.numeric(k2 <- k/(365.2425/12)), big.mark=" ", nsmall=2,
digits=2, width=weite(k2)), "months \n")
k <- difftime(abfrage, strptime(geburt, "%Y-%m-%d %H:%M"), units="d")
cat(format(as.numeric(k2 <- k/365.2425), big.mark=" ", nsmall=2,
digits=2, width=weite(k2)), "years \n")
options(oldopt)
}
alter(geburt="2005-08-11 07:05")
pretty(736435817.66)
1e9/3600/24/365.2425 # 1 billion (dt: milliarde) sekunden in 31.7 Jahren
13e6 /60/24/365.2425 # 13 Mio Minuten bei 24.7 Jahren
# Aufg5: Lücken
diff(zeitreihe1)
difftime(zeitreihe1[-1], zeitreihe1[-length(zeitreihe1)], units="h")
data.frame(Luecke_Nach=zeitreihe2[-length(zeitreihe2)],
Lueckenlaenge=diff(zeitreihe2))[diff(zeitreihe2)!=5,]
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 14: Regression (Quadratmittelapproximation) ########################
# 14.1 lineare Regression (Anpassung einer Gerade an Wertepaare)
# 14.2 Geradengleichung ermitteln und verwenden
# 14.3 logarithmische Regression
# 14.4 multivariate Regression
# 14.5 cubic regression
# 14.6 Funktion
# 14.1 lineare Regression (Anpassung einer Gerade an Wertepaare) ---------------
x <- c(5,8,10,9,13,6,2) ; y <- c(567,543,587,601,596,533,512)
plot(x,y, las=1, pch=16)
modellgerade <- lm(y~x) # y abh. von x. xy-Reihenfolge anders als sonst meist!
modellgerade
?formula
# lm(y~x): lineares Modell mit y abhängig von x
# y = a*x + b, a = Steigung (slope), b = Konstante, Achsenabschnitt (intercept)
# y: abhängige (dependent) Variable , criterion / response
# x: erklärende unabhängige (independent) Variable , predictor
summary(modellgerade)
# Zur Interpretation siehe zB:
browseURL("http://blog.yhathq.com/posts/r-lm-summary.html")
# http://stats.stackexchange.com/questions/5135/interpretation-of-rs-lm-output
# Siehe 9.1 Hypothesentests Allgemein
# Nullhypothese beim Linearen Modell: Zusammenhang zufällig ("Rauschen")
# Also ist p die Wahrscheinlichkeit, dass der Parameter (die Variable)
# - das Modell nicht verbessert (multivariate Regression)
# - und nicht zur Erklärung von y relevant ist.
# Wenn p < alpha (oft = 0.05): Variable drin lassen. Es gibt einen Grund für die
# Korrelation, auch wenn R^2 gering ist ("Signal")
# Ziel LinMod: R2 hoch, Anzahl Parameter klein (Parsimonie prinzip, aka
# Ockhams Razor: so viele Variablen wie nötig, so wenige wie möglich).
str(modellgerade) # sämtliche Info, die im Objekt enthalten ist!
abline(modellgerade, lwd=3, col=4)
cor(x, y)
# A high correlation does NOT imply correct causation!
# the strong negative correlation childbirth ~ stork abundancy does not mean
# getting more storks to breed here will help us fight demographic change 😉
cov(x,y)/(sd(x)*sd(y)) # Kovarianz durch Standardabweichungen
cor(x, y)^2 # = R^2
summary(modellgerade)$r.square # 0.6685065
p <- predict(modellgerade)
p # Achtung, für jeweilige x-Punkte, aber in aufsteigender Reihenfolge!
points(x, predict(modellgerade), col=2, pch=16)
# Beim lm wird die Summe der Distanzen zum Quadrat minimiert:
segments(x, y, x, p)
# oder umständlich: for (i in 1:7) lines( x=c(x[i], x[i]), y=c(y[i], p[i]) )
text(x, (y+p)/2, round(y-p,1), adj=-0.2)
# Die Linien stellen die Residuen (Restfehler) dar
residuals(modellgerade)
sum((y-p)^2) # 2328.884 SSE (Sum of Squared Errors) oder RSS: (Residual SS)
sum((p-mean(y))^2)/sum((y-mean(y))^2) # Variation Regresswerte / Variation Y
1 - sum((y-p)^2)/sum((y-mean(y))^2)
# Unsicherheit:
plot(x,y, las=1, pch=16) ; abline(modellgerade, lwd=3, col=4)
x2 <- seq(1, 14, length=30)
p2 <- predict(modellgerade, newdata=data.frame(x=x2), interval="confidence" )
polygon(x=c(x2, rev(x2)), y=c(p2[,2], rev(p2[,3])), col=rgb(0.3,0.5,0, 0.5))
# ANOVA # mehr zu ANOVA im Kapitel 9.5
anova(modellgerade)
summary(anova(modellgerade))
# 14.2 Geradengleichung ermitteln und verwenden --------------------------------
# bei jedem ersten Durchlauf an jedem Rechner gleiche "Zufallszahlen"
set.seed(12)
x <- c(runif(100,0,3), runif(200, 3, 25)) # random aus uniformer Verteilung
y <- 12.367*log10(x)+7.603+rnorm(300) # random aus Normalverteilung
plot(x, y, las=1, pch=20, main="Verdeutlichung lin.regr. vs. log.regr.")
# Koeffizienten des linearen Regressionsmodells bestimmen
linmod <- lm( y ~ x )
linmod
str(linmod) # Liste mit 12 Elementen
abline(linmod, col=2, lwd=3)
summary(linmod)
r <- round( summary(linmod)$r.squared ,2)
m <- round( coef(linmod)[2] ,2 ) # runde die zweite Koefiziente auf 2 Stellen
b <- round( coef(linmod)[1] ,2 )
Txt <- paste("y =", m, "* x +", b, "\nR\U00B2 =", r) #\n macht Zeilenumbruch
Txt
text(12,10, Txt, adj=0, col=2)
# Predict
predict(linmod) # Standardmäßig für Werte, an denen gefitted wurde
predict(linmod, newdata=data.frame(x=seq(0,20,by=4)))
# colnames() des Dataframes muss genau den Namen enthalten, der auch für
# lm verwendet wurde...
# 14.3 logarithmische Regression -----------------------------------------------
# Koeffizienten des logarithmischen Regressionsmodells bestimmen
# Achtung: log berechnet ln (Log zur Basis e), log10 den dekadischen Logarithmus
# Negative Werte ersetzen, sonst können wir keinen Logarithmus berechnen
# x <- replace(x, x <= 0, 0.0001) # in den Beispieldaten nicht nötig
# ist eine Veränderung der Daten und daher nur begründet anzuwenden!
set.seed(12)
x <- c(runif(100,0,3), runif(200, 3, 25)) # random aus uniformer Verteilung
y <- 12.367*log10(x)+7.603+rnorm(300) # random aus Normalverteilung
plot(x,y)
logmod <- lm( y ~ log10(x) ) ; summary(logmod)
r2 <- round( summary(logmod)$r.squared ,2)
m2 <- round( coef(logmod)[2] ,2 )
b2 <- round( coef(logmod)[1] ,2 )
xvekt <- seq(0,35,0.01)
lines(xvekt, m2*log10(xvekt)+b2, col=4, lwd=3)
text(12,5, paste("y=", m2, "* log(x) +", b2, "\nR\U00B2 =", r2), adj=0, col=4)
# 14.4 multivariate Regression -------------------------------------------------
# Regression mit mehreren erklärenden Variablen A,B,C
# Syntax (Pseudo-code): lm(y ~ A + B + C)
#mit Transformationen müsste sowas gehen: lm( y ~ A + I(log(B)) + C)
browseURL("http://dl.dropboxusercontent.com/u/4836866/Sonstiges/lakes.txt")
lakes <- read.table("Daten/14_lakes.txt", header=T )
pairs(lakes)
textcor <- function(x,y)
legend("center", legend=round(cor(x,y),2), cex=2, bty="n", adj=0.4)
pairs(lakes, lower.panel=textcor)
summary( lm(chlorophyll ~ P, data=lakes) )$r.squared # 0.77
summary( lm(chlorophyll ~ P + N + temp, data=lakes) )$r.squared # 0.79
# R^2 etwas höher, aber:
summary( lm(chlorophyll ~ P + N + temp, data=lakes) )
# Von 3 Variablen nur einer, der signifikant einen Anteil der Varianz erklärt.
# Nach Ockhams Razor (Sparsamkeitsprinzip) werden die anderen beiden weggelassen
# "Steht man vor der Wahl mehrerer möglicher Erklärungen für dasselbe Phänomen,
# soll man diejenige bevorzugen, die mit der geringsten Anzahl an Hypothesen
# auskommt und somit die 'einfachste' Theorie darstellt." [Wikipedia]
# Also sukzessive die Variable mit dem höchsten P-value rauslassen,
# dabei die Änderung des R2 etwas mit im Blick behalten.
# Siehe Hinweis Zeile 30.
cor(lakes$N, lakes$P) # erklärende Variablen sind untereinander korreliert.
pnt_mod <- lm(chlorophyll ~ P + N + temp, data=lakes)
pnt_mod
# chlorophyll = 0.34*P + 1.14*N + 2.97*t - 33.7
K <- coef(pnt_mod)
K
K[2] # ist ein Vector, auch wenn es wie ein data.frame aussieht
class(K) # numeric
# Schätze die Chlorophyllkonzentration für eine Messung mit P=100, N=15, t=9
K[1] + 100*K[2] + 15*K[3] + 9*K[4] # 44
# zeigt als Namen den ersten Namen von k an, daher
as.vector( K[1] + 100*K[2] + 15*K[3] + 9*K[4] )
View(data.frame(lakes, Chl_ESTIMATE=0.3409*P + 1.1403*N + 2.9687*temp -33.7333))
# Random Forest - wichtigkeit einzelner Prädiktoren
set.seed(131)
ozone.rf <- randomForest(Ozone ~ ., data=airquality, mtry=3,
importance=TRUE, na.action=na.omit)
print(ozone.rf)
## Show "importance" of variables: higher value means more important:
round(importance(ozone.rf), 2)
# Correlation matrix:
round(abs(cor(airquality, use="complete")), 2)
# compare cor and imp
correlation <- c(0.35,0.61,0.70,0.14,0.01)
importance <- c(11.97,23.22,45.22,2.10,2.83)
berryFunctions::linReg(correlation, importance)
# 14.5 cubic regression --------------------------------------------------------
a <- c(-5:7, 10) ; b <- 5*a^3 + 2*a^2 - 7*a + 3 + rnorm(14, 0, 100)
data.frame(a,b) ; plot(a,b)
modell <- lm(b ~ poly(a,3, raw=T)) ; modell # gleich mehr zu raw
lines(a, predict(modell))
xdraw <- seq(3, 11, length=1e4)
lines(xdraw, predict(modell, data.frame(a=xdraw)), col=2)
coef(modell)
dummy.coef(modell)
poly(a,3, raw=T)
poly(a,3) # Beim Default-Wert von Raw (FALSE) wird normiert!
i <- seq (-5,10,0.1) ; j <- 3*i^2 - 5*i + 4
lm(j ~ poly(i,2)) # kommt was falsches bei rum
lm(j ~ poly(i,2,raw=T)) # jetzt richtig
# eine andere Möglichkeit, "+" auch für multivariate Regression:
lm(b ~ I(a^3) + I(a^2) + a)
# 14.6 Funktion ----------------------------------------------------------------
# Mit Funktion Werte in die Graphik schreiben (toll für Rprofile oder Paket)
linreg <- function(x,y, pos="top", new=TRUE, pch=20, col=2, lwd=3, inset=0, ...)
{ if (new) plot(x, y, las=1, pch=pch, ...)
mod <- lm( y ~ x )
abline(mod, col=col, lwd=lwd, ...)
r <- round( summary(mod)$r.squared ,2)
m <- round( coef(mod)[2] ,2 )
b <- round( coef(mod)[1] ,2 )
Txt <- paste0("y = ", m ," * x ", if(b>0)" + ", b, "\nR\U00B2 = ", r)
legend(pos, Txt, bty="n", text.col=col, inset=inset) }
# linreg(x, y, new=FALSE) fügt die Werte einem bestehenden (!) Plot hinzu
# Wäre leicht um ein digits-Argument für die Nachkommastellen zu erweitern.
a <- 1:30
b <- a/2.345+rnorm(30,0,3)
linreg(a,b, ylab="Hallo", pch=1, col=3)
plot(a,b, xlim=c(-5,45))
linreg(a, b, pos="bottomright", new=F, inset=.1) # inset für Abstand vom Rand
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 15: Extremwertstatistik für Hochwassermaxima #######################
# 15.1 Datensatz
# 15.2 Jährlichkeit berechnen
# Exkurs: Verschiedene Methoden zur Bestimmung der PP
# 15.3 Verteilung anpassen
# 15.4 Zeichnen der Abflüsse über ihre Jährlichkeiten
# 15.5 Abflüsse für bestimmte Jährlichkeiten
# 15.6 Funktion mit mehreren Verteilungen : L-Momente
# 15.7 Funktion distLextreme
# 15.8 andere Datensätze
# Ich bin kein Experte für Verteilungsfunktionen, sämtliche Kommentare sind also
# bitte mit Vorsicht zu genießen...
# In meinem Paket sind einige hier benötigte Funktionen: movAv
install.packages("berryFunctions") ; library(berryFunctions)
# 15.1 Datensatz ---------------------------------------------------------------
# Jahresmaxima des Abflusses 1976-2010 an einem Bach in Österreich.
JM <- c(61.5, 77.0, 37.0, 69.3, 75.6, 74.9, 43.7, 50.8, 55.6, 84.1, 43.6, 81.9,
60.1, 72.4, 61.6, 94.8, 82.6, 57.2, 63.1, 73.8, 51.3, 93.6, 56.9, 52.1, 40.4,
48.9, 113.6, 35.4, 40.1, 89.6, 47.8, 57.6, 38.9, 69.7, 110.8)
# Dies müssen unbedingt unabhängig voneinander sein! Es darf keine
# Jahrestrennung an einem Extremereignis erfolgen (das dann in beiden Jahren
# als Jahresmaxima auftaucht).
# Zeitreihe anschauen
plot(1976:2010, JM, type="l", las=1, ylab="Jahresmaxima Abfluss [m\U00B3/s]")
# Moving Average (in berryFunctions):
lines(1976:2010, movAv(JM), col=2, lwd=3) # kann ein Indikator für Trends oder
# Zyklizitäten sein, aber: Fensterbreite = "researcher degree of freedom"
# Wer die Theorie hinter der Extremwertstatistik schon kennt, kann gleich zu
# 15.7 springen. Ansonsten soll das Folgende einen Eindruck vermitteln...
# Wie immer ist Feedback sehr willkommen!
# 15.2 Jährlichkeit berechnen --------------------------------------------------
quantile(JM) # 75% der Werte sind kleiner als 76.3 m3/s (oberes Quartil)
# Dieser Abfluss wird in 25 von 100 Jahren überschritten, kann also statistisch
# gesehen alle 4 Jahre auftreten (bei gleichbleibendem Abflussregime).
# Bestimmung der Plotting Positions (PP). Früher verwendet zur grafischen
# Anpassung der Extremwertverteilung als Gerade auf Gumbel-skalierter Achse
# Plotting Positions nach Weibull:
m <- sort(JM) # aufsteigend sortierte Maxima
Rang <- 1:length(m)
Pu <- Rang / (length(m)+1) # empirische Unterschreitungswahrsch. = Perzentil
par(mfrow=c(1,2))
plot(m,Pu, xlab="Abflusswert", ylab="Unterschreitungshäufigkeit", las=1)
# Vergleichend die empirische kumulierte Wahrscheinlichkeitsdichte:
plot(ecdf(JM)) # inhaltlich fast das gleiche, jedoch von 0 bis 1, also
# keine Möglichkeit für noch höhere Abflüsse, die es aber gibt.
# wieder zum ersten Plot:
plot(m,Pu, xlab="Abflusswert", ylab="Unterschreitungshäufigkeit", las=1)
# Wiederkehrintervall = Jährlichkeit = 1/Pü = 1/(1-Pu):
RP <- 1/(1-Pu) # RP: Return Period
plot(m,RP, xlab="Abflusswert", ylab="Wiederkehrintervall [Jahre]", las=1)
# Der höchste Abflusswert hat eine Jährlichkeit von 36.
# Der zweithöchste Wert wird in den 35 Jahren zweimal ge- oder übertroffen,
# tritt statistisch gesehen 2x in 35 Jahren, also alle 18 Jahre auf.
# Überschreitung des dritten Wertes passiert alle 12 Jahre. etc.
# Diese Reziproke (1/x) Beziehung bewirkt die unintuitive Abstände zwischen den
# Return Periods.
# Exkurs: Verschiedene Methoden zur Bestimmung der PP --------------------------
browseURL(paste0("http://www.uni-potsdam.de/db/fsr-ggr/data/archiv/klausuren/",
"Hydrologie/Extremwertstatistik.pdf"))
browseURL("www.tandfonline.com/doi/abs/10.1080/03610920701653094#.U0UC7FczZqM")
browseURL("http://journals.ametsoc.org/doi/abs/10.1175/JAM2349.1")
browseURL(paste0("http://www.whycos.org/fck_editor/upload/File/Pacific-HYCOS/",
"Surface_Waters/Estimating_Flood_Frequency.pdf"))
# Plotting Positions nach Hazen (1914):
# Formeln aus dem dritten Link
Pu.hazen <- (Rang-0.50)/ length(m)
# Plotting Positions nach Gringorten (1963) (wie in lmom:::evplot.default):
Pu.gring <- (Rang-0.44)/(length(m)+0.12)
# Plotting Positions nach Beard (1943):
Pu.beard <- (Rang-0.31)/(length(m)+0.38)
# Plotting Positions nach Weibull (1939):
Pu.weibu <- Rang /(length(m)+1)
RP.hazen <- 1/(1-Pu.hazen)
RP.gring <- 1/(1-Pu.gring)
RP.beard <- 1/(1-Pu.beard)
RP.weibu <- 1/(1-Pu.weibu)
legtitles <- c(
"Hazen (1914) : (Rang-0.50)/ n",
"Gringorten (1963): (Rang-0.44)/(n+0.12)",
"Beard (1943) : (Rang-0.31)/(n+0.38)",
"Weibull (1939) : Rang /(n+1)")
par(mfrow=c(1,1))
plot( m, RP.hazen, las=1, type="n", main="Methoden für Plotting Positions",
ylab="Wiederkehrintervall = RP = 1/(1-Pu)", xlab="Abfluss Jahresmaxima")
lines(m, RP.hazen, col=1, type="o", lty=1, pch=1)
lines(m, RP.gring, col=2, type="o", lty=2, pch=2)
lines(m, RP.beard, col=3, type="o", lty=3, pch=3)
lines(m, RP.weibu, col=4, type="o", lty=4, pch=4)
par(family="mono")
legend("topleft", legtitles, col=1:4, lty=1:4, pch=1:4, title="P_unterschr. =")
par(family="")
# macht bei den hohen Werten einen großen Unterschied!
# Entsprechend dann auch in der Berechnung der n-jährlichen Abflüsse!
# Hazen, Gringorton und Beard überschätzen Wiederkehrinterval und somit
# unterschätzen sie das Risiko! (verglichen mit Weibull)
# I tend to like the Weibull method for Plotting Positions more.
# With 35 years of data, Gringorton estimates a return period of 60 years.
# But if we were to take n -> Inf samples of 35, in the long run, the highest of
# these values should have a recurrence interval of ~35 years.
# 15.3 Verteilung anpassen -----------------------------------------------------
# Abflussvektor für die Graphik Qrg <- extendrange(m, f=0.1)
Q <- seq(par("usr")[1], par("usr")[2], len=500)
# Parameter der Gammaverteilung aus den Momenten der Daten bestimmen
shape <- mean(m)^2/var(m)
rate <- mean(m)/var(m)
# Die plotting positions sind für das fitten der Verteilungen irrelevant.
# Verteilung über den Bereich Q (Graphikbereich an Abflüssen)
hist(m, freq=FALSE, col=3)
lines(Q, dgamma(x=Q, shape=shape, rate=rate), lwd=2)
# Trifft halbwegs zu... Allerdings ist das Tailing der Werte ausgeprägter!
# Neben der Gammaverteilung gibt es noch eine Menge andere, zB GEV, GPD, Gumbel,
# normal, exponential, logistic, usw. Siehe Abschnitt 15.6 bis 15.8
# 15.4 Zeichnen der Abflüsse über ihre Jährlichkeiten --------------------------
plot(RP, m, ylab="Discharge HQ [m\U00B3/s]", xlab="Return Period RP [a]",
las=1, main="yearly Discharge Extrema / Return Period", pch=20)
# Return Period für Abflussvektor nach der so parametrisierten Gammaverteilung:
RPgamma <- 1/( 1- pgamma( Q, shape=shape, rate=rate) )
# Einzeichnen der Extremwertverteilung
lines(RPgamma, Q)
legend("right", legend=c("Jahresmaxima 1976-2010", "Gammaverteilung"),
lty=c(NA,1), pch=c(20,NA))
# Schreiben der Parameter
text(15,60, paste("mean =",round(mean(m),3)," sd =",round(sd(m),3)), adj=0 )
text(15,53, bquote("shape = "* mu^2/sigma*" = "*.(round(shape, 2))), adj=0 )
text(15,46, bquote("rate = " * mu/sigma * " = "*.(round(rate, 4))), adj=0 )
text(15,39, bquote("scale = "* 1/rate * " = "*.(round(1/rate,3))), adj=0 )
# Auch hier wird das Unterschätzen hoher Abflüsse deutlich!
# Die Achsen sind bewusst linear skaliert - das kann man sich visuell viel
# besser vorstellen.
# 15.5 Abflüsse für bestimmte Jährlichkeiten -----------------------------------
ReturnYears <- c(5,10,20,50)
qgamma(1-(1/ReturnYears), shape=shape, rate=rate)
# 15.6 L-Momente ---------------------------------------------------------------
if(!require(lmom)){install.packages("lmom");require(lmom)} # samlmu, pelgev
if(!require(evd)){install.packages("evd");require(evd)} # dgev
# Lineare Momente der Stichprobe
momente <- samlmu(JM, sort.data=TRUE, nmom=3)
momente
browseURL("http://en.wikipedia.org/wiki/L-moment")
# Parameter der General Extreme Value Distribution anhand der L-Momente schätzen
param <- pelgev(samlmu(JM)) # PEL: Parameter Estimation L-moments
param # Die Hilfe zu pelgev hat eine Übersicht der implementierten Verteilungen
# Angepasste General Extreme Value Distribution einzeichnen:
hist(JM, col=4, freq=FALSE, breaks=10)
lines(density(JM), lwd=1, lty=2)
x <- seq(20,130, len=100)
lines(x, dgev(x=x,loc=param["xi"],scale=param["alpha"],shape=param["k"]), lwd=3)
# Verteilung ist mäßig geeignet, auch die anderen sind auszuprobieren
# Quantil der so parametrisierten GEV Verteilung (analog zu qnorm):
quagev(f=0.99, param) # 99% der Werte dieser GEV Verteilung liegen unter 127.5
# Extreme Value plot mit nicht-linearer Gumbel-Achse:
evplot(JM, rp.axis = TRUE, ylab= "Discharge [m\U00B3/s]")
evdistq(qfunc=quagev, param, col ="black")
# Wert der Verteilung für ReturnPeriod 50 Jahre finden:
quagev(f=1-1/50, para=pelgev(momente)) # aus exp(-exp(log(-log(1-1/50))))
quagam(f=1-1/50, para=pelgam(momente))
# Unterschiedliche Datensätze brauchen unterschiedliche Verteilungen.
# Die GEV passt nicht immer am besten!
# Unsicherheit der Parameter bestimmen via Bootstrapping, hier mit Leave5Out:
q99gev <- sapply(X=5:length(JM), FUN=function(i){
param <- pelgev( samlmu(JM[-((i-4):i)], sort.data = TRUE) )
evdistq(qfunc=quagev, param, col ="black") # Verteilung mit einzeichnen
quagev(0.99, param) }) # 99% Quantil bestimmen und rausgeben
# Da sind schon Abweichungen in der Graphik erkennbar!
# Wenn man für die Berechnung des Hochwassers einer gewissen Jährlichkeit einen
# Unsicherheitsbereich angeben will, ist das eine der Methoden.
as.numeric(q99gev)
hist(q99gev, breaks=15, col=3)
# Mögliche andere Verteilungen - Generalized Pareto Distribution (GPD):
plot(ecdf(JM))
lines(x, cdfgpa(x, pelgpa(samlmu(JM))), col=2) # Naja, geht so
quagpa(0.99, pelgpa(samlmu(JM))) #99% der Werte dieser GPD sind < 113.6
# Das ist schon ein wenig anders als 127.5!
# To-Do: wie relevant werden Verteilungen durch L-momente anders bestimmt
# als durch klassische statistische Momente?
# 15.7. Funktion distLextreme --------------------------------------------------
install.packages("devtools")
devtools::install_github("brry/extremeStat")
library(extremeStat)
library(lmomco)
distLextreme(JM)
# 50-jähriges Ereignis anhand verschiedener Verteilungen abschätzen:
Q50 <- distLextreme(JM, RPs=50, plot=FALSE)$returnlev
sort(Q50[,1], decr=TRUE)
# das ist eine ganze Bandbreite an Möglichkeiten...
# 107 oder 132 m3/s kann schon einen relevanten Unterschied machen!
# Farben der Linien und Auswahl zu plottender Verteilungen:
dle <- distLextreme(JM)
distLextreme(dlf=dle, coldist=heat.colors(10), nbest=15)
distLextreme(dlf=dle, coldist=heat.colors(13), selection=c(4,9,1))
distLextreme(dlf=dle, coldist=1:3, selection=c(4,9,1))
distLextreme(dlf=dle, legarg=list(cex=0.5, x="bottom", box.col="red", col=3))
# Nach Anpassungsgüte gewichtetes Mittel der verschiedenen Verteilungen:
?distLextreme # Section Examples
# 15.8 andere Datensätze -------------------------------------------------------
# MackenzieRiver
MR <- c(26600,30300,34000,32000,29200,28300,28600,26400,28300,28800,
29000,22100,32900,31800,21600,32100,27000,24800,28000,35000,32000,25000,15800,
28800,29900,28000,25600,19700,25700,29500,26800,30000,29500)
?Nile # 100 years!
evplot(MR, ylab= "Discharge [m\U00B3/s]")
evdistq(qfunc=quagev, para=pelgev(samlmu(MR)), col ="black")
extremeStatLmom(MR)
extremeStatLmom(MR, returnParam=T)
extremeStatLmom(Nile)
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### Kapitel 16 - Geostatistik ##################################################
# 16.0 im Kapitel benoetigte Pakete
# 16.1 ESRI-Shapefiles laden und bearbeiten
# 16.2 Kriging (also die eigentliche Geostatistik)
# 16.3 Karten im Hintergrund
# 16.0 im Kapitel benoetigte Pakete ---------------------------------------------
# Runterladen der Pakete nur einmal, danach auskommentieren
install.packages("geoR")
install.packages("rgeos")
install.packages("maptools")
install.packages("tmap")
install.packages("berryFunctions")
install.packages("rgl")
# Laden der Pakete nach jedem R-Start (oder in Rprofile, siehe Kap 10)
library(maptools) # fuer ReadShapeSpatial
library(geoR) # as.geodata, variog, variofit, krige.control, krige.conv,
library(sp) # geoR bereits installiert legend.krige
library(rgeos) # gCentroid, gArea, gUnion
library(tmap) # qtm
library("berryFunctions") # colPoints, classify
help(package="geoR") # und
library(help="geoR")
data(Europe)
data(World)
# Karten mit richtiger Projektion:
qtm(World, fill = "economy", text="iso_a3", text.size = "AREA",
fill.palette="-Blues", theme = "World", fill.title="Economy")
qtm(Europe, fill="gdp_cap_est", text="iso_a3", text.size="pop_est",
fill.title="GDP per capita", fill.textNA="Non-European countries")
class(Europe) # SpatialPolygonsDataFrame entspricht ESRI shapefile
# 16.1 ESRI-Shapefiles laden und bearbeiten ------------------------------------
browseURL("srtm.csi.cgiar.org") # globales DGM
browseURL("http://diva-gis.org/Data") # Grenzen aller Staaten
browseURL("http://spatialreference.org/") # proj4strings
# IndiaDaten fuer das Beispiel unten gibt's hier:
# http://globalresourcesuae.com/www.maptell.com/index3e60.html?Itemid=159
india <- readShapeSpatial("Daten/Shapefiles/india_state.shp")
# der "Fehler in getinfo.shape(fn) : Error opening SHP file"
# bedeutet oft, dass das Arbeitsverzeichnis nicht richtig gesetzt ist -> getwd()
# oder der Dateiname falsch geschrieben ist. -> dir()
class(india) # SpatialPolygonsDataFrame
# Some states are made up of more than one polygon.
# State names are in the NAME column of the dataframe
View(india)
head(india@data)
str(india[1,])
plot(india, col = india$NAME)
plot(india[india$NAME == "Tamil Nadu",], col = 1:8)
attr(methods(plot),"info") # sp has a plot method
?"SpatialPolygons-class"
methods(class="SpatialPolygonsDataFrame")
plot(india, col=seqPal(169)); axis(1);axis(2,las=1);
colPointsLegend(1:169, nbin=169, col=seqPal(169))
plot(india[72,]); axis(1);axis(2,las=1)
head(india[72,]@polygons[[1]]@Polygons[[1]]@coords)
tail(india[72,]@polygons[[1]]@Polygons[[1]]@coords)
# Extract centroid for a State (this works):
gCentroid(india[india$NAME == "Tamil Nadu",])
# Error arises for this State:
gCentroid(india[india$NAME == "Gujarat",])
# orphaned hole, cannot find containing polygon for hole at index 2
plot(india[india$NAME == "Gujarat",], col = 1:8)
# 16.2 Kriging (also die eigentliche Geostatistik) -----------------------------
# 16.2.1 Daten Beispiel --------------------------------------------------------
x <- c(1,1,2,2,3,3,3,4,4,5,6,6,6)
y <- c(4,7,3,6,2,4,6,2,6,5,1,5,7)
z <- c(5,9,2,6,3,5,9,4,8,8,3,6,7)
plot(x,y, pch="+", cex=z/4)
colPoints(x,y,z, col=addAlpha(rainbow2(50)), cex=3)
GEODATEN <- as.geodata(cbind(x,y,z))
plot(GEODATEN)
# 16.2.2 Variogramm ------------------------------------------------------------
EXP_VARIOGRAMM <- variog(GEODATEN)
plot(EXP_VARIOGRAMM)
FIT_VARIOGRAMM <- variofit(EXP_VARIOGRAMM)
lines(FIT_VARIOGRAMM)
# 16.2.3 Kriging (Interpolation unter Beruecksichtigung der Semivarianz) --------
aufloesng <- 0.1
gitter <- expand.grid(seq(min(x),max(x),aufloesng),seq(min(y),max(y),aufloesng))
krico <- krige.control(type.krige="OK", obj.model=FIT_VARIOGRAMM)
krobj <- krige.conv(GEODATEN, loc=gitter, krige=krico) # KRigingObjekt
# 16.2.4 Farbliche Darstellung -------------------------------------------------
colPoints(x,y,z, add=F)
image(krobj, col=rainbow2(100))
legend.krige(col=rainbow2(100), x.leg=c(6.2,6.7), y.leg=c(2,6), vert=T, off=-0.5,
values=krobj$predict)
contour(krobj, add=T)
colPoints(x,y,z, col=rainbow2(100)) # verstreute Punkte nach Hoehe gefaerbt
points(x,y)
# 16.2.5 3D Darstellung --------------------------------------------------------
Farben <- rainbow2(100)[classify(krobj$predict)$index]
library(rgl)
open3d() # Fenster in Sichtweite ziehen : Maus ziehen und scrollen
plot3d(gitter[,1], gitter[,2], krobj$predict, type = "p",
xlab = "X", ylab = "Y", zlab = "Z", col=Farben)
plot3d(gitter[,1], gitter[,2], krobj$predict, type = "s", size=1,
xlab = "X", ylab = "Y", zlab = "Z", col=Farben)
# Letzteres mit allen Schatten zu rendern ist ein Heidenaufwand, hab also Geduld
# 16.2.6 3D Mesh (durchgaengige Oberflaeche) -------------------------------------
# aus ?mesh3d:
vertices <- c( # jeweils x,y,z,winkel; siehe ?rgl::matrices
-1, -1, 0, 1,
1, -1, 3, 1,
1, 1, 0, 1,
-1, 1, 0, 1 )
indices <- c( 1, 2, 3, 4)
open3d()
shade3d( qmesh3d(vertices,indices), col=2:5 )
open3d()
wire3d( qmesh3d(vertices,indices) , col=2:5, lwd=3)
# jetzt fuer Beispieldaten:
dummy <- data.frame(x=gitter[,1], y=gitter[,2], z=krobj$predict, w=1)
vertices <- unlist(t(as.matrix(dummy))) ; rm(dummy)
round(vertices[1:30],2)
open3d()
shade3d( qmesh3d(vertices, rep(1:4, length(vertices)/4)) )
# nimmt nur die ersten 4 Punkte...
# 16.2.7 3D perspective plot ---------------------------------------------------
nx <- length(seq(min(x),max(x),aufloesng))
ny <- length(seq(min(y),max(y),aufloesng))
krmat <- matrix(krobj$predict, nrow=nx)
Farben <- rainbow2(100)[classify(krmat)$index]
persp(1:nx, 1:ny, krmat, col=Farben) # farben stimmen nicht
plot(1,1, xlim=c(1,51), ylim=c(1,61), type="n") # hier stimmts aber:
for (x in 1:51) for (y in 1:61) points(x,y, col=Farben[x,y], pch=16)
# mehr zu 3D-Darstelluungen in Kapitel 4.14
# 16.3 Karten im Hintergrund ---------------------------------------------------
browseURL("www.java.com") # install JAVA in same bit-version as R (eg 64bit):
browseURL("http://www.java.com/de/download/manual.jsp")
if(!requireNamespace("OpenStreetMap")) install.packages("OpenStreetMap")
library("OpenStreetMap")
d <- read.table(header=TRUE, text=
"lon lat elevation
13.01480 52.39975 66.93333
13.01233 52.40032 67.00000
13.00975 52.40110 64.76667")
plot(lat~lon, d)
bbox <- c(extendrange(d$lon), extendrange(d$lat, f=0.5))
map <- openmap(bbox[c(4,1)], bbox[c(3,2)], type="osm")
par(mar=c(0,1,3,0))
plot(map, removeMargin=FALSE)
points(projectMercator(d$lat,d$lon), pch=3, lwd=3)
lines(projectMercator(d$lat,d$lon) )
berryFunctions::scaleBarOSM()
library("mapmisc")
proj4string(map)
mapmisc::scaleBar(map) # doesn't work...
#SDMTools::Scalebar(1448400, 6872900, distance=0.5, unit="km")
map_aer <- openmap(bbox[c(4,1)], bbox[c(3,2)], type="bing")
plot(map_aer, removeMargin=FALSE)
points(projectMercator(d$lat,d$lon), pch=3, lwd=3, col=2, cex=1.5)
# I find type="maptoolkit-topo" to work very fine as well. For all types, see
browseURL("http://www.r-bloggers.com/the-openstreetmap-package-opens-up")
# Project map to UTM Coordinates
map_utm <- openproj(map, "+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84")
plot(map_utm)
axis(1, line=-2); axis(2, line=-2)
# Reading GPX files, eg from apps like sportractive or OSMtracker:
#install.packages("plotKML")
d <- plotKML::readGPX("Filename.gpx")
options(warn=0)
d <- d$tracks[[1]][[1]]
d$ele <- as.numeric(d$ele)
#d$time <- as.POSIXct(strptime(d$time, format="%Y-%m-%dT%H:%M:%S."))
library(berryFunctions)
colPoints(lon, lat, ele, data=d, add=FALSE, asp=1, lines=TRUE, pch=NA, lwd=3)
d2 <- equidistPoints(lon, lat, ele, data=d, n=25)
points(d2$x, d2$y, pch=3, lwd=3) # 25 segments equally along track
browseURL("http://www.r-bloggers.com/stay-on-track-plotting-gps-tracks-with-r")
# Reading KML files, eg from GoogleMaps:
##install.packages("maptools")
d <- maptools::getKMLcoordinates("Filename.kml",ignoreAltitude=TRUE)
d <- as.data.frame(do.call(rbind, d[which(sapply(d, class)=="matrix")]))
colnames(d) <- c("lon","lat")
head(d)
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 17: Character Stings (Zeichenketten) ###############################
# 17.1 Wilde Einführung
# 17.2 Suchen und ersetzen
# 17.3 Anwendungen Suchen und ersetzen
# 17.4 Zeichen regelmäßig ersetzen, Wildcards, Regex
# 17.5 Ausgabe in die Konsole
# 17.6 Anführungszeichen
# 17.7 in Buchformat drucken
# 17.8 Google PDF Links extrahieren
# 17.9 Wildes Ende (Verschiedenes)
# 17.1 Wilde Einführung --------------------------------------------------------
paste("WunderbareSache", 3:5) # macht einen Vektor mit Modus character
paste(2, "hot", 4, "U", sep=" ") # Standardeinstellung seperator: Leerzeichen
einVektor <- 5:9
paste("SeeYouL", einVektor[4], "ter", sep="") # Auf bestehende Objekte zugreifen
Buchstaben <- c("einige verschiedene", "Buch", "staben")
Buchstaben
length(Buchstaben)
nchar(Buchstaben)
paste(Buchstaben, "tolle", sep="_")
paste(Buchstaben, sep="")
paste(Buchstaben, collapse="")
# alle Leerzeichen ersetzen (mehr in 17.2 und 17.3)
Buchst <- paste(Buchstaben, collapse=""); Buchst
sub(" ", "", Buchst) # substitute (ersetze) " " in Buchst mit ""
# in R vordefiniert ist schon das alphabet als Objekt
letters
LETTERS
nchar("Monster Party") # Anzahl Zeichen; " " (Leerzeichen) ist auch ein Zeichen!
substring("Monster Party", first=c(1,3,5), last=c(2,5,9) ) # "t" zweifach drin
substr("Monster Party", start=c(1,3,5), stop=c(2,5,7)) # kann nicht vektorisiert
strsplit("Die Ärzte - Monster Party", " - ")[[1]]
s <- "am Donnerstag oder Freitag"
s1 <- strsplit(s, " ")[[1]] # split string into single words
s2 <- sapply(strsplit(s1, ""), rev) # reverse letters in each word
s3 <- sapply(s2, paste, collapse="") # combine to new words
paste(s3, collapse=" ") # paste into one single string
"ma gatsrennoD redo gatierF"
# Backslash "\" : Escape Character
# "\n" für Zeilenumbruch (newline), zB in Graphiküberschriften und -Texte
# "\t" für Tabstopp, zB für write.table("path/file.txt", yourdata, sep="\t")
# "\\" für Backslash als Zeichen
# "\"" für " als Zeichen
# "\U" für Unicode, zB.:
plot(1) # Unicode symbols without "expression"
mtext("\U00D8: 5.4 [\U00B0 C]", col=2) # average and degree symbol (with space)
mtext(expression(paste("Temp [",degree,"C]")), adj=0)
mtext("2 * 2\U00B2 = 2 \U00B3", col=2, adj=1)
browseURL("http://unicode-table.com/de/")
# control characters \a, \b, \f, \n, \r, \t, \v
# octal and hexadecimal representations like \040 and \0x2A
num <- c(3.2, 12.03, 789)
num
paste(num)# schneller zu tippen als as.character(num)
sprintf("%5.1f", num) # C-style String Formatting: 5 Zeichen, 1 nach dem Komma
sprintf("%6.1f", num)
sprintf("%7.2f", num)
sprintf("%07.2f", num) # pad with zeros (mit führenden Nullen füllen)
sprintf(" %07.2f", num) # Leerzeichen am Anfang
# 17.2 Suchen und ersetzen -----------------------------------------------------
schoko <- c("Lindor", "Ritter Sport", "Lindor","Frigor", "Milka", "More Lindor")
schoko
match("Milka", schoko) # wo kommt "Milka" in schoko vor? (exakter Match)
grep("or", schoko) # Elemente 1,2,3 und 6 von Schoko enthalten "or"
match("Lindor", schoko) # liefert nur das erste Element, das "Lindor" ist.
"Lindor" %in% schoko # logische Antwort
grep("Lindor", schoko) # liefert alle Teilmatches
# Ersetzen (SUBstitution)
sub("or", "__", schoko) # Letztes Lindor unverändert (First matches)
gsub("or", "_", schoko) # alle Matches, Global SUBstitution
sapply(gregexpr("or", schoko), function(x) sum(x>0)) # Wie oft komt "or" vor?
Kette <- "Lange ZeichenKette"
grep("g", Kette)
grep("u", Kette) # schwierig in Bedingungen etc zu verwenden, daher
grepl("g", Kette)
# grepl vektorisiert:
grepl("h|e", Kette)
grepl("u", Kette)
"u" %in% Kette # geht schneller
grepl("t", "D:/A_Data") # liefert TRUE, wie erwartet.
"t" %in% "D:/A_Data" # FALSE: %in% ist binärer Operator für match,
# match liefert nur die erste exakte (komplette) Übereinstimmung.
# an welchen Stellen einer Zeichenkette ist ein bestimmter Buchstabe:
unlist(gregexpr("e", Kette))
# 17.3 Anwendungen Suchen und ersetzen -----------------------------------------
# Eine Funktion soll den Speicherort einer Datei nur in die Console schreiben,
# wenn der Dateispeicherort nicht als kompletter Pfad genannt wurde, also
# das Working Directory (Arbeitsverzeichnis) als Ort verwendet wurde:
Datei <- "Beispiel.txt"
if( ! ( grepl("/", Datei) | grepl("\\", Datei, fixed=T) ) ) #"
cat("Hier ist deine Datei gespeichert:",getwd(),"\n")
Datei <- "D:/A Data/Beispiel.txt"
# analog zum obigen Vorgehen
Datei <- "D:\\A Data\\Beispiel.txt"
# analog zum analogen Vorgehen
# Number of occurences of character
# wie oft ein Buchstabe in einer Zeichenkette vorkommt
# elegantere Lösungen willkommen!
nocc <- function(pattern, x, ignore.case=FALSE)
sum(nchar(x)-nchar(gsub(pattern,"",x,ignore.case=ignore.case)))/nchar(pattern)
#
nocc(pattern="M", x="Moin an diesem schoenen Morgen") # 2
nocc("m", "Moin an diesem schoenen Morgen") # 1
nocc("M", "Moin an diesem schoenen Morgen", ign=T) # 3
nocc("m", "Moin an diesem schoenen Morgen", ign=T) # 3
nocc("n", "Moin an diesem schoenen Morgen") # 5
nocc(c("M","n"), "Moin an diesem schoenen Morgen") # nicht vektorisierbar, also
sapply( c("M","n"), nocc, x="Moin an diesem schoenen Morgen")
nocc("n", c("Moin an diesem", "schoenen Morgen"))
# gibt nur Gesamtsumme statt 2 und 3
# Lesetipp: Die Beispiele in
?chartr
# Quersumme (digit sum)
quersumme <- function(dat) # vektor mit (Komma-)Zahlen
sapply(strsplit(sub(".","",dat,fixed=T),""), function(x)sum(as.numeric(x)) )
#
quersumme(3.14159) # = 23 = 3+1+4+1+5+9
quersumme( c(9.37, 2.1) )
quersumme(1:30)
plot(quersumme(1:80), type="o") # ansteigende Periodizität
plot(quersumme(1:800), type="o") # Selbstähnlichkeit (Fraktalität)
plot(quersumme(1:8000), type="l")
# 17.4 Zeichen regelmäßig ersetzen, Wildcards, Regex ---------------------------
browseURL("http://www.regular-expressions.info/rlanguage.html")
x1 <- c("1.5 mg" , "50 mg" , "100 mg" , "0.5 mg") ; x1
x2 <- gsub(" mg", "", x1) ; x2 # ersetze alle " mg" in x1 mit "" (nichts)
as.numeric(x2)
# mit fixed=T in gsub kann man auch "[" etc (Sonderzeichen) ersetzen:
x3 <- c("Konz. [mg/l]" , "Masse [g]" , "Porosität [%]" , "Anzahl [-]") ; x3
gsub(" [g]", "", x3, fixed=TRUE)
# Ausflug: Introduction von Wildcards:
x4 <- c("Konz." , "Masse" , "Porosi tät" , "Anzahl", "Akt ") ; x4
grep("[t]$", x4) # die Indizes aller Elemente in x4 die auf "t" enden
# alternative Methode mit Wildcard im bekannten Stil (Platzhalter: "*")
grep(glob2rx("*t"), x4)
grep(glob2rx("*ss*"), x4) # Welche Elemente von x4 beinhalten "ss"
grep("ss", x4) # macht das gleiche, nur weniger kompliziert 😉
gsub(glob2rx("*t"), "tt", x4) # mit Wildcard ganzes Element ersetzen
gsub(".*t$", "tt", x4) # Hier steht das gleiche als Regular expression (?regex)
gsub("[t]$", "tt", x4) # nur t zu tt ändern
# zurück zu x3: Alle Einheiten entfernen (zB. um Objekte als Dateinamen zu
# verwenden); ohne glob2rx, also mit Wildcard in reg exp
x3
gsub(" \\[.*", "", x3) # "\\" um aus der Escape Sequenz zu escapen und "["
# analog zu fixed=T für exact matching zu nehmen; ".*" als Wildcard
gsub("ss.*", "lange,so sieht_manwas", x3)
# Alternativ mit strsplit
x5 <- strsplit(x3, split=" [", fixed=T) ; x5
unlist(x5)[(1:length(x5))*2-1]
# Sofern in den Zeichenketten keine weiteren "[" sind:
# Remove trailing zeros
x <- "0.710000"
sub("0+$", '', x)
# Remove leading and trailing white spaces
x <- " middle of the string "
sub("^[[:space:]]*(.*?)[[:space:]]*$", "\\1", x, perl=TRUE)
sub(" +$", '', x) # nur Leerzeichen am Ende der Zeichenkette entfernen
sub("^\\s+", "", x) # nur am Anfang
x <- c(" -some string" , "- leading minus sign")
sub("^[[-]]*", "\\1", x, perl=TRUE)
# 17.5 Ausgabe in die Konsole --------------------------------------------------
a <- print(c(1:3,"Hallo")); a
b <- cat(1:3,"Hallo") ; b
c <- paste(1:3,"Hallo") ; c
# print schreibt in die Konsole und gibt den Wert auch zurück
# cat schreibt nur in die Konsole, der zugeordnete value ist NULL
# cat zeigt keine Anführungsstriche!
# paste schreibt, wenn zugeordnet, nichts, sondern gibt einen zusammengefügten
# character Vektor zurück. Wenn nur der Befehl ausgeführt wird, gibt es das
# Ergebnis an die Konsole (wie sonst bei Funktionen auch immer der Fall).
# 17.6 Anführungszeichen -------------------------------------------------------
a <- "Ein \"schöner\" Text" ; a ; cat(a, "\n")
plot(1, main=a) # alles Roger
b <- "Ein 'noch schönerer' Text"; b ; cat(b, "\n") ; plot(1, main=b)
cc <- "Ein \"toller' Text"; cc ; cat(cc, "\n") #'
# http://tolstoy.newcastle.edu.au/R/help/02b/5198.html
# 17.7 Dokument in Buchformat drucken ------------------------------------------
# Broschürendruck zum Falten und Durchblättern wie im Buch.
# Dazu müssen die Seiten in der richtigen Reihenfolge (siehe Funktion) und
# Anordnung gedruckt werden (siehe Kommentar unten).
# Beispiel für 12 Seiten: 12,1,2,11,10,3,4,9,8,5,6,7
# Ungelesene Links:
# http://www.helpster.de/a5-buchdruck-in-word-vorbereiten-so-gelingt-s_54969
# http://www.ehow.com/how_4915094_print-booklet-adobe.html
Buchformatdruck <- function(n, sep=",") # n: number of pages in document
{ # sep: seperator required by printer
p <- ceiling(n/4)*4 # p: number of pages including empty pages at the back
d <- p/4 # d: number of double pages
o <- "" # o: order of pages
for(i in 1:d)
o <- paste(o, p-2*(i-1), 2*i-1, 2*i, p-2*(i-1)-1, sep=",")
cat(substr(o, start=2, stop=nchar(o)), "\n")
}
Buchformatdruck(12)
Buchformatdruck(13) # hier wären leere Seiten 14 bis 16 anzuhängen
# Beidseitig drucken, Drehung an langer Kante
# Word 2010:
# Seiten wie von Funktion angegeben
# 2 Seiten pro Blatt, Seite Einrichten -> gegenüberliegende Seiten
# Foxit Reader 4.3:
# Printer Properties - Querformat
# Print Range: Pages [Ergebnis der Funktion]
# Print Handling: # Margins 0.1 in # Page Scaling fit to paper
# Page Arrange multiple pages per sheet # rotate normal
# page order vertical # Pages per sheet custom 2 by 1 # Autocenter on
# 17.8 Google PDF Links extrahieren --------------------------------------------
# PDF-Links aus der Google-Suche sind ganz schön lang - ich hab noch keinen
# anderen Weg gefunden, nur den tatsächlichen Link zu erhalten.
# Die Links sehen zB so aus:
googlink <- paste0("http://www.google.de/url?sa=t&rct=j&q=semilog&source=web&c",
"d=7&ved=0CHgQFjAG&url=http%3A%2F%2Fwww.math.neu.edu%2Ffries%2Fsemilog-ex-hwk.",
"pdf&ei=8qnZT-qUKInCtAbU6-DFCA&usg=AFQjCNEwix4dJhFfWbKj86cBv5QeysgbVA&cad=rja")
library(berryFunctions)
googleLink2pdf # in berryFunctions
googleLink2pdf(googlink)
googleLink2pdf("http://www.google.de/url?sa=t&rct=j&q=gr%C3%BC%C3%9Fe%20pdf
&source=web&cd=1&ved=0CFkQFjAA&url=http%3A%2F%2Fwww.fischerverlage.de%2Fsixcms
%2Fmedia.php%2F308%2FLP_978-3-10-022303-6.pdf&ei=hq_ZT7TOKI7GtAaVs5HfBw&usg=
AFQjCNEfynh6ICrjiqI6D2phge7puMQzxQ&cad=rja")
# 17.9 Wildes Ende (Verschiedenes) --------------------------------------------
strReverse <- function(x)
sapply(lapply(strsplit(x, NULL), rev), paste, collapse="")
strReverse(c("abc", "Statistics"))
plot(1, type="n")
text(1,1, expression(rho ~ theta))
a <- c("A_A","A_B","C_A","BB","A_Asd")
a[grepl("^[AB]", a)]
# Interessantes aus der grep-und-co Hilfe mit "[" und "]":
grep("[a-k]", letters) # siehe Abschnitt character class in
?regex
# Thema Format:
as.character
l <- sprintf("%02.0f",1:12)
l
cat(sprintf("%s is %6.2f feet tall\n", "Sven", 7.1) )
## re-use one argument three times, show difference between %x and %X
xx <- sprintf("%1$d %1$x %1$X", 0:15)
xx <- matrix(xx, dimnames=list(rep("", 16), "%d%x%X"))
xx
noquote(format(xx, justify="right"))
sprintf("%s %d", "test", 1:3)
## Double all 'a' or 'b's; "\" must be escaped, i.e., 'doubled'
gsub("([ab])", "\\1_\\1_", "abc and ABC")
txt <- c("The", "licenses", "for", "most", "software", "are",
"designed", "to", "take", "away", "your", "freedom",
"to", "share", "and", "change", "it.",
"", "By", "contrast,", "the", "GNU", "General", "Public", "License",
"is", "intended", "to", "guarantee", "your", "freedom", "to",
"share", "and", "change", "free", "software", "--",
"to", "make", "sure", "the", "software", "is",
"free", "for", "all", "its", "users")
grep("[gu]", txt, value = TRUE)
sub("[b-e]",".", txt)
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel 18: Optimierung ####################################################
# 18.1 Exp Abnahme grob anpassen
# 18.2 Zu minimierende Größe: RMSE
# 18.3 Zu minimierende Funktion
# 18.4 Ergebnis darstellen
# 18.5 RMSE abhängig von k
# 18.6 Verschiedene Startwerte
# 18.7 verschiedene Methoden
# Parameter mit Optimierungsverfahren genauer schätzen -------------------------
# Exponentielle Temperaturabnahme im Kaffeebecher
# Anmerkung: mehrere Parameter mit optim, 1 einzelner mit optimize
# Daten
Zeit <- c( 0.0, 2.5, 6.0, 10.0, 12.5, 13.6, 17.3, 20.3, 23.9) # Minuten
Temp <- c(43.5, 41.6, 39.6, 37.6, 36.4, 36.2, 35.0, 34.1, 32.9) #°C
plot(Zeit, Temp, las=1, pch=16, ylab="Wassertemperatur [°C]", xlab="Zeit [min]")
# 18.1 Exp Abnahme grob anpassen -----------------------------------------------
T0 <- 44 # Anfangstemperatur
Te <- 25 # Endtemperatur
diff(range(Temp)) / diff(range(Zeit)) / ( mean(Temp)-Te)
k <- -0.036 # Abnahmerate [1/minute]
t <- seq(0, 30, by=0.1)
lines(t, (T0-Te)*exp(k*t)+Te ) # gar nicht mal so schlecht
param <- c(T0=44, Te=25, k=-0.036) # Anfangsparameter, später für optim
# Temperaturfunktion definieren:
temperaturfun <- function(p, t) { (p["T0"]-p["Te"])*exp(p["k"]*t)+p["Te"] }
temperaturfun(p=param, t=5) # nach 5 Minuten
lines(t, temperaturfun(p=param, t=t) , col=2 )
# gehen die gleichen Parameter wie vorhin rein
# 18.2 Zu minimierende Größe: RMSE ---------------------------------------------
# root mean square error: berechnet Abweichung von x1 zu x2, wird auch zur
# Berechnung des R2 in einer linearen Regression verwendet
rmse <- function(x1,x2)
{ if(length(x1) != length(x2)) stop("vectors not of equal length")
sqrt( sum((x1-x2)^2)/length(x2) )
} # Funktionsende
# Man könnte auch als andere Größe die Korrelation maximieren:
?cor()
rmse(Temp, temperaturfun(p=param, t=Zeit) ) # 0.53
cor(Temp, temperaturfun(p=param, t=Zeit) ) # 0.9986 - schon ziemlich gut
# graphisch zu sehen: geht aber noch viel besser
# 18.3 Zu minimierende Funktion ------------------------------------------------
# optim braucht eine Funktion mit par als erstem Input
minfun <- function(p) rmse(Temp, temperaturfun(p=p, t=Zeit))
minfun(param) # 0.53
optimized <- optim(par=param, fn=minfun)
optimized
# $value 0.11 Wir hatten angefangen mit 0.53, ist also deutlich verbessert
data.frame(initial=param, optimized=optimized$par)
# So viel hat sich nicht geändert - wir hatten ja auch sehr gute Startwerte
cor(Temp, temperaturfun(p=optimized$par, t=Zeit) ) # 0.9994, vorher 0.9986
# ist nicht so die ganz geeignete Größe dafür, liefert aber auch das Richtige:
maxfun <- function(p) cor(Temp, temperaturfun(p=p, t=Zeit))
optim(par=param, fn=maxfun, control=list(fnscale=-1) ) # fnscale<0: Maximierung
# 18.4 Ergebnis darstellen -----------------------------------------------------
plot(Zeit, Temp, las=1, pch=16, ylab="Wassertemperatur [°C]", xlab="Zeit [min]")
lines(t, temperaturfun(p=param, t=t) , col=2 )
lines(t, temperaturfun(p=optimized$par, t=t), col=1)
# deutlich sichtbar verbessert
legend("topright", c("observed", "initial", "optimized"), lty=c(NA,1,1),
pch=c(16,NA,NA), col=c(1,2,1), inset=0.05)
text(15 ,40.5, expression("y = (T0" - "Te) "%.%" "*e^("k "%.%" x")*" "+" Te"))
legend(14 ,40, paste(names(optimized$par), "="), bty="n", adj=1)
legend(14.5,40, round(optimized$par,3), bty="n", adj=0)
# 18.5 RMSE abhängig von k -----------------------------------------------------
k_values <- seq(-0.9, 0.025, len=200)
RMSerror <- rep(NA, length(k_values))
# mit den bereits optimierten anderen beiden Werten
for(i in 1:length(k_values))
RMSerror[i] <- rmse(Temp, (43.419-28.372)*exp(k_values[i]*Zeit)+28.372)
plot(k_values, RMSerror, type="l", las=1)
# Ist jetzt aber nur eine eindimensionale Betrachtung!
# 18.6 Verschiedene Startwerte -----------------------------------------
vgl_startw <- function(T0=44, Te=25, k=-0.036)
{param <- c(T0=T0, Te=Te, k=k)
optimized <- optim(par=param, fn=minfun)
list(data.frame(initial=param, optimized=optimized$par), RMSE=optimized$value)
} # Funktionsende
vgl_startw()
vgl_startw(k=-0.2)
vgl_startw(k=-0.9) # k optim = -10.6 -> geht nicht mehr auf das richtige k zu
# Zeigt, dass Startwerte sinnvoll gesetzt werden müssen!!
# Systematisch Grenzen der Methode ausloten
k_test <- seq(0.5, -1, length=5000)
sense <- matrix(0, nrow=length(k_test), ncol=2); colnames(sense)<- c("k","RMSE")
head(sense)
total <- length(k_test)
pb <- winProgressBar(title="progress bar", min=0, max=total)
# Fortschrittsbalken hilft in langsamen for-Schleifen, evtl in Sichtweite ziehen
# Linux/Mac: tkProgressBar im Paket tcltk
for(i in 1:total)
{sense[i,1] <- k_test[i]
sense[i,2] <- vgl_startw(k=k_test[i])[[1]][3,2]
setWinProgressBar(pb, i, title=paste( round(i/total*100, 0), "% done"))
if(i==total) {Sys.sleep(2) ; close(pb) }
} # Ende for Schleife
windows()
plot(sense, type="l", las=1, xlab="initial k", ylab="optimized k")
abline(h=-0.049, col=2, lty=2)
text(-0.8, -0.049, "real k", col=2)
# 18.7 verschiedene Methoden ---------------------------------------------------
# parametersatz mit nur k-Startwert falsch
optim(par=c(T0=43.42, Te=28.37, k=-0.9), fn=minfun, method="Nelder-Mead")$par
# Jetzt richtig, obwohl vorhin mit -0.9 falsch. (Te besser gewählt)
# Findet echten Wert?
optim(par=c(T0=44, Te=29, k=-0.9), fn=minfun, method="Nelder-Mead")$par # nein
optim(par=c(T0=44, Te=29, k=-0.9), fn=minfun, method="BFGS")$par # ja
optim(par=c(T0=44, Te=29, k=-0.9), fn=minfun, method="CG")$par # nein
optim(par=c(T0=44, Te=29, k=-0.9), fn=minfun, method="SANN")$par # nein
# "L-BFGS-B" und "Brent" gehen hier nicht
# Bei mehr Daten ist die Wahl der Startwerte weniger problematisch
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
### RclickHandbuch.wordpress.com ### Berry Boessenkool ### berry-b@gmx.de ######
### Kapitel V (Verschiedenes): Rprofile und weitere Tipps ######################
# V.1 Rprofile.site
# V.2 Verschiedene Tipps
# V.3 Workspace speichern
# V.4 Pakete # ToDo: Info aus dem gesamten Handbuch hier bündeln
# V.5 Logische Operatoren
# V.6 Gleichheit Kommazahlen
# V.7 Fehlermeldungen
# ToDo: Sammlung Fehler und Lösungen (vielleicht in ein früheres Kapitel?)
# V.1 Rprofile.site ------------------------------------------------------------
# TinnR kann mehrere Zeilen Code (bis zur nächsten Leerzeile) mit einem Klick
# an R senden (R send: contiguous), braucht aber besondere Pfade dazu:
.trPaths <- paste(paste(Sys.getenv('APPDATA'), '\\Tinn-R\\tmp\\', sep=''),
c('', 'search.txt', 'objects.txt', 'file.r', 'selection.r',
'block.r','lines.r'),sep='')
# Praktisch ist es, die Zuordnung dauerhaft an R zu übergeben. Kopiere ihn dazu
# in "Rprofile.site" (im R-Programmordner (wo R installiert ist) unter "etc").
# Die Datei kann mit Editor oder Äquivalenten geöffnet werden.
# Nach dem nächsten R Start können mehrere Zeilen auf einmal gesendet werden.
# Siehe auch RclickHandbuch.wordpress.com/install-r/rprofile
# Auch eigene Zuordnungen, zB Funktionen, lassen sich im Rprofile
# speichern, und müssen dann nicht mehr jedes Mal an R geschickt werden.
# Das alles geht leider nur, wenn Admin-Rechte vorliegen.
# Falls nicht: erstelle ein Skript mit allen Zuweisungen mit einem sinnvollen
# Namen (zB Funktionen.r), und speichere das auf der Festplatte. Dann muss nur
# am Anfang jeder Sitzung einmal source("D:/Funktionen.r") ausgeführt werden.
?Rprofile
# Wichtig: Sachen in Rprofile überschreiben ls(.BaseNamespaceEnv)!
# Schreibe nichts da rein, was es im base von sich aus gibt!
# Praktisch: Dateipfade in Rprofile speichern
desktop <- ("C:/Dokumente und Einstellungen/berry/Desktop")
# Danach in jedem Skript möglich:
setwd(desktop)
# Setze den CRAN mirror (Wo Pakete runtergeladen werden)
local({r <- getOption("repos")
r["CRAN"] <- "http://mirrors.softliste.de/cran" # Berlin
options(repos=r)})
# leute richtig verwirren (heimlich in deren Rprofile schreiben):
pi <- 3
library(fortunes); fortune(352)
# oder richtig dissen:
q("no") # Bitte nicht ausprobieren, wenn grad viele Objekte zugeordnet sind...
# Graphiken standardmäßig mit den y-Achsenwerten aufrecht erstellen
setHook("plot.new", function(...) par(las=1), "append")
# Geht aber noch nicht. Ideen vor!
# Note to self:
?windows.options # Examples
# Mehr unter RclickHandbuch.wordpress.com/install-r/rprofile
# V.2 Verschiedene Tipps -------------------------------------------------------
# Viele Tipps zum angenehmen Programmieren via TinnR finden sich hier:
# https://rclickhandbuch.wordpress.com/install-r/tinn-r/
# Lange Berechnungen in der R-Konsole können mit ESC abgebrochen werden.
# Diese sind daran sichtbar, das noch kein ">" in der letzten Zeile erscheint.
round(17.8)
as.integer(17.8) # hättst jetzt 18 erwartet, wa?
#
x <- -1:12 # %/%: ganze Teiler; %%: Rest (modulo operation)
data.frame(x, Teiler_5=x%/%5, Rest=x%%5 )
# Spaltennamen, die mit einer Zahl beginnen, wird ein X vorgesetzt
data.frame("2.7.2010"=12:18, ZweiteSpalte=18:12)
# Das passiert auch bei mit read.table geladenen Dataframes!
# Es kann aber ersetzt werden (praktisch zB für Subs von Boxplots):
# names(Objekt) <- c("1.Tag", "2. Messung") # sogar Leerzeichen sind OK!
# sourcecode von Funktionen:
lm
summary
methods(summary)
summary.table
# Nachkommastellen (der Ausgabe in der RConsole) einstellen
d = data.frame( a=rnorm(10) , b=rnorm(10) )
options(digits=4) # Übersichtlichkeit der Datendarstellung einstellen
d
print(d[1,1], digits=18) # Die Daten selbst haben noch viele Nachkommastellen.
options(digits=7) # Standardeinstellung
# Umgang mit NAs (Not available, unbekannt)
v <- c(7,9,NA,5,6)
cumsum(v) # 7 16 NA NA NA
# Inhaltlich falsch, aber wenn man davon ausgeht, dass NAs 0 entsprechen:
# Ersetze in v alle Werte == NA mit 0 (replace)
cumsum( replace(v, is.na(v), 0) ) # 7 16 16 21 27
# Alternative: Mittelwerte der nicht-NA-Werte:
cumsum(replace(v,is.na(v),mean(v, na.rm=T))) # 7.00 16.00 22.75 27.75 33.75
# V.3 Workspace speichern ------------------------------------------------------
# Beim Schließen von R braucht das WorkSpace normalerweise nicht gespeichert zu
# werden, es kann ja mittels des Skripts in TinnR erneut hergestellt werden.
# Spätestens nach vierstündigen Berechnungen sollte es aber gespeichert werden.
# siehe dazu
save(list="Objekt", file="Dateiname.RData")
?save.image # Für das gesamte Workspace, also kurz für save(list=ls(), ".RData")
# mit rm("Objekt7") kann man Objekte einzeln löschen und dann save.image nutzen.
load("ewaigeSubordner/Dateiname.RData")
?history
# Nach solchen großen Berechnungen kann es sinnoll sein, R zu schließen und neu
# zu öffnen. Siehe dazu die relative Größe der Auslagerungsdatei im Taskmanager
# (STRG + ALT + ENTF), (unter Linux in der Systemüberwachung)
# V.4 Pakete -------------------------------------------------------------------
library() # zeigt alle runtergeladenen und mitinstallierten packages
library("datasets")
library(help="datasets")
# Etwas zu Paketen, Vgl auch Kapitel 4.12
library(alr3) # geht noch nicht
install.packages("alr3") # nur einmal abschicken! Nicht in RProfile stecken!
# jetzt wurde das Paket gedownloaded und installiert.
# Aus Recheneffizenzgründen musss es aber nach jedem R-neustart geladen werden.
library(alr3)
# in eigenen Funktionen sollte require statt library benutzt werden.
# Dann werden Fehlermeldungen a la "Wurde unter R 2.13 erstellt" unterdrückt.
# V.5 Logische Operatoren ------------------------------------------------------
2>4
!2>4 # NICHT
A <- c(TRUE, TRUE, FALSE, FALSE) ; B <- c(TRUE, FALSE, TRUE, FALSE)
?"&"
data.frame(A, B, Nicht_A = !A, A_und_B = A&B, A_oder_B = A|B)
dat <- rnorm(300) ; hist(dat)
dev.off()
hist(dat, plot=F)
F
F <- 3
F
hist(dat, plot=F)
# F geht nur als Abkürzung für FALSE, solange es nicht überschrieben wird.
rm(F)
F
# Viele Anwendungen von TRUE/FALSE auch in Bedingungen --> Kap 11.1: if
dat <- sample(-6:10)
dat > 4
which(dat > 4) # An welchen Stellen eines logical Vectors steht TRUE?
sum(dat > 4) # Anzahl TRUE. Funktioniert, da TRUE = 1 und FALSE = 0
dat < -2 # Leerzeichen wichtig! dat <-2 würde "dat" mit der Zahl 2 überschreiben
dat<(-2) # ist schlecht menschenlesbar, aber eindeutig.
# V.6 Gleichheit Kommazahlen ---------------------------------------------------
0.5-0.2 == 0.3 # ja nee, is klar. Denkst du?! Dann check mal das:
0.4-0.1 == 0.3 # Grundgedanke aus "The R-Inferno" entnommen - sehr lesenswert!
# Zweiteres zeigt FALSE, obwohl es richtig ist! Folgendes ist der Grund:
# ein Rechner kann Fließkommazahlen nur endlich genau verarbeiten - arbeitet
# mit binären Zahlen - abhängig von bit-Betriebssystem, Prozessor-Architektur,
# definition des Zahlenformats (single = 1 byte = 8bit, double = 2 byte), usw.
# Am Rande für Informatiker: R verwendet nur Double, siehe Note zu ?numeric.
print(0.4, digits=20)
print(0.1, digits=20) # Schlagwort: arithmetische Präzission
print(0.3, digits=20)
print(0.4-0.1 , digits=20)
# Folgendes ist die Lösung für einzelne Zahlen:
all.equal( 0.4 - 0.1 , 0.3 ) # TRUE (lässt PC-Ungenauigkeit zu)
# Aber nicht für ganze Vektoren:
set.seed(12); zahlen <- sample( 1:8+3:6/10, 30, replace=T) ; zahlen
all.equal(zahlen, 4.6) # geht nicht
# Wenn man z.B. das Ergebnis einer anderen Rechnung verwenden will:
zahlen == 4.9-0.3 # Alles FALSEs, obwohl es 4.6 gibt!
4.9-0.3 == 4.6 # Absolute Gleichheit nicht gegeben
# Lösung für Zahlenvergleich von Vektoren:
# Komplizierte (aber allgemeingültige) Lösung:
TF_Werte_einzeln <- by(zahlen, 1:length(zahlen), FUN=all.equal, current=4.9-0.3)
welcheT <- which(TF_Werte_einzeln==T); welcheT
str(TF_Werte_einzeln)
str(welcheT) # welcheT hat names, aber das stört ja nicht
names(welcheT) <- NULL ; welcheT # Falls es doch mal stört 😉
# einfache Lösung:
which( zahlen == round(4.9-0.3, 8) ) # Problematisch, wenn mehr
# Nachkommastellen wichtig sind, und man das nicht bemerkt!
# Randnotiz: Runden auf signifikante Stellen, also Zahlen nach Nullen:
round( 0.00314869, digits=3)
signif(0.00314869, digits=3)
# ANzahl notwendiger digits beim Ausgangsbeispiel:
round(0.4-0.1, 20) == round(0.3, 20) # noch FALSE
round(0.4-0.1, 15) == round(0.3, 15) # TRUE
# Merke: niemals Fließkommazahlen einfach so auf absolute Gleichheit prüfen!
# V.7 Fehlermeldungen ----------------------------------------------------------
# Unterschied Fehler und Warnung:
# Fehlermeldungen treten bei Fehlern auf (Ach nee). --> Abbruch der Funktion.
# Warnmeldungen sollten gelesen werden, weil die Funktion nicht abgebrochen,
# sondern wahrscheinlich anders ausgeführt wird, als du willst!
# häufige Fehlermeldungen:
# "Unerwartete(s) Symbol"
# Meistens eine Klammer noch nicht wieder geschlossen, ein Komma vergessen, etc.
# Sammlung to be ausgedehnt.
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
################################################################################
R-click Handbuch 
0 Antworten to “AllesZurSuche”