#......Εντολές από το αρχείο R-guide.doc.............. #---------Σύνταξη εντολών στην R---------- x <- 3.1 #απόδοση τιμής x = 3.1 #απόδοση τιμής 3.1 -> x #απόδοση τιμής x 3 * x x <- 3 * x # αντικατάσταση τιμής x a <- 2.45 ; b <- 6.315 # initial values 1-pi+sin(pi/3) g <- (1-pi)*cos(pi/3); 2 - .Last.value # χρήση τελευταίας τιμής .Last.value gamma(1:10) # υπολογισμός διανύσματος γνωστής συνάρτησης Γ(n+1)=n! factorials <- .Last.value facorials #-----------Διανύσματα και πίνακες----------------- grades <- c(29.3,34.5,22.7,31.8,27.3,28.5) # δοθέν διάνυσμα grades grades[2] # η δεύτερη συνιστώσα grades > 30 # λογικό διάνυσμα F αν δεν ικανοποιείται η συνθήκη T αλλιώς students.names <- c("Alfred", "Denis", "Barbara", "John", "Mary", "Nick") # αλφαριθμητικό διάνυσμα ονομάτων names(grades) <- students.names # προσθήκη ονομάτων στο διάνυσμα grades # το διάνυσμα με ονόματα συνιστωσών names(grades) # έλεγχος και εξαγωγή των ονομάτων grades["John"] # έλεγχος στοιχείου με το όνομά του grades[grades > 30] # συρρίκνωση διανύσματος με συνθήκη (επιλογή στοιχείων) length(grades) # μήκος διανύσματος class(grades) # έλεγχος τύπου αντικειμένου class(names(grades)) # έλεγχος τύπου αντικειμένου class(c(29.3, 34.5, 22.7, 31.8, 27.3, 28.5)) # έλεγχος τύπου αντικειμένου class(grades > 30) # έλεγχος τύπου αντικειμένου let <- letters[1:5] # τα 5 πρώτα πεζά γράμματα ως διάνυσμα # συνένωση της λέξης letter με τους αριθμούς 1,2,...,5 χωρίς κενά # και ονοματοδοσία των στοιχείων του διανύσματος let names(let) <- paste("letter", 1:5, sep = "") ; let num <- 1:5 # αριθμητικό διάνυσμα με στοιχεία 1,2,...,5 # συνένωση της λέξης number με τους αριθμούς 1,2,...,5 χωρίς κενά # και ονοματοδοσία των στοιχείων του διανύσματος num names(num) <- paste("number", 1:5, sep = "") ; num grades <- as.vector(grades) ; grades # αφαίρεση ιδιοτήτων όπως ονόματα # απόδοση ονομάτων της μορφής student.1, student.2, ... names(grades) <- paste("student", 1:6, sep = ".") grades #έλεγχος με ονόματα dim(grades) <- c(2,3) # μετατροπή σε πίνακα δοθείσης διάστασης. #Τα 2 πρώτα στοιχεία αποτελούν την πρώτη στάθμη (1η στήλη) της δεύτερης διάστασης #Τα 2 επόμενα στοιχεία αποτελούν τη δεύτερη στάθμη (2η στήλη) της δεύτερης διάστασης grades # έλεγχος # [,1] [,2] [,3] #[1,] 29.3 22.7 27.3 #[2,] 34.5 31.8 28.5 grades[2, ] # η δεύτερη γραμμή του πίνακα grades[, 3] # η τρίτη στήλη του πίνακα grades[2, 3] # το στοιχείο (2,3) του πίνακα dim(grades) <- NULL # Η ακύρωση της προηγούμενης μετατροπή matrix(grades, 2, 3) # άλλος τρόπος μετατροπής σε πίνακα matrix(grades, 2, 3, byrow=T) # ή αν θέλουμε κατά γραμμές # [,1] [,2] [,3] #[1,] 29.3 34.5 22.7 #[2,] 31.8 27.3 28.5 c(T, F, pi, 7) # παράθεση ομοειδών αντικειμένων και μετατροπή λογικών τιμών #------------- Λίστες ------------------------------- sc1 <- c(27,31,30) ; sc2 <- c(28,38,32,40) # ορισμ΄ός διανυσμάτων sc3 <- c(28,20,20) ; sc4 <- c(30,35,35,27) sc5 <- c(36,30,16) ; sc6 <- c(27,32,25,30) Class <- list(students=students.names, age=c(12,13, 12,12,13,12), tests=c(3, 4, 3, 4, 3, 4), scores=c(sc1, sc2, sc3, sc4, sc5, sc6)) # ορισμός λίστας με 4 συνιστώσες Class # έλεγχος Class$students # η συνιστώσα students Class$st # ίδιο με προηγούμενο, αρκεί τα st να καθορίζουν μόνο μία συνιστώσα Class$age[3] # το 3ο στοιχείο της συνιστώσας age places <- c("Thes","Kav","Thes","Thes","Kat","Gian") Class <- c(Class, list(Place.birth=places)) # προσθήκη συνιστώσας σε λίστα unlist(Class) # μετατροπή της λίστας σε διάνυσμα unlist(Class, use.names=F) # το ίδιο χωρίς ονόματα dim(grades) <- c(2,3); grades dimnames(grades) <- list(paste("row",letters[1:2]), paste("col",LETTERS[1:3])) ; grades # ονόματα διαστάσεων πίνακα dimnames(grades) <- list(NULL, paste("Τάξη",c("Α","Β","Γ"))) # ομοίως grades #--------------Πλαίσια δεδομένων (Data Frames) ------------------- iris # 150 επί 5 πλαίσιο δεδομένων iris[c(1:3,147:150), , ] # πρώτες και τελευταίες γραμμές του names(iris) # ονόματα στηλών-μεταβλητών #[1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species" z<-iris$Sepal.Width # η δεύτερη μεταβλητή z<-iris[[2]] # το ίδιο c(mean=mean(z),st_dev=sd(z)) # μέση τιμή και την τυπική απόκλιση της δεύτερης μεταβλητής # mean st_dev #3.0573333 0.4358663 table(iris$Species) #πίνακας συχνοτήτων ως προς το είδος των λουλουδιών attach(iris) # κάνει τα ονόματα του πλαισίου ενεργά x1=Sepal.Length[1:50] # οι τιμές της Sepal.Length για το είδος setosa x1=Sepal.Length[Species=="setosa"] # άλλος τρόπος x2=Sepal.Length[51:100] # οι τιμές της Sepal.Length για το είδος versicolor x3=Sepal.Length[101:150] # οι τιμές της Sepal.Length για το είδος virginica summary(x1) # περιγραφικά μέτρα θέσης # Min. 1st Qu. Median Mean 3rd Qu. Max. # 4.300 4.800 5.000 5.006 5.200 5.800 summary(x2) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 4.900 5.600 5.900 5.936 6.300 7.000 summary(x3) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 4.900 6.225 6.500 6.588 6.900 7.900 #-------------------Πολλαπλοί πίνακες (arrays) --------------------- Aarray <- array(c(1:8, 11:18, 111:118), dim = c(2, 4, 3)) ;Aarray # για την στάθμη 1 της τρίτης διάστασης (1ος όροφος) τα πρώτα 2 επί 4=8 στοιχεία # του διανύσματος σχηματίζουν πίνακα διπλής εισόδου. # Η πρώτη στάθμη της δεύτερης διάστασης (1η στήλη) αποτελείται από τα # πρώτα 2 στοιχεία, η 2η στήλη από τα άλλα δύο στοιχεία κ.ο.κ. # Το ίδιο για το 2ο όροφο με τα επόμενα 8 στοιχεία και # τέλος ο 3ος όροφος με τα τελευταία 8 στοιχεία Aarray[, , 2] # ο δεύτερος όροφος # [,1] [,2] [,3] [,4] #[1,] 11 13 15 17 #[2,] 12 14 16 18 Aarray[1, , ] # οι πρώτες γραμμές και στους τρεις ορόφους # [,1] [,2] [,3] #[1,] 1 11 111 #[2,] 3 13 113 #[3,] 5 15 115 #[4,] 7 17 117 x<-as.matrix(iris[,1:4]) # μετατροπή σε πίνακα που διαβάζεται κατά στήλες # δηλαδή το x ως διάνυσμα περιέχει 150*4=600 τιμές με πρώτες τις 150 τη; # πρώτης στήλης, μετά τις 150 της δεύτερης στήλης κλπ y<-array(x,dim=c(50,3,4)) # μετατροπή σε array, κάθε στήλη = 50*3 y[1:2,,] # Οι δύο πρώτες (από τις 50) γραμμές, όλων των στηλών και ορόφων # αλλαγή ονομάτων διαστάσεων ns=names(table(Species)) nm=names(iris)[-5] # αφαιρέθηκε το όνομα της πέμπτης στήλης y<-array(x,dim=c(50,3,4),dimnames=list(NULL,ns,nm)) y[1:2,,] apply(y,2:3,FUN=mean) # τ 2¨3 καθορίζει τις διαστάσεις #---------- Παράγοντες ------------------------ smok <- factor(c("Sm","No-Sm","No-Sm","NA","Sm","NA","Sm")) #ορισμός παράγοντα smok print.default(smok) # εκτύπωση με κωδικούς attributes(smok) # εκτύπωση ιδιοτήτων smoka <- factor(c("Sm", "No-Sm", "No-Sm","NA", "Sm","NA","Sm"), levels=c("No-Sm","Mid-Sm","Sm")) # καθορισμός σταθμών smoka print.default(smoka) income <- ordered(c("Mid","Lo","Lo","Hi","Mid", "Hi","Hi","Lo"), levels=c("Lo","Mid","Hi")) # διαταγμένος παράγοντας income #---------------------------------------------------------- attach(Class) # ελέγχουμε τις υπάρχουσες συνιστώσες students age tests scores Place.birth # Υπολογισμός βαθμών από τα τεστ g1=round(mean(scores[1:3]),1);g1 g2=round(mean(scores[4:7]),1);g2 g3=round(mean(scores[8:10]),1);g3 g4=round(mean(scores[11:14]),1);g4 g5=round(mean(scores[15:17]),1);g5 g6=round(mean(scores[18:21]),1);g6 grades=c(g1,g2,g3,g4,g5,g6);grades pass <- factor(grades > 30, labels = c("FAIL","PASS")) ; pass # παράγοντας για επιτυχία stds <- data.frame(tests, grades, pass, age, Place.birth); stds # ορισμός πλαισίου δεδομένων stds <- data.frame(n.of.tests = tests, grades, pass, age, birth.places = Place.birth) ; stds # με καθορισμό ονομάτων στηλών row.names(stds) <- students ; stds # καθορισμός ονομάτων γραμών-περιπτώσεων # n.of.tests grades pass age birth.places #Alfred 3 29.3 FAIL 12 Thes #Denis 4 34.5 PASS 13 Kav #Barbara 3 22.7 FAIL 12 Thes #John 4 31.8 PASS 12 Thes #Mary 3 27.3 FAIL 13 Kat #Nick 4 28.5 FAIL 12 Gian #--- Αρχεία, Workspace και Working Directory --- ls() # καταγραφή ενεργών αντικειμένων rm(foo, bar) # διαγράφονται τα αντικείμενα foo, bar remove(income, smok, smoka) # διαγράφονται τα αντικείμενα income, smok, smoka ls() # έλεγχος διαγραφής x=c(1.2,2,5,7,8) dat <- data.frame(x=c(1:10,1:10), y=1:20) attach(dat) x+y rm(x) x+y getwd() #[1] "c:/Users/Chronis/Desktop/My R Dir" setwd("f:/temp") getwd() #[1] "f:/temp" > bathmoi<-read.table("bathmoi.dat");bathmoi sex school g1 g2 g3 g4 pref stud001 M 15 6.2 12.5 19.6 17.6 3 stud002 M 12 12.8 15.0 17.2 19.9 1 stud003 F 5 7.2 9.5 11.6 13.1 9 stud004 M 11 15.5 14.3 10.9 16.4 2 stud005 F 11 19.2 18.6 15.9 18.6 1 stud006 M 5 17.6 14.0 10.5 8.2 6 library(foreign) ?read.spss read.spss(file.sav) library(gdata) ?read.xls read.xls(lile.xls) if(mode(x)!="character" && min(x)>0) log(x) # εκτελειίται αν όλα τα x>0 και αριθμοί x x <- cos(((1:10) * pi)/3) ; x if(mode(x) != "character" && x > 0) log(x) # εκτελειίται μόνο για x>0 και αριθμοί # ...... παράδειγμα χρήσης του τελεστή || ............... k=0 for (i in 1:1000) { x <- round(runif(100),2) y <- round(runif(100),2) sort(x) sort(y) if (any(x == 0.55) || any(y == 0.55)) { ".55 encountered" k=k+1} } k # ............................... x <- 1:16 ; x y <- 1:10 ; y x + y x <- 1:16 ; x y <- 1:8 ; y x + y 2 ^ y - 1 #---- Οι δικές μας Συναρτήσεις --------------------- msq <- function(x){ x0 <- sqrt(x) ; x1 <- sqrt (x+0.5) y <- asin(x1)-asin(x0) mean(y) } msq(c(.25,.32,.37,.47,.5)) msq <- function(x) { if(min(x) < 0 || max(x) > 0.5) { stop("vector x is outside (0,.5)") } x0 <- sqrt(x) x1 <- sqrt(x + 0.5) y <- asin(x1) - asin(x0) mean(y) } msq(c(.25,.32,.37,.47,-.5)) TwoIndPop <- function(x1, x2, dif=0) { n1 <- length(x1); n2 <- length(x2) mx1 <- mean(x1); mx2 <- mean(x2) vx1 <- var(x1); vx2 <- var(x2) vpooled <- ((n1-1)*vx1 + (n2-1)*vx2)/(n1+n2-2) tstat <- (mx1 - mx2-dif)/sqrt(vpooled*(1/n1 + 1/n2)) tstat } x1=c(100,105, 102,96,106,95,110,120, 115,112,112,90);x1 x2=c(104,98,100,90,110,116,120,115,95,88);x2 TwoIndPop(x1,x2) qt(.95,length(x1)+length(x2)-2) [1] 1.724718 # ---- Συστηματικές ακολουθίες αριθμών ------------------ k <- 6; 1:k-1; 1:(k-1) # οι πέντε επόμενες εντολές δίνουν το ίδιο διάνυσμα seq(5,11,2) seq(5,11,,4) seq(,11,2,4) seq(5,,2,4) seq(5,,2,,10:13) x <- c(1,3,4,7) # τίθεται x=(1, 3, 4, 7) y <- rep(x,2) # τίθεται y=(1, 3, 4, 7, 1, 3, 4, 7) w <- rep(x,,6) # τίθεται w=(1, 3, 4, 7, 1, 3) i <- rep(2,4) # τίθεται i=(2, 2, 2, 2) z <- rep(x,i) # τίθεται z=(1, 1, 3, 3, 4, 4, 7, 7) u <- rep(x,1:4) # τίθεται y=(1, 3, 3, 4, 4, 4, 7, 7, 7, 7) y <- c(0.31, 0.45, 0.46, 0.43, 0.82, 1.1, 0.88, 0.72, 0.43, 0.45, 0.63, 0.76, 0.45, 0.71, 0.66, 0.62, 0.36, 0.29, 0.4, 0.23, 0.92, 0.61, 0.49, 1.24, 0.44, 0.35, 0.31, 0.4, 0.56, 1.02, 0.71, 0.38, 0.22, 0.21, 0.18, 0.23, 0.3, 0.37, 0.38, 0.29, 0.23, 0.25, 0.24, 0.22, 0.3, 0.36, 0.31, 0.33) z1 <- rep(1:3, each = 16) ; z1 z2 <- rep(rep(1:4, 3), each = 4) ; z2 y[(z1 == 2) & (z2 == 3)] x1 <- ceiling(1:48/16) x2 <- 1+(ceiling(1:48/4)-1) %% 4 # ---- Τυχαίοι Αριθμοί ------------------ runif(100,1,3) #100 τυχαίους αριθμούς από την ομοιόμορφη κατανομή στο διάστημα (1,3) mean(rchisq(50,12)) #μέση τιμή 50 τυχαίων τιμών από την κατανομή χ^2 με 12 βαθμούς ελευθερίας set.seed(1949) # τίθεται συγκεκριμένο τυχαίος σπόρος mean(rchisq(50, 12)) # μετά την προηγούμενη θα δώσει "13.14979" dnorm(1.25,0,2) # η πιθανότητα P(Y=1.25)=0.1640805, όπου Y έχει την Ν(0,22) κατανομή pnorm(1.25) # η πιθανότητα P(Z<1.25), όπου η Z έχει την τυπική κανονική κατανομή qt(.85,10) # το 85ο εκατοστιαίο σημείο της κατανομής t_10 με 10 βαθμούς ελευθερίας # ---- Διάταξη ------------------ x <- c(3,5,9,1,7,4,15,4,3,9,7) ; x # ένα διάνυσμα sort(x) # διατάσσει τα στοιχεία διανύσματος x σε αύξουσα σειρά sort(names(grades)) # διατάσσει τα στοιχεία διανύσματος αλφαριθμητικά names(x) <- paste("x",1:length(x),sep="") ; x # ονοματίζουμε τα στοιχεία του x sort(x) # τα ονόματα ακολουθούν τη διάταξη order(x) # δίνει τους δείκτες του αρχικού x στη διάταξη x[order(x)] # Παρατηρείστε ότι ισχύει x[order(x)]=sort(x) x[order(x)] == sort(x) # καλύτερος έλεγχος του προηγούμενου Ένα άλλο διάνυσμα y ίδιου μήκους με το x διατάσσεται σύμφωνα με τη διάταξη στο x με την εντολή y[order(x)]. Παραδείγματος χάριν, αν y είναι το διάνυσμα μήκους όσο το x, που διαμερίζει το [1,10] σε ίσου εύρους διαστήματα, τότε: y <- seq(1, 10, , length(x)) ; y # διάταξη του y σύμφωνα με τη διάταξη στο x y[order(x)] rev(sort(x)) # διατάσσει τα στοιχεία του x σε φθίνουσα σειρά. rank(x) # διατακτικοί αριθμοί (αν ο 2ος και 3ος είναι ίσοι αντιστοιχούν στο 2.5) # [1] 2.5 6.0 9.5 1.0 7.5 4.5 11.0 4.5 2.5 9.5 7.5 sort(rank(x)) # [1] 1.0 2.5 2.5 4.5 4.5 6.0 7.5 7.5 9.5 9.5 11.0 diff(x) # πρωτες διαφορές # Ισχύει diff(x)=x[-1]-x[-length(x)] z <-c("L","M","M","H","M","H","L","H","M") match(z,c("L","M","H")) # συγκρίνει τιμές δύο διανυσμάτων # ----- Ανάγνωση Δεδομένων ------------------------------------ x <- scan() # περιμένει το πρόγραμμα να διαβάσει δεδομένα που θα τα αποδώσει στο x 2.3 12 3.5 456 34 23.1 2 12 # τελειώνουμε με μία κενή γραμμή x # έλεγχος x <- scan(what="") # εισαγωγή αλφαριθμητικού διανύσματος 1: f F mm 3F 4 F1 7: > x heart.list <- scan("heart.dat", skip = 1, what = list(name = "", group = 0, x1 = 0, x2 = 0)) # ανάγνωση λίστα από αρχείο heart.list # έλεγχος heart1.list <- scan( what=list(name="", group=0,x1=0,x2=0)) # εισαγωγή λίστας # μετά την προηγούμενη δίνουμε π.χ. John 1 46 18.2 Nick 1 49 21.3 Mary 1 32 16.4 Helen 2 58 20.7 Teo 2 25 18.4 Peter 2 38 20.5 heart1.list # έλεγχος mat <- matrix(scan( ), byrow=T, ncol=2) # εισαγωγή πίνακα # μετά την προηγούμενη δίνουμε π.χ. 450 54.6 760 73.4 325 50.3 285 58.1 mat # έλεγχος names<-c("John","Nick","Mary","Helen","Teo", "Peter") colmns<-c("Group","x1","x2") heart.df <- matrix(scan(),ncol=3,dimnames=list(names,colmns)) # εισαγωγή πλαισίου δεδομένων # μετά την προηγούμενη εντολή περιμένει τα τρία διανύσματα για τις στήλες 1 1 1 2 2 2 46 49 32 58 25 38 18.2 21.3 16.4 20.7 18.4 20.5 heart.df # έλεγχος Group x1 x2 John 1 46 18.2 Nick 1 49 21.3 Mary 1 32 16.4 Helen 2 58 20.7 Teo 2 25 18.4 Peter 2 38 20.5 fix.d <- as.list(read.fwf("fixwidth.dat", widths=c(1,3,4), col.names=c("a","b","c"))) # ανάγνωση με συγκεκριμένο format fix.d # η λίστα με συνιστώσες a, b, c fix1.d <- read.fwf("fixwidth.dat", widths=c(1,3,4), col.names=c("a","b","c")) # ανάγνωση με συγκεκριμένο format fix1.d # πίνακας με στήλες a, b, c #----- Έξοδος δεδομένων – Εκτύπωση -------------------------- ??(airquality) # πληροφορίες για αρχείο δεδομένων data(airquality) # κατέβασμα του αρχείου attach(airquality) # ώστε οι στήλες του να είναι διαθέσιμες με το όνομά τους airquality[1:3,] # έλεγχος των 3 πρώτων γραμμών # Ozone Solar.R Wind Temp Month Day #1 41 190 7.4 67 5 1 #2 36 118 8.0 72 5 2 #3 12 149 12.6 74 5 3 sink("airquality.res") # ανοίγουμε αρχείο στον τρέχοντα κατάλογο # για αποθήκευση των αποτελεσμάτων cat("Το Αρχείο airquality \n", "από το data sets \n", " ---------------------------- ", "\n") # επικεφαλίδα για τα επόμενα στοιχεία airquality # εμφάνιση του αρχείου cat("\n \n summary(airquality) \n") # επικεφαλίδα summary(airquality) # περιγραφικά μέτρα θέσης plot(Temp,Ozone) # γραφική παράσταση sink() # κλείσιμο αρχείου sink("airquality.res ",append=T) # άνοιγμα υπάρχοντος αρχείου για προσθήκη cat("\n \n πίνακας διασπορών-συνδιασπορών του airquality \n") # επικεφαλίδα var(airquality, na.rm=T) # πίνακας διασπορών- συνδιασπορών χςρίς ελλείπουσες τιμές sink() # κλείσιμο αρχείου #............. dump και dput ........................... larger <- function(x, y) { y.is.bigger <- y > x x[y.is.bigger] <- y[y.is.bigger] x } # ορισμός συνάρτησης x<-c(12.3, 32, 15, 21.2, -2.3, 0) y<-c(22.5, 25, -24, 2.7, 1, 2.4) larger(x,y) # εφαρμογή της συνάρτησης dump("larger","larger.func") # μεταφορά σε αρχείο larg.df <- data.frame(x,y,larger(x,y)); larg.df dump(c("larger","larg.df"),"larger.func") # μεταφορά σε αρχείο και της εφαρμογής dump("state.x77","state1.dat") # μεταφορά δεδομένων σε αρχείο με την dump dput(state.x77,"state2.dat") # μεταφορά δεδομένων σε αρχείο με την dput state <- dget("state2.dat") state==state.x77 # έλεγχος όλα TRUE #............. write.table ..................... write.table(airquality, file="air.dat") # αποθηκεύεται πλαίσιο δεδομένων σε αρχείο read.table("air.dat") # έλεγχος # εξαγωγή αποτελεσμάτων write.table(coef(lm(Ozone~Solar.R+Wind+Temp)), file="myresults1.csv") read.table("myresults1.csv") # έλεγχος write.table(cor(na.omit(airquality)[,1:4]), file="myresults2.csv") read.table("myresults2.csv") # έλεγχος # ...... Πράξεις με Διανύσματα και Πίνακες ---------------------- #-----Οι συναρτήσεις crossprod και outer -------------------- salary.df=read.table("salary.df") attach(salary.df) X <- matrix(cbind(1,exper,educ,manag),ncol=4);X Y <- matrix(salary) B <- crossprod(X,Y) M <- crossprod(X);M x <- 1:5 y <- 6:1 outer(x,y) f <- function (x, y) exp(-x-y)* sin(y) z <- outer(x, y, f); z outer(month.name, 2001:2003, paste) p <- c(0.30, 0.15, 0.20, 0.25) ; p #[1] 0.30 0.15 0.20 0.25 VC <- -3*outer(p,p) ; VC # [,1] [,2] [,3] [,4] #[1,] -0.270 -0.1350 -0.18 -0.2250 #[2,] -0.135 -0.0675 -0.09 -0.1125 #[3,] -0.180 -0.0900 -0.12 -0.1500 #[4,] -0.225 -0.1125 -0.15 -0.1875 diag(VC) <- 3*p*(1-p) ; VC # [,1] [,2] [,3] [,4] #[1,] 0.630 -0.1350 -0.18 -0.2250 #[2,] -0.135 0.3825 -0.09 -0.1125 #[3,] -0.180 -0.0900 0.48 -0.1500 #[4,] -0.225 -0.1125 -0.15 0.5625 #------------ Οι συναρτήσεις apply και sweep --------------- # μέσες τιμές array ως προς συνδυασμό σταθμών του, εδώ των επιπέδων 1 και 3 meanΑar13 <- apply(Aarray, c(1,3),mean); meanΑar13 # [,1] [,2] [,3] #[1,] 4 14 114 #[2,] 5 15 115 apply(X, 2, mean, trim=.25) # αποκομμένες μέσες τιμές στηλών του Χ #[1] 1.0000000 6.5172414 2.0000000 0.4482759 library(AER) # άνοιγμα του πακέτου AER # αν του πακέτου AER δεν υπάρχει στον υπολογιστή σας τρέξτε την επόμενη εντολή # install.packages("AER", repos = "http://cran.r-project.org", # destdir = NULL, dependencies = TRUE) data("OECDGas") OECDGas # εμφάνιση των δεδομένων grGas <- OECDGas[OECDGas$country=="Greece",] # οι εγγραφές για την Ελλάδα grGasm <- as.matrix(grGas[-c(1:2)]);grGasm # μετατροπή σε πίνακα των 4 τελευταίων στηλών grm <- apply(grGasm,2,mean);grm # οι μέσες τιμές των στηλών grv <- apply(grGasm,2,var);grv # οι διασπορές των στηλών grdevs <- sweep(grGasm, 2, grm);grdevs # αφαίρεση από κάθε στήλη της μέσης τιμής της grstnd <- sweep(grdevs, 2, sqrt(grv), "/");grstnd # διαίρεση κάθε στήλης με την τυπική της απόκλιση round(apply(grstnd,2,mean),10) # έλεγχος μέσης τιμής "τυπικών τιμών" apply(grstnd,2,var) # έλεγχος διασποράς "τυπικών τιμών" # ------- Συναρτήσεις που εφαρμόζονται σε πίνακες ------------ set.seed(150) A <- matrix(10*round(runif(9),1),ncol=3);A # τυχαίος 3x3 πίνακας Α b <- matrix(10*round(runif(3),1),ncol=1);b # τυχαίος 3x1 πίνακας b solve(A) # ο αντίστροφος του Α solve(A,b) # το διάνυσμα-λύση του συστήματος Ax=b R=A;L=A R[2,1]=0;R[3,1]=0;R[3,2]=0;R # μετατροπή σε άνω τριγωνικό L[1,2]=0;L[1,3]=0;L[2,3]=0;L # μετατροπή σε κάτω τριγωνικό br<-backsolve(R,b) # το διάνυσμα-λύση του συστήματος Rx=b R %*% br # έλεγχος bl<-forwardsolve(L,b) # το διάνυσμα-λύση του συστήματος Lx=b L %*% bl # έλεγχος eigen(A, only.values=T)$values # ιδιοτιμές eigen(A) # ιδιοτιμές και ιδιοδιανύσματα det(A) # ορίζουσα exp(determinant(A)$modulus) # ορίζουσα prod(eigen(A)$values) # ορίζουσα sum(diag(A)) # ίχνος sum(eigen(A)$values) # ίχνος set.seed(200) # ορίζουμε ένα τυχαίο 3x3 πίνακα x <- matrix(sample(-3:3, size = 9, replace = T), nrow = 3, ncol = 3) ; x M <- x %*% t(x) ; M # ο συμμετρικός πίνακας Μ=x*x΄ U <- chol(M) # Choleski decomposition round(U, 3) t(U) %*% U # έλεγχος set.seed(200) x <- matrix(sample(-3:6,size = 12,replace = T), nrow = 3) ; x svdx <- svd(x) ; svdx # singular value decomposition x1 <- svdx$u %*% diag(svdx$d) %*% t(svdx$v) round(x1 - x, 3) # έλεγχος round(t(svdx$u) %*% svdx$u, 6) # έλεγχος ορθογωνιότητας u round(t(svdx$v) %*% svdx$v, 6) # έλεγχος ορθογωνιότητας v L <- qr(x) ; L # QR decomposition q <- qr.Q(L);q # ο πίνακας Q r <- qr.R(L);r # ο πίνακας R round(q%*%r,6) # έλεγχος ισότητας QR=x round(t(q)%*%q, 6) # έλεγχος ορθογωνιότητας Q # ------ Γραφήματα ----------------- plot(rnorm(10)) plot(1:10,runif(10)) x11() # ανοίγει μία συσκευή (η 2 που είναι ενεργός) plot(1:10) # γράφημα στη συσκευή 2 x11() # ανοίγει μία άλλη συσκευή (η 3 που είναι ενεργός, ενώ η 2 ανενεργός)) plot(rnorm(10)) # γράφημα στη συσκευή 2 dev.set(dev.prev()) # κάνουμε ενεργό τη συσκευή 2 και ανενεργό την 3 abline(0,1) # through the 1:10 points # προσθέτουμε την ευθεία y=x dev.set(dev.next()) # κάνουμε ενεργό τη συσκευή 3 και ανενεργό την 2 abline(h=0, col="gray") # προσθέτουμε την ευθεία y=0, με χρώμα γκρι dev.set(dev.prev()) # κάνουμε ενεργό τη συσκευή 2 και ανενεργό την 3 dev.off(); dev.off() # κλείνουμε και τις δύο συσκευές # ----- συσκευές Postscript --------- postscript(file="Plots.ps") # ανοίγει η συσκευή Postscript που # θα τυπώνει τα επόμενα στο αρχείο Plots.ps plot(1:10,runif(10)) # ένα γράφημα dev.off() # κλείνουμε τη συσκευή library(Cairo) # αν δεν υπάρχει τρέξτε την επόμενη εντολή #install.packages("Cairo", repos = "http://cran.r-project.org", # destdir = NULL, dependencies = TRUE)cairo_ps(file="RPlots.ps") # ανοίγει η συσκευή Postscript του πακέτου # Cairo που θα τυπώνει τα επόμενα στο αρχείο RPlots.ps cairo_ps(file="RPlots.ps") # ανοίγει η συσκευή Postscript του Cairo plot(rnorm(1000),rnorm(1000)) # ένα γράφημα dev.off() # κλείνουμε τη συσκευή cairo_pdf(file="Rplots.pdf") # ανοίγει η συσκευή Postscript του Cairo plot(rnorm(1000),rnorm(1000)) # ένα γράφημα dev.off() # κλείνουμε τη συσκευή # ......... Οι γάτες του Fisher .................... library(MASS) # άνοιγμα βιβλιοθήκης MASS attach(cats) # ενεργοποίηση αρχείου δεδομένων cats (help(cats)) catsF <- lm(Hwt~Bwt,data=cats, subset=Sex=="F") # γραμμική παλινδρόμηση για τις θηλυκιές γάτες catsM <- update(catsF, subset = Sex == "M") # γραμμική παλινδρόμηση για τις αρσενικιές γάτες plot(Bwt, Hwt, xlab="Body Weight (kg)", ylab="Heart Weight (gm)",type="n") # άδειο γράφημα με άξονες text(Bwt, Hwt, c("+","-")[Sex]) # πρόσθεση σημείων με "+" οι θηλυκιές γάτες, με "-" οι αρσενικιές legend(2.0,20,c("Females"," ","Males"), pch="+ -", bty="n") # προσθήκη κατάλληλης λεζάντας lines(Bwt[Sex == "F"], fitted(catsF), lty=1,col=8,lwd=2) # προσθήκη ευθείας παλινδρόμησης για τις θηλυκιές γάτες lines(Bwt[Sex == "M"], fitted(catsM), lty=4,col=9,lwd=2) # προσθήκη ευθείας παλινδρόμησης για τις αρσενικιές γάτες # ... ένα παράδειγμα υποδιαίρεσης οθόνης και καθορισμού κλίμακας ........................ set.seed(1000); ynorm <- rnorm(1000) # 1000 κανονικοί τυχαίοι αριθμοί par(mfrow=c(3,1)) # υποδιαίρεση της οθόνης (3 γραμμές, 1 στήλη) plot(1:1000, ynorm) # το γράφημα plot(1:1000, ynorm, log="x") # το γράφημα με λογαριθμημένο άξονα x plot(log10(1:1000), ynorm,tck=1) # το γράφημα με λογαριθμημένες τιμές x par(mfrow=c(1,1)) # επαναφορά της οθόνης στην κανονική της μορφή # -------- Γράφημα του αρχείου salary.df -------------- salary.df=read.table("salary.df") attach(salary.df) plot(exper, salary) # ... τροποποιημένο γράφημα ........... plot(exper,salary,type="n", xlab="εμπειρία (exper)", ylab= "μισθός (salary)") points(exper[educ==1],salary[educ==1],pch=16,col=3) points(exper[educ==2],salary[educ==2],pch=17,cex=1.3,col=4) points(exper[educ==3],salary[educ==3],pch=18,cex=1.6,col=5) title(main="Μισθός ανά εμπειρία (και εκπαίδευση)") legend(locator(1), legend=c("educ=1", "educ=2", "educ=3"), pch=16:18, col=c(3,4,5)) # ....... εναλλακτικά...... plot(exper,salary,type="n", xlab="εμπειρία (exper)", ylab= "μισθός (salary)") points(exper[educ==1],salary[educ==1],pch=16,col=3) points(exper[educ==2],salary[educ==2],pch=17,cex=1.3,col=4) points(exper[educ==3],salary[educ==3],pch=18,cex=1.6,col=5) title(main="Μισθός ανά εμπειρία (και εκπαίδευση)") legend(15,600, legend=c("educ=1", "educ=2", "educ=3"), col=c(3,4,5), pch=16:18) x1<-seq(-pi,0,length=500) x2<-seq(0,1,length=500) x3<-seq(1,4,length=500) y1<-cos(x1) y2<-exp(-x2) y3<-sqrt(x3) x<-c(x1,x2,x3) y<-c(y1,y2,y3) plot(x,y,type="n",axes=F,xlab="",ylab="",pch=".") points(x1,y1,type="l",col=2,lwd=2) points(x2,y2,type="l",col=3,lwd=2) points(x3,y3,type="l",col=4,lwd=2) axis(1,pos=0, at=c(-pi,-pi/2,0,1:4), labels=c("-pi","-pi/2","0","1","2","3","4")) axis(2,pos=0) #--------- εναλλακτικά ----------------- plot(x1, y1, type = "o", axes = F, xlab = "", ylab = "", pch = ".", xlim = c( - pi, 4), ylim = c(-1, 2), col = 2) par(new = T) plot(x2, y2, type = "o", axes = F, xlab = "", ylab = "", pch = ".", xlim = c( - pi, 4), ylim = c(-1, 2), col = 3) par(new = T) plot(x3, y3, type = "o", axes = F, xlab = "", ylab = "", pch = ".", xlim = c( - pi, 4), ylim = c(-1, 2), col = 4) axis(1,pos=0, at=c(-pi,-pi/2,0,1:4), labels=c("-pi","-pi/2","0","1","2","3","4")) axis(2, pos = 0) # ------ Γραφική Παράσταση κατανομής Γ(1.7) drawgamma <- function(shape){ x <- qgamma(seq(.001,.999,length=200),shape) plot(x,dgamma(x,shape), type="l",xlab="",ylab="",axes=F) axis(1,pos=0);axis(2,pos=0) } # ορισμός σχετικής συνάρτησης drawgamma(sh<-1.7) # εφαρμογή text(locator(1),paste("Γάμμα(",sh,")",sp="")) # λεζάντα # ------ Γραφική Παράσταση Κανονικής κατανομής N(mean, sd^2) drawnorm <- function(mean=0,sd=1){ x <- qnorm(seq(.001,.999,length=200),mean,sd) plot(x,dnorm(x,mean,sd), type="l",xlab="",ylab="",axes=F) axis(1,pos=0);axis(2,pos=0) } # ορισμός σχετικής συνάρτησης drawnorm(m<-3,s<-1) # εφαρμογή text(locator(1),paste("Normal N(",m,",",s^2,")",sp="")) # λεζάντα