Podstawy R 1) zmienne nazwy do 32 znaków, na początku litera lub kropka jeżeli nazwa zaczyna się kropką to drugim znakiem nie może byc cyfra 2) Typy obiektów: a) numeryczne (liczby) - trzy klasy: - całkowite +- 1.6E+9 - rzeczywiste 12-13 cyfr, +- 1.5E+-308 zapis +-mE+-c, m - mantysa (nie musi zawierać kropki) c - cecha - zespolone - pary liczb rzeczywistych, stale: z = Re + Imi np.: z = 7.5 + 3.14i b)Typ czynnikowy (wyliczeniowy, kategoryczny). wektory wartości na kilku poziomach (kategoriach) do tworzenia obiektów tego typu służy funkcja factor >oceny=factor(c("ndst","dost","db","db","dost","bdb")) > oceny [1] ndst dost db db dost bdb Levels: bdb db dost ndst > summary(oceny) bdb db dost ndst 1 2 2 1 > b) tekstowe (do 500 znaków): stałe: "tekst" c) logiczne stałe: TRUE FALSE d) wyrażenia (expression): teksty zawierające zapisane wyrażenia (wzory) 3) Struktury danych: a) wektory (podstawowy typ obiektów): uporządkowany zbiór wielkości tego samego typu) przykłady: c(1 ,3.5, 2.4, 7.8) <== wektor liczb c("ABC","ANNA","Zorro","Woda") <== tekstów c(TRUE,FALSE,TRUE,TRUE) <== wektor logiczny b) listy Uporządkowany zbiór obiektów róznego typu list("Jan","Kowalski",1991,"Kraków","TRUE") elementy list (i wektorów) moga być nazwane) list(imie="Jan",nazwisko="Kowalski",ur=1991,zam="Kraków",stud="TRUE") zamiast poj. wartości można wprowadzać wektory d1 = list(imie=c("Jan","Piotr"),nazwisko=c("Kowalski","Nowak"),ur=c(1991,1995), + zam=c("Kraków","Opole"),stud=c(TRUE,FALSE)) $imie [1] "Jan" "Piotr" $nazwisko [1] "Kowalski" "Nowak" $ur [1] 1991 1995 $zam [1] "Kraków" "Opole" $stud [1] TRUE FALSE > d1$stud <==== mozemy pytać o elementy listy [1] TRUE FALSE c) ramki danych zbiory nazwanych wektorów o tej samej liczbie elementów, data.frame(d1=1:5,d2=6:10,nazw=c("a","b","c","d","e")) df = data.frame(d1=1:5,d2=6:10,nazw=c("a","b","c","d","e")) > df d1 d2 nazw 1 1 6 a 2 2 7 b 3 3 8 c 4 4 9 d 5 5 10 e > df$nazw #<=== jak w przypadku list można wywołać kolumne ramki [1] a b c d e Levels: a b c d e jezeli nie nadamy nazw to R zrobi to "po swojemu" df = data.frame(1:5,d2=6:10,c("a","b","c","d","e")) > df X1.5 d2 c..a....b....c....d....e.. 1 1 6 a 2 2 7 b 3 3 8 c 4 4 9 d 5 5 10 e d) tablice - struktury dwu- lub wiecej wymiarowe tworzenie tablic: a) z istniejącego wektora, przez przypisanie wymiarów tab = 1:20 dim(tab)=c(5,4) - tworzy tablice o 5 wierszach i 4 kolumnach sama funkcja dim zwraca informacje o organizacji tablicy dim(tab) wymiary tablicy moga byc zmieniane, pod warunkiem, że iloczyn wymiarów pozostanie równy liczbie elementów dim(tab)=c(2,10) dim(tab)=c(4,5) dim(tab)=c(10,2) b) użycie funkcji matrix mm = matrix(data=1:24,nrow=3,ncol=8) mm = matrix(data=1:24,nrow=4) # liczba kolumn zostanie obliczona mm = matrix(data=1:24) # utworzy macierz kolumnową tablica jest wypełniana "kolumnami", jęzeli chcemy zmienić kolejność uzywamy parametru "byrow=TRUE" mm = matrix(data=1:24,nrow=3,ncol=8,byrow=TRUE) c) użycie funkcji array (tablica może mieć więcej niz dwa wymiary) ma = array(data=1:81,dim=c(3,3,9)). Wektory: Konstruktor: funkcja c(v1,v2,....,vn), vi - stała, zmienna, wektor składowe muszą bycć tego samego typu. wekt=c(1,2,5,6,7) > # elementy wektora moga miec nazwy > c(pierwszy = 12, drugi = 10, trzeci = 18) pierwszy drugi trzeci 12 10 18 sekwencje (ciągi): w1 = 1:10 w2 = 10:(-2) w3=seq(from=0,to=2*pi,by=0.01) <== zadany skok wartości w3=seq(0,2*pi,length.out=400) <== przedzial [0,2pi] dzielimy na 399 czesci w4=seq(length=51, from=-5, by=0.2) w5=seq(from=-5,by=0.1,along=w4) <== ma tyle samo elementów ile w4 powtórzenia: x = 1:3 w6=rep(x,times=5) w6 [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 w7=rep(x,each=5) w7 [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 0/0 NaN # <=== (Not a Number) brakujące dane: NA (Not Available) braki = c(1,3.5,2.7,NA,7.5,24) Operacje arytmetyczne i logiczne +, - , *, /, ^ Standardowe operatory arytmetyczne. %% Reszta modulo z dzielenia. %/% Dzielenie całkowite. wektory lub macierze powinny mieć ten sam rozmiar lub jeden z nich jest pojedynczą liczbą! ??? > x = 1:3 x*5 [1] 5 10 15 > y=6:10 > x*y [1] 6 14 24 9 20 <=== ????? jak to sie liczy Warning message: In x * y : longer object length is not a multiple of shorter object length %*% Iloczyn dwóch macierzy. > m1=matrix(data=1:6,ncol=3) > m1 [,1] [,2] [,3] [1,] 1 3 5 [2,] 2 4 6 > m2=matrix(data=1:9,nrow=3) > m2 [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > m1 %*% m2 [,1] [,2] [,3] [1,] 22 49 76 [2,] 28 64 100 > operatory logiczne <, ==, >, <=, >=, != Standardowe operatory porównywania wartosci liczbowych (relacje). ! Operator negacji. &, &&, |, || Logiczny iloczyn oraz logiczna suma. any(), all() Logiczna suma(iloczyn) wszystkich elementów wektora. > 2 > 3 [1] FALSE > > !2 > 3 <=== negacja [1] TRUE > > x = 1:10 > x >= 5 [1] FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE > x >= 5 & x > 7 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE > > x >= 5 | x < 3 [1] TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE > > x >= 5 && x < 7 # iloczyn i suma logiczna dla wszystkich elementów [1] FALSE # obliczane łącznie > x >= 5 || x < 7 [1] TRUE > any(x>5) [1] TRUE > all(x>5) [1] FALSE > UWAGA: Jeżeli zamiast zmiennej czy stałej logicznej wstawimy do wyrażenia logicznego liczbe (stałą lub zmienną) to będzia traktowana jako FALSE gdy ma wartość zero a jako TRUE dla wartści różnych od zera !!!!! Indeksowanie >x = 1:10 > x[x==6] [1] 6 > x[x>=7] [1] 7 8 9 10 > x[-3] [1] 1 2 4 5 6 7 8 9 10 > y=c(1,2,3,NA,5,6) > z=y[!is.na(y)] <=== kopiujemy jedynie zdefiniowane elementy wektora > z [1] 1 2 3 5 6 > indeksowanie macierzy > x = array(data=1:20,dim=c(4,5)) > x [,1] [,2] [,3] [,4] [,5] [1,] 1 5 9 13 17 [2,] 2 6 10 14 18 [3,] 3 7 11 15 19 [4,] 4 8 12 16 20 > > > x[3,4] [1] 15 > x[2,] [1] 2 6 10 14 18 > x[,3] [1] 9 10 11 12 > x[1:3,1:3] [,1] [,2] [,3] [1,] 1 5 9 [2,] 2 6 10 [3,] 3 7 11 > x[,2:4] # wybieramy kolumny 2:4 z macierzy x [,1] [,2] [,3] [1,] 5 9 13 [2,] 6 10 14 [3,] 7 11 15 [4,] 8 12 16 > macierze indeksów > i = array(data=c(1:3,3:1),dim=c(3,2)) > i [,1] [,2] [1,] 1 3 [2,] 2 2 [3,] 3 1 > x[i] [1] 9 6 3 > x[i] = 0 > x [,1] [,2] [,3] [,4] [,5] [1,] 1 5 0 13 17 [2,] 2 0 10 14 18 [3,] 0 7 11 15 19 [4,] 4 8 12 16 20 > a=outer(1:10,1:10) > a [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 2 3 4 5 6 7 8 9 10 [2,] 2 4 6 8 10 12 14 16 18 20 [3,] 3 6 9 12 15 18 21 24 27 30 [4,] 4 8 12 16 20 24 28 32 36 40 [5,] 5 10 15 20 25 30 35 40 45 50 [6,] 6 12 18 24 30 36 42 48 54 60 [7,] 7 14 21 28 35 42 49 56 63 70 [8,] 8 16 24 32 40 48 56 64 72 80 [9,] 9 18 27 36 45 54 63 72 81 90 [10,] 10 20 30 40 50 60 70 80 90 100 d = outer(1:10,1:10,"/") <=== można uzyć dowolnego operatora > d [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 0.5 0.3333333 0.25 0.2 0.1666667 0.1428571 0.125 0.1111111 0.1 [2,] 2 1.0 0.6666667 0.50 0.4 0.3333333 0.2857143 0.250 0.2222222 0.2 [3,] 3 1.5 1.0000000 0.75 0.6 0.5000000 0.4285714 0.375 0.3333333 0.3 [4,] 4 2.0 1.3333333 1.00 0.8 0.6666667 0.5714286 0.500 0.4444444 0.4 [5,] 5 2.5 1.6666667 1.25 1.0 0.8333333 0.7142857 0.625 0.5555556 0.5 [6,] 6 3.0 2.0000000 1.50 1.2 1.0000000 0.8571429 0.750 0.6666667 0.6 [7,] 7 3.5 2.3333333 1.75 1.4 1.1666667 1.0000000 0.875 0.7777778 0.7 [8,] 8 4.0 2.6666667 2.00 1.6 1.3333333 1.1428571 1.000 0.8888889 0.8 [9,] 9 4.5 3.0000000 2.25 1.8 1.5000000 1.2857143 1.125 1.0000000 0.9 [10,] 10 5.0 3.3333333 2.50 2.0 1.6666667 1.4285714 1.250 1.1111111 1.0 > macierze specjalne > v = c(1,2,3,4,5) > diag(v) <=== funkcja diag, argument to wektor [,1] [,2] [,3] [,4] [,5] [1,] 1 0 0 0 0 [2,] 0 2 0 0 0 [3,] 0 0 3 0 0 [4,] 0 0 0 4 0 [5,] 0 0 0 0 5 > m=array(data=1:36,dim=c(6,6)) > m [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 7 13 19 25 31 [2,] 2 8 14 20 26 32 [3,] 3 9 15 21 27 33 [4,] 4 10 16 22 28 34 [5,] 5 11 17 23 29 35 [6,] 6 12 18 24 30 36 > diag(m) <=== argumentem diag jest macierz [1] 1 8 15 22 29 36 > Rownania liniowe A x = b A - macierz kwadratowa, x wektor(kol), b= wektor(kol) x - wektor rozwiązań > A [,1] [,2] [,3] [,4] [,5] [1,] 1.0 6.00 2.5 16.00 21.00 [2,] 2.0 7.00 12.0 17.00 11.00 [3,] 3.0 8.00 13.0 18.00 23.00 [4,] 4.0 9.00 14.0 19.00 24.00 [5,] 13.3 7.37 2.5 4.22 8.37 > b [1,] 248.9200 [2,] 219.3700 [3,] 318.3100 [4,] 337.2800 [5,] 122.0009 >x = solve(A,b) # <== rozwiazalismy układ równań liniowych >x # macierz A to współczynniki, b = wyrazy wolne x # x - wektor rozwiązań [1,] 1.20 [2,] 2.30 [3,] 3.70 [4,] 4.50 [5,] 7.27 # funkcja solve moze byc wywolana tylko z macierza wspólczynników # bez wektora wyrazów wolnych - jako wynik daje macierz odwrotną do A # przykład: > m = matrix(data=rnorm(100, 0, sd=2),ncol=10) # m - macierz liczb losowych > m [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1.67498430 -1.5170697 0.30019641 -3.5189839 2.2784411 -0.4594726 [2,] 0.24030595 -0.8180661 0.07293386 -1.3187860 6.0678052 -0.8293906 [3,] -1.11005852 2.2719240 4.70903775 -1.7290213 4.0782852 2.5049563 [4,] 2.93213341 1.6039369 -2.29256001 -2.0532593 -0.6077552 2.6265419 [5,] -3.30305162 -1.1887949 0.78310686 -2.2550469 -0.3167670 -1.0849581 [6,] 0.15185438 0.6774267 -1.55520657 1.7834675 -1.6478389 2.9086296 [7,] -0.80882763 1.4700707 -3.05373079 -2.6725921 1.8700990 2.9815190 [8,] -0.93559044 -0.8837732 0.69308657 -0.6584932 -1.3380378 0.2220822 [9,] 0.09792175 2.9666755 -2.30412833 1.5342386 2.0285016 0.1196257 [10,] 0.12474493 0.1470531 -0.64126313 -2.0266899 1.1951560 -1.5479236 [,7] [,8] [,9] [,10] [1,] 0.1505598 0.0986578 -0.7340547 2.0014618 [2,] -0.8220253 0.7045624 1.6821123 -2.6861366 [3,] 0.8063396 1.6133942 -4.9330954 -0.4192431 [4,] 2.3571510 -5.2713044 -1.5165878 -3.6336524 [5,] 0.1609639 3.1814710 2.1746334 -0.8024235 [6,] 2.0798974 2.5115399 -2.0973668 2.3223904 [7,] -2.7353438 -0.3138051 -4.2432918 1.3157193 [8,] -0.9525206 -0.3938316 -1.8213786 -3.1362713 [9,] -0.4009397 -0.2689624 -2.2088725 1.1064806 [10,] 1.1955689 -0.3615500 2.7510395 3.7773370 > mo = solve(m) # <=== mo to macierz odwrotna do m > mo [,1] [,2] [,3] [,4] [,5] [1,] 0.480581424 -0.25459276 -0.0630699211 0.14630438 0.32944433 [2,] 0.205183916 -0.32698988 0.0377465896 0.17974327 0.47230556 [3,] -0.005780227 -0.08488879 0.1085879300 -0.01746754 -0.01299851 [4,] -0.234627568 0.16642747 -0.0007795198 -0.11975649 -0.29724368 [5,] -0.145028500 0.23491355 0.0386968037 -0.06898931 -0.21548312 [6,] -0.223199748 0.08980998 0.0483446512 0.02132207 -0.10454026 [7,] -0.112819657 0.15797866 0.0430218569 0.01573817 -0.14865790 [8,] 0.351870346 -0.18044240 -0.0416880966 0.07181509 0.43706614 [9,] -0.066875523 -0.02697441 -0.0191867020 0.06286077 0.13761729 [10,] -0.078211681 0.01087005 0.0225637928 -0.09237168 -0.17695792 [,6] [,7] [,8] [,9] [,10] [1,] -0.29996289 -0.0573407934 -1.0105590 0.11728786 -0.90097687 [2,] -0.39625202 -0.0008404219 -0.9998334 0.24382007 -0.72147175 [3,] -0.12931457 -0.0172971252 -0.2036870 -0.12228541 -0.11258319 [4,] 0.18516853 -0.0399159112 0.3985484 -0.04694553 0.30895687 [5,] 0.14758581 0.0066297187 0.3609360 -0.03502608 0.35294294 [6,] 0.09272156 0.1780128980 -0.1689993 -0.34084635 0.02631187 [7,] 0.31347677 -0.1584823249 0.6252794 0.05195133 0.52686827 [8,] -0.15026811 -0.0439189866 -0.7500500 0.16681074 -0.72138895 [9,] -0.13909376 0.0613772491 -0.5122975 -0.14641715 -0.21449819 [10,] 0.08172985 0.0496806590 0.1715072 -0.08174829 0.28875519 > wyn = m %*% mo # <== sprawdzamy czy iloczyn m * mo = 1 > wyn [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 9.367507e-17 -1.387779e-17 0.000000e+00 3.885781e-16 [2,] -1.110223e-16 1.000000e+00 4.163336e-17 -2.775558e-17 0.000000e+00 [3,] -5.551115e-17 1.734723e-18 1.000000e+00 6.245005e-17 1.665335e-16 [4,] -2.220446e-16 4.371503e-16 5.551115e-17 1.000000e+00 1.110223e-16 [5,] -4.302114e-16 5.898060e-17 3.122502e-17 1.387779e-17 1.000000e+00 [6,] 8.326673e-17 -2.775558e-17 -4.163336e-17 8.326673e-17 5.551115e-17 [7,] 8.326673e-17 -3.469447e-18 2.428613e-17 -9.714451e-17 2.220446e-16 [8,] -1.110223e-16 6.938894e-18 0.000000e+00 -1.110223e-16 -1.110223e-16 [9,] 5.551115e-17 2.550044e-16 -4.857226e-17 1.526557e-16 1.110223e-16 [10,] -1.110223e-16 1.110223e-16 8.326673e-17 0.000000e+00 -2.220446e-16 [,6] [,7] [,8] [,9] [,10] [1,] -3.608225e-16 1.387779e-17 0.000000e+00 -3.608225e-16 3.330669e-16 [2,] -1.110223e-16 0.000000e+00 4.996004e-16 0.000000e+00 -2.220446e-16 [3,] -3.330669e-16 -1.318390e-16 -1.665335e-16 -7.632783e-17 -2.081668e-16 [4,] 1.665335e-16 -8.326673e-17 2.220446e-16 5.551115e-17 2.220446e-16 [5,] 1.110223e-16 -1.387779e-17 -8.326673e-17 -1.665335e-16 -7.494005e-16 [6,] 1.000000e+00 -1.387779e-17 1.110223e-16 3.885781e-16 -1.110223e-16 [7,] -2.220446e-16 1.000000e+00 -8.326673e-17 3.330669e-16 -6.661338e-16 [8,] -5.551115e-17 -2.775558e-17 1.000000e+00 -5.551115e-17 2.220446e-16 [9,] -2.775558e-17 -6.938894e-18 3.053113e-16 1.000000e+00 2.775558e-16 [10,] 5.551115e-17 0.000000e+00 2.220446e-16 -3.885781e-16 1.000000e+00 > # poza diagonalą macierzy wyn są bardzo małe liczby, # które są wynikiem dokładności obliczeń (12-13 cyfr > wyn[wyn < 1e-10] = 0 # <== wyrzucamy z wyn wszystkie elementy mniejsze niż > wyn # <== 10^-10 i sprawdzamy... [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 0 0 0 0 0 0 0 0 0 [2,] 0 1 0 0 0 0 0 0 0 0 [3,] 0 0 1 0 0 0 0 0 0 0 [4,] 0 0 0 1 0 0 0 0 0 0 [5,] 0 0 0 0 1 0 0 0 0 0 [6,] 0 0 0 0 0 1 0 0 0 0 [7,] 0 0 0 0 0 0 1 0 0 0 [8,] 0 0 0 0 0 0 0 1 0 0 [9,] 0 0 0 0 0 0 0 0 1 0 [10,] 0 0 0 0 0 0 0 0 0 1 > # rzeczywiście otrzymaliśmy macierz jednostkową Wprowadzanie danych - część I Elementarne wprowadzanie danych i korekty 1) wpisywanie bezpośrednio do wektorów v = c(1,2, 2.5, 3.14) # z użyciem funkcji c(.... ) 2) z użyciem funkcji scan( ) # czyta dane z klawiatury koniec danych <== wciśnięcie Enter w pustej linii > v1 = scan() 1: 1.3 2: 2.4 3: 5.6 4: <=== tutaj naciśnięto Enter Read 3 items > v1 [1] 1.3 2.4 5.6 3) z użyciem funkcji scan(nazwa pliku) przykład: plik danev2.txt w kartotece "dane" na dysku d: (Windows) zawiera tekst: 1.5 4.6 7.5 11.2 3.5 2.7 1.5 czytany przez scan: > v2 = scan("d:/dane/danev2.txt") # <=== Read 7 items > v2 # sprawdzamy co przeczytano: [1] 1.5 4.6 7.5 11.2 3.5 2.7 1.5 Jeżeli pomiędzy liczbami występują np. przecinki czy inne znaki niż odstępy to należy poinformować o tym funkcję scan plik danev3.txt zawiera tekst: 1.5, 4.6, 7.5, 11.2, 3.5, 2.7, 1.5, czytany przez scan: v3 = scan("d:/dane/danev3.txt", sep=",") # <== sep oznacza "separator" > v3=scan("d:/danev3.txt",sep=",") Read 10 items # <=== ma więcej danych????? skąd, # dodaliśmy jedynie przecinki! > v3 [1] 1.5 4.6 7.5 NA 11.2 3.5 2.7 NA 1.5 NA # są trzy NA: ^ ^ ^ # Wniosek: po ostatniej liczbie w linijce nie należy stawiać # znaku separatora, te rolę pełni wtedy znak nowej linii funkcja scan może czytać rozmaite dane: rządzi tym parametr "what" xl = scan("nazwapliku",what="logical",nmax=20) # czyta 20 wartości logicznych x2 = scan("nazwapliku",what="numeric",nmax=20) # czyta 20 liczb (domyslne) x3 = scan("nazwapliku",what="integer",nmax=20) # czyta 20 liczb całkowitych x4 = scan("nazwapliku",what="complex",nmax=20) # czyta 20 liczb zespolonych x5 = scan("nazwapliku",what="character",nmax=20) # czyta 20 tekstów l1 = scan("nazwapliku",what=c("","","")) # wczytujemy listę w kazdej linii wczytujemy trzy elenty listy przyklad: plik lista1.txt zawiera tekst: Jan Kowalski 340 Piotr Nowak 730 Paweł Kot 600 wczytany rozkazem: lis1 = scan("d:/lista1.txt",what=list("","","")) Read 3 records > lis1 [[1]] [1] "Jan" "Piotr" "Paweł" [[2]] [1] "Kowalski" "Nowak" "Kot" [[3]] [1] "340" "730" "600" >names(lis1) # elenty listy nie mają nazw! NULL >names(lis1)=c("imie","nazwisko","wplata") >lis1 # <== sprawdzamy czy się coś zmieniło $imie [1] "Jan" "Piotr" "Paweł" $nazwisko [1] "Kowalski" "Nowak" "Kot" $wplata [1] "340" "730" "600" # Czy wpłaty to liczby??? >lis1$wplata[1]+lis1$wplata[2] Error in lis1$wplata[1] + lis1$wplata[2] : # scan wczytala wpłaty jako non-numeric argument to binary operator # teksty >#Mozemy je zamienic na liczby jedn a z funkcji konwersji typu lis1$wplata=as.numeric(lis1$wplata) > lis1 $imie [1] "Jan" "Piotr" "Paweł" $nazwisko [1] "Kowalski" "Nowak" "Kot" $wplata [1] 340 730 600 ========================================== drobne korekty wartości danych: > data.entry(x) # tylko, jeżeli x już wcześniej zdefiniowano > v=c(NA) data.entry(v) # <== w okienku data entry wpisano w kolumnie v # kolejne liczby 1 : 4 > v # sprawdzamy zawartośc v [1] 1 2 3 4 > data.entry(x=c(NA)) # jeżeli tworzymy x # w oknie funkcji data.entry należy # w nagłówku pierwszej kolumny wpisać x # a poniżej liczby lub teksty (wybór zaznaczamy # w otwierającym sie okienku dilaogowym) > #ta sama funkcja (data.entry) występuje pod nazwą "de" # zmiany wartości ( i rozmiaru) wektora można dokonać z uzyciem edytora x = edit(x) # nie wolno użyć samego "edit(x)" ! x # można także uzyć komendy fix(x) # automatycznie zapamiętuje poprawki (nie da sie wycofać zmian) ========================================================================= Czytanie tablic danych (ramek danych) użycie: xx = read.table("nazwapliku",header=TRUE) # read.table() czyta tylko liczby i dane kategoryzowane !!! przykład: plik("d:/dane/dane.txt) zawiera tekst: pomiar temp cisnienie zw1 105 24.5 zw2 110 11.6 zw3 98 57.65 zw4 152 34.7 > wyniki = read.table("d:/dane/dane.txt",header=TRUE) > wyniki pomiar temp cisnienie 1 zw1 105 24.50 2 zw2 110 11.60 3 zw3 98 57.65 4 zw4 152 34.70 > wyniki$pomiar # < [1] zw1 zw2 zw3 zw4 Levels: zw1 zw2 zw3 zw4 # read.table może również użyć dodatkowych parametrów: nrow = liczba linii, które nalezy przeczytać sep="znak separatora" , np. sep=";" dec="znak ułamka dziesiętnego", np. dec="," skip=liczba linii na poczatku pliku, linii które należy opuścic przed czytaniem danych podobną rolę pełni funkcja: df = read.csv("nazwa pliku",header=TRUE,...) # czyta dane oddzielane przecinkami, # parametry jak read.table pliki .csv można wygenerować np. w Excelu funkcje: read.xls write.xls pozwalaja, odpowiednio, czytać i pisać pliki Excela. do uzycia potrzebny jest pakiet xlsReadWrite, który należy sprowdzić z Internetu. O B L I C Z E N I A na wektorach (i tablicach) ======================================================================== Funkcje statystyki opisowej, działające na całych wektorach : wartości skrajne max(v) - maksymalna wartość składowej wektora min(v) - minimalna wartość składowej wektora średnie mean(v)- wartość średnia, dodatkowe parametry: trim=0-50, oblicza średnią po odrzuceniu 200%*trim skrajnych wartości z v, na.rm=FALSE, czy wyrzucać wartości nieznane? weighted.mean(v,w) - średnia ważona, w - wektor wag geometric.mean(psych) - z pakietu "psych" harmonic.mean(psych) - z pakietu "psych" median(v) - wartość środkowa quantile(x) - wyznacza kwantyle ------------- przyklady > quantile(x <- rnorm(1001)) # Extremes & Quartiles by default 0% 25% 50% 75% 100% -2.960752743 -0.647673124 0.004523922 0.632159008 3.100949014 > quantile(x, probs = c(0.1, 0.5, 1, 2, 5, 10, 50, NA)/100) 0.1% 0.5% 1% 2% 5% 10% -2.920614903 -2.551016732 -2.206166415 -1.908824109 -1.605590363 -1.221637594 50% 0.004523922 NA -------------------------------------------------------------- IQR(v) - rozstęp międzykwartylowy range(x) zakres zmienności x [1] -2.960753 3.100949 var(x) - wariancja w próbie sd(v) - odchylenie standardowe cor(v), cov(v) macierz korelacji i kowariancji sd(v)/mean(v) - tzw współczynnik zmienności (CV - Coefficient of Variance) length(v) - liczba elementów wektora kurtosis(v) - (wym. pakiet e1071), Kurtoza, miara spłaszczenia rozkładu skewness(v) - (wym. pakiet e1071), Skośność, miara asymetrii rozkładu moda(v) - pakiet dprep, Dominanta, wartość występująca najczęściej użycie: ponizsze dane pochodzą ze strony autora podręcznika Przewodnik po pakiecie R, aut. P. Biecek daneSoc = read.table("d:/daneSoc",header=TRUE,sep=";") head(daneSoc) # <== oglądamy poczatek ramki wiek wyksztalcenie st.cywilny plec praca 1 70 zawodowe w zwiazku mezczyzna uczen lub pracuje 2 66 zawodowe w zwiazku kobieta uczen lub pracuje 3 71 zawodowe singiel kobieta uczen lub pracuje 4 57 srednie w zwiazku mezczyzna uczen lub pracuje 5 45 srednie w zwiazku kobieta uczen lub pracuje 6 48 srednie w zwiazku mezczyzna nie pracuje cisnienie.skurczowe cisnienie.rozkurczowe 1 143 83 2 123 80 3 167 80 4 150 87 5 130 83 6 138 75 attach(daneSoc) #<== nie chcemy każdorazowo wypisywać nazwy ramki # np. daneSoc$wiek, komenda attach włacza składniki # ramki do naszej przestrzeni roboczej i możemy traktować # nazwy kolumn jako nazwy zmiennych lokalnych badamy dane nt. wieku > range(wiek) # przedział wartości [1] 22 75 > IQR(wiek) # rozstęp międzykwartylowy [1] 23 > mean(wiek) # średni wiek [1] 43.16176 > mean(wiek,trim=0.2) # średnia po odrzuceniu 40% skrajnych wartości [1] 42.58065 > geometric.mean(wiek) # chcemy wyznaczyć średnią geometryczną Error: could not find function "geometric.mean" # nie znamy tej funkcji > library(psych) # ładujemy pakiet "psych" > geometric.mean(wiek) # i już umiemy to liczyć [1] 40.89757 > harmonic.mean(wiek) # przy okazji z pakietu psych załadowano [1] 38.66 # funkcję, która potrafi obliczyć średnią # harmoniczną > median(wiek) # mediana [1] 45 > moda(wiek) # jaki wiek występuje najliczniej w próbie? Error: could not find function "moda" > library(dprep) # potrzebny będzie pakiet dprep Loading required package: MASS #ładując się wprowadza dodatkowe Loading required package: nnet # biblioteki Loading required package: lattice Loading required package: class > moda(wiek) # <===== tu już potrafimy wyznaczyć wartość moda() [1] 26 > sd(wiek) # odchylenie standardowe [1] 13.8471 > kurtosis(wiek) # czy rozkład jest skoncentrowany? Error: could not find function "kurtosis" > library(e1071) # znowu trzeba załadować bibliotekę > kurtosis(wiek) # wynik odległy od zera (kurtoza = 0 dla rozkładu normalnego) [1] -0.955848 > > skewness(wiek) # czy rozkład jest skośny? [1] 0.233151 Podsumowania: summary(wektor) > summary(wyksztalcenie) # <== dane kategoryczne podstawowe srednie wyzsze zawodowe 93 55 34 22 summary(wiek) # <== dane liczbowe Min. 1st Qu. Median Mean 3rd Qu. Max. 22.00 30.00 45.00 43.16 53.00 75.00 > Tablice wystąpień (kontyngencji) funkcja table > table(wyksztalcenie,praca) praca wyksztalcenie nie pracuje uczen lub pracuje podstawowe 22 71 srednie 16 39 wyzsze 6 28 zawodowe 8 14 > > table(wyksztalcenie,wiek,praca) , , praca = nie pracuje wiek wyksztalcenie 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 podstawowe 0 0 1 0 3 1 0 1 0 1 1 0 0 0 0 0 0 1 0 0 srednie 0 0 0 0 0 1 0 0 1 0 0 1 1 1 0 1 1 0 0 0 wyzsze 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 zawodowe 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 wiek wyksztalcenie 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 podstawowe 1 0 0 0 0 1 1 2 1 0 0 1 1 0 2 1 1 0 0 0 srednie 0 0 0 0 0 1 1 3 1 0 0 0 0 0 0 0 2 1 0 0 wyzsze 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 zawodowe 0 0 0 0 0 0 0 0 1 0 0 0 2 0 0 0 2 0 0 0 wiek wyksztalcenie 63 64 66 67 68 70 71 72 74 75 podstawowe 0 0 0 1 0 0 0 0 0 0 srednie 0 0 0 0 0 0 0 0 0 0 wyzsze 0 0 0 0 0 0 0 0 0 0 zawodowe 0 0 0 0 0 0 0 0 0 0 , , praca = uczen lub pracuje wiek wyksztalcenie 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 podstawowe 3 6 2 3 7 3 3 2 2 1 2 1 4 0 0 1 0 0 1 1 srednie 0 0 0 0 0 0 1 1 4 1 2 1 0 1 2 1 1 0 0 0 wyzsze 0 0 0 0 1 0 0 0 1 0 0 2 1 0 0 0 1 1 2 0 zawodowe 1 0 1 0 0 0 0 0 1 0 2 0 0 0 0 0 0 0 0 0 wiek wyksztalcenie 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 podstawowe 1 0 1 1 1 1 5 1 2 3 0 2 2 3 1 0 0 0 0 0 srednie 0 1 0 4 2 0 1 2 1 2 2 2 0 0 1 1 0 0 0 0 wyzsze 1 0 0 1 0 2 0 0 0 1 0 4 1 0 1 1 0 1 0 1 zawodowe 0 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 wiek wyksztalcenie 63 64 66 67 68 70 71 72 74 75 podstawowe 1 1 0 1 1 0 1 0 0 0 srednie 3 1 0 0 0 0 0 1 0 0 wyzsze 0 0 0 1 0 0 2 0 1 1 zawodowe 0 0 2 0 0 1 2 0 0 0 > # UWAGA - tablice wystąpięń (kontyngencji) dla wielu zmiennych > #o wielu różnych wartościach potrafią przyjąć ogromne rozmiary wykład 3 próby losowe > sample(1:49,6) [1] 3 19 33 15 34 13 <== mozemy wyręczyć totalizator > sample(1:49,6) [1] 47 25 5 18 49 14 >sample(1:49,6) > dziesiec losowych liter > sample(letters,10,T) [1] "u" "q" "x" "s" "q" "f" "c" "f" "l" "x" > # wylosujmy wektor cyfr od 1 do 3, z zadanymi prawdopodobienstwami wylosowania > sample(1:3,20,replace=TRUE, prob=c(0.6,0.3,0.1)) [1] 2 3 1 1 3 3 1 1 1 1 1 2 1 1 1 2 1 1 2 2 combn(n,k) - kombinacje n nad k > combn(5,3) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 1 1 1 1 1 2 2 2 3 [2,] 2 2 2 3 3 4 3 3 4 4 [3,] 3 4 5 4 5 5 4 5 5 5 >choose(n,k) - liczba tych kombinacji > choose(5,3) [1] 10 > > Liczby losowe 1) funkcja rnorm > rnorm(n,mean,sd) # <== losuje n liczb z rozkładem normalnym, # wokół wartości średniej mean # i średnim odchyleniu standardowym sd > r=rnorm(5) # losujemy 5 liczb, domyślnie mean = 0, sd = 1 > r [1] -0.60301666 0.01874362 0.71110645 -2.32731271 0.30780569 > mean(r) [1] -0.3785347 > sd(r) [1] 1.189994 # dla 5 liczb wyniki są kiepskie > r=rnorm(5000) > mean(r) [1] -0.01162106 > sd(r) [1] 0.9939866 # dla 5000 widać poprawę > r=rnorm(50000) > sd(r) [1] 1.003997 > mean(r) [1] 0.00353229 # dla 50000 jest wyrażnie lepiej (prawo wielkich liczb) > r=rnorm(500000) # dla pół miliona liczb losowych > mean(r) # widać dalszą poprawę wyników [1] 0.0004333034 # prawie o rząd wielkości > sd(r) [1] 1.000585 > > # po zdefiniowaniu wektora r o rozmiarze 500 000 elementów > # lepiej posprzątać pamięć... rm(r) # <== funkcja rm usuwa obiekt z pamięci 2) funkcja runif - wyznacza liczby o jednakowym prawdopodobieństwie wystąpienia w zadanym przedziale wywołanie: runif(n,min,max), domyslnie min=0, max=1 > runif(5) [1] 0.9111863 0.3747355 0.4957358 0.1021286 0.2929222 # 5 liczb losowych [0,1] runif(1,0,2) # <== ile poczekamy na zmiane swiatla na przejsciu [0,2] min [1] 1.615157 > runif(5,0,2) # <== a na pięciu przejściach ? [1] 1.596893 1.858550 1.911551 1.215455 0.266690 R Y S O W A N I E I. Rozkazy "wysokiego poziomu" Każde uzycie otwiera nowy rysunek (okno graficzne) hist, barplot, pie, boxplot, … (1-D plot) <== wykresy jednowymiarowe plot, curve, matplot, pairs,… ((x-y)plots) <== wykresy zalezności x,y image, contour, filled.contour,... (2-D surface plots) <== powierzchnie 2 wym. persp, scatterplot3d(4) (3-D plots) <== rysunki 3-wymiarowe II. Rozkazy "niższego poziomu" dodają nowe elementy do rysunku już wykonanego przez rozkaz z grupy "wysokiego poziomu" lines, points, segments, polygon, rect, text, arrows, legend, abline, locator, rug, itd. III. Parametry definiujące: a) wygląd obiektów na rysunku cex (rozmiar tekstu i symboli), col (kolor), font (czcionka), las (orientacja etykiet na osiach), lty (rodzaj linii), lwd (szerokość linii), pch (rodzaj znaków reprezentujących punkty), asp (stosunek długości odcinków na osiach x i y), itd. b) wygląd okna graficznego mar (wielkość marginesów), mfrow (liczba obrazów w wierszu), mfcol(liczba obrazów w kolumnie) przykład 1: x=1:100 y=runif(100) plot(x,y) # rysujemy wykres x-y z=rnorm(100) # generujemy 100 liczb o rozkładzie normalnym points(x,z,col="red",pch=3) # i dorysowujemy je do istniejącego rysunku przykład 2: data(Orange) head(Orange) #rysujemy zależność obwodu pnia od wieku drzewka plot(Orange$age, Orange$circumference,xlab="wiek, dni", ylab="obwód pnia, mm", main= "Wzrost drzewek pomarańczowych") #ponownie rysujemy, rozróżniając dane dla kolejnych pięciu drzewek plot(Orange$age, Orange$circumference,xlab="wiek, dni", ylab="circumference, mm", main= "Orange tree growth", pch=(15:20)[Orange$Tree], col=(1:5)[Orange$Tree], cex=1.3) legend("bottomright",pch=15:20,col=1:5,legend=1:5) # warto dodać legendę Rysunki funkcji matematycznych 1) funkcja curve(wyrażenie) np. curve(sin(3*pi*x)) # zakres x [0,1], na podstawie 101 wartości curve(sqrt(x)) # j.w. curve(sqrt) # mozna podać tylko nazwe funkcji # bardziej skoplikowany przypadek curve(sin(3*pi*x),from=0,to=2,col="blue", xlab="x",ylab="f(x)",main=" funkcja curve") curve(cos(3*pi*x),add=TRUE,col="red",lty=2) # dorysujemy jeszcze cosinus # add=TRUE <== dodaj do istniejącego # rysunku abline(h=0,lty=2) # dorysujemy poziomo na wysokosci 0 (h=0) # linie przerywaną (lty=2) # i dodamy opis legend("bottomleft",legend=c("sin","cos"),text.col=c("blue","red"),lty=1:2) # Kilka rysunków obok siebie funkcja par zmienia sposób rysowania par(mfrow=c(3,2)) # rysunki w trzech wierszach i dwu kolumnach, # rysowane "wierszami par(mfcol=c(3,2)) # rysunki w trzech wierszach i dwu kolumnach, # rysowane "kolumnami przykład par(mfrow=c(2,2)) for ( i in 1:4) curve(sin(i*pi*x),0,1,main=i) #Programowanie #Składnia instrukcji warunkowych if (war) instrukcja # to moze byc tzw. instrukcja złożona: # ciąg instrukcji ujętych w nawiasy # klamrowe {.... } if (war) # tutaj trzeba pilnować by "else" instrukcja1 else instrukcja2 # else nie rozpoczynalo linijki ifelse(war, instrukcja1,instrukcja2) # alternatywna wersja poprzedniej # instrukcji, stosowana zwykle # w przypadku gdy instrukcje 1 i 2 # są krótkie # Przykłady ################################################################## # Instrukcje warunkowe - przykłady ################################################################## Wyznaczanie pierwiastków trójmianu kwadratowego # równanie: a*x^2 +b*x + c = 0 a = b = c = delta = b^2 - 4*a*c # obliczamy wartośc delta if (delta != 0.0) { #jeżeli delta różni sie od zera to #mamy dwa różne pierwiastki delta = as.complex(delta) # zamieniamy wielkość delta na liczbe # zespoloną x1 = (-b + sqrt(delta))/(2*a) # pierwiastek nr 1 x2 = (-b - sqrt(delta))/(2*a) # pierwiastek nr 2 cat("x1 = ",x1,"\n","x2 = ",x2,"\n") # możemy wyświetlić wyniki } else { # jeżeli delta = 0 to x1=x2 (pierwiastek podwójny) x1 = -b/(2*a) # wytarczy obliczyć tylko jeden pierwiastek cat("x1 = x2 = ",x1,"\n") # } # <== ten nawias kończy drugą instrukcje złozoną # W powyższym przykładzie wykorzystano fakt, że liczby rzeczywiste # są podzbiorem liczb zespolonych, gdy delta > 0 część urojona # liczby zespolonej będzie miała wartość 0 # # Jeżeli nie chcemy żeby maszyna wypisywała dla delta > 0 zbędne # informacje o częściach urojonych wyników wystarczy niewielka " zmiana w powyższym kodzie: if(delta > 0.0) { x1 = Re(x1); x2 = Re(x2) } # funkcja Re pozostawia # z liczby zespolonej # tylko część rzeczywistą # ^--- krótkie instrukcje możemy pisać # w jednej linii oddzielając je od siebie średnikami cat("x1 = ",x1,"\n","x2 = ",x2,"\n") # możemy wyświetlić wyniki, # P e t l e for (zmienna in wektor) instrukcja for (i in 1:4) cat(i) # "cat" drukuje na ekranie kolejne wartości # zmiennej i > for (i in 1:4) cat(i) 1234 # <== nie wyszło tak jak sobie to wyobrażaliśmy > for (i in 1:4) cat(i,"\n") # <== tutaj po wydrukowaniu każdej wartości i 1 # funkcja cat dodawała znak nowej linii "\n" 2 3 4 > while (war) instrukcja # wykonywana instrukcja powinna móc zmienić # wartość war (wielkość logiczna - jak długo # war ma wartość TRUE petla bedzie wykonywana) repeat instrukcja # wyjscie z pętli jedynie przez break! # To jest instrukcja, która łatwo może generować # petle nieskończone next - skok do nastepnego kroku petli break - przerywa pracę pętli, maszyna wykonuje następny rozkaz po petli ################################################################# # Pętla for (wiemy ile kroków zamierzamy wykonać) suma = 0 dane = rnorm(50) # bierzemy 50 liczb losowych for (i in 11:35) suma = suma + dane[i] # sumujemy ich wartości # zaczynajac od 11-tej # kończąc na 35-tej > suma [1] 7.925669 > ################################################################## # Pętla while (nie wiemy ile kroków będzie wykonane) # Stary (2500 lat) grecki sposób na obliczanie pierwiastka kwadratowego P = 1000 # liczba pierwiastkowana x1 = 60 # odgadnięta wartość pierwiastka (raczej bardzo kiepsko) x2 = 1000/x1 # dzielimy P przez odgadniętą wartość nie_udalo_sie = abs(x2 - x1) >= 1.0e-7 # patrzymy co wyszło while (nie_udalo_sie) { #<== początek pętli x1 = (x1+x2)/2 # obliczamy średnią arytmetyczną x1 i x2 x2 = P/x1 nie_udalo_sie = abs(x2 - x1) >= 1.0e-7 # cat(x1," ", x2," ", nie_udalo_sie,"\n") } # <== koniec pętli x1 = (x1+x2)/2 # jeszcze możemy poprawić wynik... cat(" pierwiastek z",P," =",x1) Regresja liniowa funkcja lm() najprostsze użycie: model = lm(zm_zal ~ zm_niez + 1) 1 reprezentuje stałą część modelu - funkcję definiującą wyraz wolny Jeżeli wiemy, że linia regresji przechodzi przez punkt (0,0) to usunięcie wyrazu wolnego z modelu zapisujemy symbolicznie: model = lm(zm_zal ~ zm_niez - 1) # minus oznacza tutaj eliminacje z modelu # a nie operację odejmowania parametry dodatkowe: dane = nazwa ramki danych, z której pochodzą zmienne modelu weights = wektor wag (jeżeli jest potrzebny) informacje, jakie możemy uzyskać z modelu: coef(model) - współczynniki regresji resid(model) - różnice pomiędzy danymi a wynikami modelu fitted(model) - wartości Yi wynikające z modelu w punktach Xi summary(model) - podsumowanie wyników predict(model,newdata=wektor nowych wartości x) - wartości Yi dla nowych Xi (innych niż uzyte do budowy modelu) deviance(model) - suma kwadratów różnic pomiędzy danymi a wynikami modelu na wykładzie były pokazane przykłady użycia modeli liniowych do danych z wartościami odbiegającymi od trendu pokazałem usuwanie takich danych z użyciem funkcji which(), pozwalającej znaleźć wskażniki elementów wektora, które spełniaja zadana relację dane1 = dane[-which(dane < 100)] # z wektora dane eliminuje sie wszystkie # elementy, które są mniejsze niż 100 pokazałem również przykład użycia funkcji lm do znalezienia współczynników wielomianu, którego stablicowane wartości weszły jako dane do modelu x = seq(-5,6,by=0.5) y = 2.4523*x^2 + 1.7252*x + 5.3421 plot(x,y) model = lm(y ~ 1 + x + I(x^2)) summary(model) instalujemy pakiet (polynom): install.packages() wybieramy repozytorium (Wiedeń, Wrocław, Oświęcim) wyszukujemy pakiet polynom w liście, zaznaczamy i wskazujemy myszą na OK linuks sciąga pakiet i kompiluje w lokalnej bibliotece uzytkownika zeby uzyc pakietu nalezy teraz zaladowac go do pamieci: library(polynom) # i mozymy juz użyć funkcji polynomial budujacej wielomian o zadanych # wspólczynnikach p = polynomial(coef(model)) # bierzemy wspólczynniki z modelu liniowego p = as.function(p) # p było tekstem, zamieniamy je na funkcje curve(p,-5,6, add=TRUE, col="red) # dorysowujemy do punktów z komendy plot # krzywą reprezentująca wielomian # komenda add=TRUE dodaje krzywą do rysunku # nie wymazując poprzedniego obrazu Przykład 2 - praca z danymi zawierjącymi błędy i nieścisłości czytamy dane z pliku ufc.csv (Upper Flat Creek - las testowy Uniwersytetu Idaho), plik mozna skopiowac z adresu: http://www.ms.unimelb.edu.au/~andrewpr/r-users/ Przykład pochodzi z podrecznika Robinson-Icebreaker.pdf (na krążku CD ) > getwd() [1] "C:/Users/mrozek/Documents" > setwd("d:/usr/mrozek/zajecia/osr_2010") > ufc <- read.csv("ufc.csv") > dim(ufc) # patrzymy jaki jest rozmiar danych [1] 637 5 > names(ufc) # jak sie nazywaja dane? [1] "plot" "tree" "species" "dbh" "height" # działka, # drzewo, # gatunek, # diameter at breast height (137 cm od ziemi), [mm] # wysokość [dm] # zmieniamy jednostki dbh z mm na cm i wysokość z dm na m ufc$dbh.cm <- ufc$dbh/10 ufc$height.m <- ufc$height/10 # wyznaczamy liczbe drzew kazdego gatunku table(ufc$species) DF ES F FG GF HW LP PP SF WC WL WP 10 77 3 1 2 185 5 7 4 14 251 34 44 ufc$species[is.na(ufc$dbh)] <- NA # 10 drzew bez nazwy oznaczamy jako NA ufc$species <- factor(ufc$species) # ponownie klasyfikujemy drzewa # F i FG??? to prawdopodobnie jeden gatunek (Fir i Grand Fir - jodła, zamieniamy na GF) ufc$species[ufc$species %in% c("F", "FG")] = "GF" ufc$species = factor(ufc$species) table(ufc$species) > ufc <- ufc[!is.na(ufc$height.m), ] # wyrzucamy drzewa z brakującymi danymi nt. wysokości # rysujemy relacje pomiedzy wysokością drzewa a jego średnicą dla kilku gatunków > plot(ufc$dbh[ufc$species=="DF"],ufc$height[ufc$species=="DF"]) > plot(ufc$dbh[ufc$species=="ES"],ufc$height[ufc$species=="ES"]) > plot(ufc$dbh[ufc$species=="GF"],ufc$height[ufc$species=="GF"]) > plot(ufc$dbh[ufc$species=="WC"],ufc$height[ufc$species=="WC"]) stosując lm() możemy sprawdzić, czy użycie regresji liniowej jest uzasadnione w każdym przypadku... ===================================================================== szukanie zer funkcji (funkcja uniroot() ) szukalismy zer trojmianu kwadratowego przy uzyciu funkcji uniroot, porównujac wyniki z rozwiązaniami dokładnymi do określenia przedziałów, w których szukano pierwiastków uzyto wykresu funkcji x = seq(-5,6,by=0.5) y = 2.4523*x^2 + 1.7252*x - 11.3475 plot(x,y) > plot(x,y) > abline(h=0) > library(polynom) > p = polynomial(coef=c(-11.3475,1.7252,2.4523)) > p=as.function(p) > x1 = uniroot(p,lower=-3, upper=-2) > x1 $root [1] -2.531413 $f.root [1] -0.0002312473 $iter [1] 4 $estim.prec [1] 6.103516e-05 > x2 = uniroot(p,lower=1, upper=2) > x2 $root [1] 1.827933 $f.root [1] 1.049456e-05 $iter [1] 4 $estim.prec [1] 6.103516e-05 porównujemy z danymi "analitycznymi" a=2.4523 b=1.7252 c=-11.3475 delta=b^2-4*a*c x1a = (-b + sqrt(delta))/(2*a) x2a = (-b - sqrt(delta))/(2*a) > x1a [1] 1.827932 > x2a [1] -2.531434 ># do 5 miejsca po przecinku wyniki sa poprawne Na ćwiczeniach warto zrobić przykład z podręcznika Soetaert ze str 37 - na zdolność buforową wody morskiej i ze strony 38 zadania 8.3.1 ####### Interpolacja wielomianowa 26.04.2010 library(polynom) # wprowadzamy do pamięci pakiet polynom x=c(1,2,3,4,5,6) # ustalamy współrzędne punktów y=c(2,3.5,5.75,6.4,8.2,7.7) plot(x,y) # rysujemy punkty p=poly.calc(x,y) # znajdujemy wielomian interpolacyjny p # możemy go obejrzeć # chcemy obliczyć wartość dla x=2 p(2) # nie działa! p jest tekstem p=as.function(p) # musimy zamienić tekst na funkcję p(2) # teraz możemy wyznaczać wartości p(x) curve(p,1,6,add=TRUE) # dorysowujemy wielomian do punktów # czy warto znajdować jeden wielomian dla wszystkich punktów? Runge = function(x) 1/(1+25*x^2) curve(Runge,-1,1,ylim=c(-0.5,1),col="red") abline(h=0) x=seq(-0.8,0.8,by=0.4) y=Runge(x) plot(x,y,pch=2,cex=0.5,col="green") curve(Runge,-1,1,ylim=c(-0.5,1),col="red") abline(h=0) x=seq(-0.8,0.8,by=0.4) # bierzemy 5 wartości x y=Runge(x) # obliczamy dla nich funkcję Rungego points(x,y,pch=2,cex=0.5,col="green") # dorysowujemy węzły dpo wykresu p5=poly.calc(x,y) # znajdujemy wielomian interpolacyjny p5 = as.function(p5) # zamieniamy go na funkcję curve(p5,-1,1,add=TRUE,col="green") # i rysujemy zielona linią x=seq(-0.8,0.8,by=0.2) # może trzeba wziąć więcej węzłów interpolacji... # teraz jest ich 9 y=Runge(x) points(x,y) p9=poly.calc(x,y) p9 = as.function(p9) curve(p9,-1,1,add=TRUE) # czy dało to zamierzony skutek? # Funkcje sklejane (spline functions) # używamy wtedy, gdy chcemy jedynie przewidzieć wartości np. danych # eksperymentalnych pomiędzy punktami pomiarowymi podstawowe narzędzie to funkcja spline s = spline(x,y,n=abc) # x,y - dane pomiarowe (węzły), abc - liczba punktów # w przedziale range(x), dla których chcemy wyznaczyć # wartości interpolowanej funkcji # zastosujmy funkcje spline do funkcji Rungego x= seq(-0.8,0.8,by=0.1) y = Runge(x) plot(x,y,ylim=c(-0.1,1.1)) z = spline(x,y,n=50) lines(z,col="green") curve(Runge,-0.8,0.8,col="red",add=TRUE) # dorysujemy wykres "prawdziwej" # funkcji Rungego ############# przykład danych pomiarowych # co trzy godziny mierzono prędkość wiatru w okresie doby czas = seq(0,24,by=3) wiatr = c(5,6,7,9,4,6,3,7,9) range(wiatr) # sprawdzamy zakres zmienności plot(czas,wiatr,ylim=c(2,11)) s24 = spline(czas,wiatr,n=24) # wyznaczyliśmy wartości prędkośći wiatru # w odstepach godzinnych points(s24,col="red",pch=2) # rysujemy interpolowane dane dla 24 godzin lines(spline(czas,wiatr,n=100),col="blue") # rysujemy linie oparta na 100 # punktach wyznaczonych przez # funkcje spline # oprócz rysunku możemy potrzebować wartości liczbowych names(s24) # sprawdzamy nazwy w s24 wiatr24=as.data.frame(s24)$y # z ramki danych wybieramy kolumne o nazwie y # funkcja poly.calc p = poly.calc(x) # buduje wielomian o pierwiastkach zawartych w wektorze x p123 = poly.calc(c(1,2,3)) > poly.calc(c(1,2,3)) > -6 + 11*x - 6*x^2 + x^3 > p = poly.calc(x,y) # buduje wielomian interpolacyjny przechodzący przez # węzły o wspólrzędnych z wektorów x, y # różniczkowanie i całkowanie wielomianów p <- poly.calc(1:5) # budujemy wielomian o węzłach 1:5 p ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 deriv(p) # wyznaczamy pochodna tego wielomianu ## 274 - 450*x + 255*x^2 - 60*x^3 + 5*x^4 integral(deriv(p)) - 120 # całkujemy pochodną i dodajemy wyraz wolny # aby odtworzyć oryginalny wielomian p ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 # Obliczanie całek i pochodnych # całka oznaczona - pole powierzchni pod krzywą # Dwa przypadki: # a) znamy funkcję podcałkową # b) znamy jedynie tabele wartości funkcji (dane pomiarowe) # a) - stosujemy wzory przyblizonego całkowania (interpolujemy wartości funkcji) # metoda trapezów, Simpsona, Gaussa # procedury adaptacyjne integrate(function(x) x^2,0,1) # funkcja integrate integrate(sin,0,pi) integrate(sin,-pi,pi) integrate(function (x) exp(-2*x^2),0,Inf) # granice mogą być nieskończone #b) Do całkowania danych pomiarowych możemy sie posłuzyć funkcją # integral z pakietu DierckxSpline library(DierckxSpline) # ładujemy pakiet do pamięci # powinien być zainstalowany # funkcją install.packages() # pozwala na instalacje pakietów z Internetu x = seq(0,pi,length=20) y = sin(x) plot(x,y,ylim=c(-1,1)) abline(h=0) z = curfit(x,y,knots=seq(0,pi,length=10)) # funkcja curfit dopasowuje f. sklejaną # opartą na węzłach "knots) # wzięliśmy tylko 10 węzłów lines(x, fitted(z), col = "blue") # fitted(z) oblicza wartości funkcji sklejanej # w punktach x yd = deriv(z,order=1,at=x) # yd zawiera wartości pochodnej (pierwsza) # w puntach zadanych wektorem at= # w naszym przypadku w punktach x lines(x,yd,col="green") # rysujemy pierwsza pochodna na zielono... integral(z) # całkujemy funkcje sklejana z, # nie podalismy granic więc całkujemy # w granicach [0,pi] w1 = integral(z, from=pi/4, to=0.75*pi) # teraz zadaliśmy granice w1 w2 = integrate(sin, pi/4, 0.75*pi) # liczymy b. dokładną procedurą numeryczną w2 w1 - w2 # o ile sie róznią? w1 - w2$value # integrate zwraca nie tylko wartośc całki ale i ocene błedu # całkowania numerycznego # musimy odjąć jedynie wartość całki