1 Wstęp i linki do laboratoriów z 1 stopnia

1.1 Wstęp

Laboratorium 1 zawiera podstawowe informacje dotyczące tematów omawianych na 1 stopniu.
Sekcje O pakiecie R i RStudio w zasadzie można pominąć, umieszczone są w celu powtórkowym. Warto zajrzeć w razie problemów z pamięcią lub instalacją środowiska programistycznego R.
Sekcje Struktury danych, Funkcje i pętle oraz Grafika podstawowa i pakiet ggplot2 zawierają sporo przykładów, warto do nich zajrzeć przy wykonywaniu zadań
Sekcja Zadania zawiera kilka dość nietrywialnych zadań powtórkowych łączących ze sobą elementy struktur danych, programistycznych i graficznych. No i trochę matematyki.
Sekcja Wybrane pakiety numeryczne zawiera wybrane przykłady pakietów numerycznych

1.2 Linki do laboratoriów z 1 stopnia

Poniżej lista laboratoriów z 1 stopnia w wersji archiwalnej.

Laboratorium 1: wstęp, wektory, testy logiczne, wyodrębnianie podzbiorów danych
Laboratorium 2: macierze, ramki danych, listy
Laboratorium 3: elementy programowania: instrukcje warunkowe, pętle, funkcje
Laboratorium 4: grafika podstawowa w R, zapis i odczyt plików
Laboratorium 5: Pakiet ggplot2

2 O pakiecie R*

R jest obiektowym, interpretowanym, interaktywnym językiem programowania. Jest obecnie podstawowym językiem wykorzystywanym w statystycznej analizie danych. Jest całkowicie darmowy (oparty o licencję GNU).
Pierwsza wersja R została stworzona na początku lat 90-tych przez Roberta Gentlemana i Rossa Ihakę (Wydział Statystyki, Uniwersytet w Auckland, Nowa Zelandia). Od 1997 roku rozwojem projektu kierował zespół ponad 20 ekspertów z różnych dziedzin tj. statystyka, analiza numeryczna czy informatyka. Od tamtego czasu grupa ta znacząco się powiększyła. Ponadto swoje funkcje piszą naukowcy z całego świata.

Na stronie projektu można znaleźć nie tylko najnowszą wersję R, ale również bogatą dokumentację oraz wiele darmowych podręczników: http://www.r-project.org.

2.1 R i RStudio są darmowe - skąd ściągnąć?

R i RStudio są darmowe (miło), można bez przeszkód korzystać z nich w domu. W przypadku zastosowań komercyjnych sugeruję dokładne zapoznanie się z licencją.

Najłatwiej, choć nie jedyny sposób, to użycie strony:

https://posit.co/download/rstudio-desktop/

Najpierw instalujemy R, a później RStudio.

3 RStudio*

RStudio jest edytorem, który ułatwia pracę z R zintegrowanym środowiskiem programistycznym dla języka R. Posiada prosty interfejs użytkownika zbudowany wokół paska menu i kilku okien, którego głównym celem jest zapewnienie skrótów do niektórych z najczęściej używanych poleceń. Oczywiście te same funkcje można wywołać, wpisując odpowiednie polecenia w oknie konsoli. RStudio jest bardzo przydatnym narzędziem ułatwiającym debugowanie, tworzenie pakietów, aplikacji i raportów.

Przykładowy wygląd okna RStudio

Okno edycji skryptów umożliwia tworzenie i zapisywanie skryptów.

W oknie Konsoli wyświetlane są wyniki, można tam też wpisywać instrukcje/polecenia.

W prawym dolnym oknie wyświetlane są rysunki (o tym na innych zajęciach) i pomoc.

Prawe górne okno Enviroment - środowisko zawiera m.in. informacje na temat utworzonych zmiennych.

3.1 Konsola

Pojedyncze polecenia (lub ich grupy) można wpisywać w konsoli, zatwierdzając je Enterem.

W konsoli pojawiają się też wyniki działania poleceń.

W konsoli, za pomocą strzałek góra i dół na klawiaturze, można wybrać z historii wcześniej wpisywane instrukcje.

W przypadku większej ilości instrukcji wygodniejsze jest utworzenie skryptu - ciągu poleceń.

3.2 Tworzenie nowego skryptu i zapisanie go

W menu File wybieramy New file a następnie R Script. Plik ten zapisujemy sobie pod wybraną nazwą w wybranym miejscu (tak żeby nie utracić zapisanej pracy).

3.3 Praca ze skryptami

Pojedyncze komendy lub ich grupy wykonuje się poprzez ich zaznaczenie i naciśnięcie Ctrl+Enter lub naciśnięciu przycisku Run (znajdującego się w górnej prawej części okna Edycji skryptów).

3.4 Wybrane ustawienia RStudio

3.4.1 Sprawdzanie pisowni

Warto włączyć sprawdzanie pisowni, czy to w języku polskim czy angielskim. W tym celu wybieramy w menu Tools pozycje Global Options… a następnie:

3.4.2 Wygląd edytora

Wybieramy w menu Tools pozycje Global Options… a następnie:

3.5 Instalowanie i ładowanie dodatkowych pakietów

3.5.1 Polecenie w R do instalowania i ładowania bibliotek

Istotną zaletą R jest duża liczba wyspecjalizowanych bibliotek (pakietów).

Do instalowania bibliotek służy polecenie install.packages(), gdzie w nawiasie podajemy nazwę instalowanej biblioteki. Daną bibliotekę instalujemy tylko raz na danym urządzeniu. Przykład:

install.packages("ggplot2")

Po zainstalowaniu, gdy chcemy skorzystać z biblioteki, musimy za każdym razem dołączyć wybraną bibliotekę poleceniem library():

library(ggplot2)
#lub
library("ggplot2")

3.5.2 Alternatywny sposób instalowania pakietu

Być może wygodniejszy sposób to wybranie z menu Tools polecenia Install Packages... i poprzez okno wyszukiwarki znaleźć odpowiedni pakiet:

Gdy otworzymy skrypt lub plik Markdown, który zawiera niezainstalowany pakiet, pojawia się też monit z informacją o braku pakietu i sugestią jego zainstalowania:

3.6 Zarządzanie środowiskiem pracy

Możliwe jest usuwanie utworzonych obiektów za pomocą funkcji rm() lub z menu okna ‘środowisko/Enviroment’ (wyszukując obiekt, a nast. używając przycisku z miotłą). Wykonaj polecenie rm(a2) i sprawdź, czy zmienna a2 zniknęła ze środowiska pracy. Spróbuj usunąć inną zmienną ustawiając w oknie ‘środowisko’ Grid (kliknij strzałkę przy List i zmień na Grid), a następnie zaznacz zmienną do usunięcia i wykorzystaj przycisk miotły. Naciśnięcie przycisku miotły przy opcji List spowoduje usunięcie wszystkich zmiennych. Natomiast kombinacja Ctrl+L czyści konsolę.

Za pomocą polecenia rm(list = ls()) można wyczyścić całe środowisko pracy. Polecenie to można umieścić na początku skryptu R, aby mieć pewność, że pracujemy w “czystym” środowisku.

3.7 Komentarze

Aby dodać komentarz w kodzie/skrypcie R używamy znaku #. Komentarzem jest wtedy cała linijka występująca po tym znaku.

Aby zakomentować blok kodu (kilka linijek) w skrypcie, zaznaczamy wybrane linie kodu, a następnie korzystamy z kombinacji Ctrl+Shift+c. Ponowne użycie tej kombinacji usuwa komentarze z zaznaczonych linijek.

4 Struktury danych - krótki przegląd

4.1 Wektory

Wektory są podstawową struktura danych. Zawierają elementy tego samego typu. Każda liczba traktowana jest jako wektor jednoelementowy. Wektory budowane są przy użyciu funkcji c(). Funkcja ta służy też do łączenia wektorów:

a1=c(11,5,8,41,5,5)
a1
a2=c(-2,3,4,9,0)
a2
a3=c(a1,a2)
a3

4.1.1 Tworzenie wektorów

Inny sposób na tworzenie wektorów to polecenia rep, seq oraz losowania:

b1=numeric(10) # tworzy pusty (zerowy) wektor liczbowy długości 10
b1
b2=rep(c(1,2,3), each=4) # kopia każdego elementu
b2
b3=rep(c(1,2,3), times=4) # kopie wektora
b3
b4=seq(0,10,by=0.5) # podział odcinka co podaną jednostkę
b4
b5=seq(0,10,length.out=7) # podział odcinka na zadaną liczbę węzłów
b5
b6=sample(b5,size=3,replace=FALSE) # wylosowanie 3 elementów z wektora b5 bez powtórzeń. Można dodać argument prob  (help) zadający prawdopodobieństwo wystąpienia danego elementu
b6
b7=rnorm(10, mean=3, sd=7) # wylosowanie 10 liczb z rozkładu normalnego N(3,7)
b7

4.1.2 Zadanie

Wyjaśnić różnicę w działaniu poleceń 3:n-4 oraz 3:(n-4). Oczywiście, w składni poleceń różnica to nawias, ale dlaczego dostajemy inne wyniki?

n=15
3:n-4

##  [1] -1  0  1  2  3  4  5  6  7  8  9 10 11

3:(n-4)

## [1]  3  4  5  6  7  8  9 10 11

Ważne: Warto zapamiętać powyższą różnicę, bywa ona powodem błędów na przykład przy pisaniu pętli.

4.1.3 Wyodrębnianie podzbiorów danych

Odwołanie wprost do elementów wektora

c1=c("A","B","C","D","E","F","G","H","I","J","K","L","M","N")
c1[7] # wybór 7-mego elementu wektora
c1[3:6] # wybór elementów od 3 do 6
c1[c(4,1,3,8)] # wybór konkrtetnych elementów wektora wymaga funkcji `c`

4.1.4 Testy logiczne:

set.seed(1234) # aby losowanie było powtarzalne
d1=rnorm(25,mean=0,sd=2) # losujemy 25 liczb z N(0,2)
d1>1 # test logiczny, czy elementy wektora są większe od 1
which(d1>1) # podaje numery pozycji, gdzie elementy wektora są większe od 1 (czyli, gdzie było "TRUE" w poprzednim poleceniu)
d1[which(d1>1)] # wypisze te wartości, które były większe od 1
d1[d1>1] # działa tak samo jak wyżej
d1[d1>2 | d1<(-2)]=NA # nadpisanie "skrajnych" elementów wektora przez "NA" - not available
d1

4.1.5 Brakujące wartości w testach logicznych

Uwaga na brakujące wartości:

e1 <- c(1, 4, NA, 8)
e1[e1 > 5]

## [1] NA  8

traktuje brak danych NA jakby spełniał dowolne kryterium.

Zamiast tego należy użyć:

e1 = c(1, 4, NA, 8)
e1[!is.na(e1) & e1 > 5]

## [1] 8

4.1.6 Dodatkowe operatory logiczne

f1=c(-2,3,4,5)
any(f1>0) # alternatywa po wszystkich elementach
all(f1>0) # koniunkcja po wszystkich elementach
5 %in% f1 # czy element należy do wektora

4.2 Macierze, listy

4.2.1 Macierze

Macierze, podobnie jak wektory, muszą mieć elementy tego samego typu.

Tworzenie macierzy:

g1=1:20
g2=matrix(g1,nrow=5) #Tworzy macierz 5x4 z wektora `a` (wartości uporządkowane są kolumnami)
g2
g3=matrix(g1,nrow=5,byrow=T) #Tworzy macierz 5x4 z wektora `a` (wartości uporządkowane są wierszami)
g4=matrix(nrow=5, ncol=3) # Utworzenie pustej macierzy (jej elementy to `NA`)
g5=matrix(0,nrow=5, ncol=3) # Utworzenie macierzy wypełnionej zerami (można i inną liczbę użyć)

Dostęp do elementów macierzy:

g5[2,3] # wybór elementu
g5[,c(2,3)] # wybór kolumn
g5[c(1,2,5),] # wybór wierszy
g5[c(1,2,5),c(2,3)] # podmacierz

Uwaga: Operator mnożenia macierzowego ma postać

g2 %*% g3[1:4,] # mnożenie macierzowe "matematyczne"
g2 * g3 # mnożenie macierzowe "skalarne" - iloczyny odpowiadających sobie elementów

4.2.2 Listy

Lista jest uporządkowanym zbiorem obiektów, nie koniecznie tego samego typu. Elementami list mogą być dowolne obiekty zdefiniowane w R, którym można nadać nazwy.

Tworzenie listy i dostęp do jej elementów

h1=list(numerki=rnorm(6),y="blue",a=TRUE) # 3 elementowa lista, wektor, zmienna tekstowa i logiczna
h1[1] # dostęp poprzez indeks obiektu
h1$y # dostęp poprzez nazwę obiektu

4.3 Ramki danych

Ramki danych (data.frame) są tablicami, których kolumny mogą różnić się między soba typami danych. W ramach kolumny dane muszą być tego samego typu. Wszystkie kolumny muszą być tej samej długości.

4.3.1 Tworzenie ramki danych

Tworzenie ramki danych z kolumn danych

j1 = seq(0,10,length.out=100) # pierwszy wektor do ramki
j2 = runif(100) # drugi
j3 = rnorm(100)<0 # trzeci - co zwraca?
j4 = sample(colors(),size=100) # czwarty, co zwraca `colors()`?
j5 = data.frame(kol1=j1, kol2=j2, kol3=j3, kol4=j4) # budowa ramki
head(j5) # ramka jest duża, to wyświetlmy kilka pierwszych wierszy
summary(j5) # podsumowanie danych w ramce
colnames(j5) = c("zmienna_determinist.", "zmienna_losowa", "zmienna_logiczna", "zmienna_znakowa") # zmiana nazw kolumn na "rozsądniejsze" - można było to zrobić od razu
head(j5)
colnames(j5) # samo wyświetlenie nazw kolumn

#row.names można użyć do zmiany nazw wierszy

Ramka danych z macierzy:

j6 = matrix(rnorm(200), ncol=5) # losowa macierz
is.data.frame(j6) # czy jest to ramka danych?
j7=as.data.frame(j6) # zamiana na ramkę danych
is.data.frame(j7) # czy jest to ramka danych?
summary(j7) # takie tam podsumowanie

4.3.2 Dostęp do kolumn ramki danych

j5[,"zmienna_losowa"] # wybór 1 kolumny przez nazwę
j5[, c(1,4)] # wybór dwóch kolumn przez ich indeksy
j5[,c("zmienna_losowa", "zmienna_logiczna")] # wybór dwóch kolumn poprzez nazwy
j5$zmienna_determinist. # wybór kolumny poprzez operator `$`

4.3.3 Pakiet `dplyr`

Więcej informacji na temat operacji na ramkach danych zostanie omówione w laboratorium na temat pakietu dplyr

5 Funkcje i pętle

5.1 Funkcje

Podstawowy schemat budowy funkcji

nazwa_funkcji = function(argumenty){
  ciało funkcji
}

gdzie

argumenty - podajemy nazwy argumentów wykorzystywanych w funkcji, jeśli chcemy, możemy nadać im wartości domyślne,
ciało funkcji - ciąg instrukcji wykonywanych przez funkcję,
zwracana jest zazwyczaj wartość wyliczona w ostatniej linijce ciała funkcji,
inną zwracaną wartość możemy wymusić poleceniem return
jeżeli ma być zwracana więcej niż jeden obiekt, to wykorzystujemy listę zwracanych obiektów.

5.2 Instrukcje warunkowe i pętle

Poniżej umieszczone zostały schematy instrukcji warunkowych oraz pętli

5.2.1 `if`…`else`

Podstawowa wersja

if(warunek logiczny){
  blok instrukcji 1
} else{        # UWAGA: `else` musi być za nawiasem, a nie w nowej linijce
  blok instrukcji 2
}

Przy większej ilości warunków:

if(warunek logiczny 1){
  blok instrukcji 1
} else if(warunek logiczny 2){
  blok instrukcji 2
} else if(warunek logiczny 3){
  blok instrukcji 3
} else{
  blok instrukcji 4
}

Przy wielokrotnych warunkach można wykorzystać switch:

aa=sample(1:3,1)
switch(aa,
       '1' = "jedynka",
       '2' = "dwójka",
       '3' = "trójka")

## [1] "dwójka"

lub case_when z pakietu dplyr (już niedługo).

5.2.2 Pętla `for`

Poniżej kolekcja może być dość dowolnym wektorem, w szczególności nie koniecznie całkowitym, uporządkowanym, ani nawet numerycznym (może być to wektor ciągów znaków).

for (iterator in kolekcja) {
  blok instrukcji
}

5.2.3 Pętla `while`

while (warunek logiczny) {
  blok instrukcji
}

5.2.4 Pętla `repeat`

repeat {
  blok instrukcji
}

Zauważmy, że powyższa pętla będzie działać tak długo aż…. powiemy stop wykorzystując break .

6 Grafika podstawowa i pakiet `ggplot2`

6.1 Podstawowa grafika przy wykorzystaniu `plot`

Poniższy kod służy przypomnieniu tworzenia podstawowej grafiki w R oraz przeglądowi podstawowych argumentów.

par(mfrow=c(1,3)) #podział okna graficznego na 1*3 podokna
#Najbardziej podstawowy wykres - tutaj funkcja kwadratowa

#Dane do wykresu - wektory współrzędnych "x" oraz "y" punktów
x=seq(-2,2,by=0.1)
y=x^2

#Podstawowy wykres
plot(x,y,
     type='p', # typ wykresu, tutaj punktowy
     col='red') #kolor to col, nic innego chyba niedziała

#Etykiety
plot(x,y,
     type='l', # typ liniowy
     col='green',
     main="Więcej szczegółów", #nagłowek
     xlab="argumenty x", #etykiety osi
     ylab="funkcja x^2")

#dodatki
plot(x,y,
     type='p',
     col='blue',
     pch='*', #rodzaj znacznika
     main="dodatkowe elementy",
     xlab="argumenty x",
     ylab="funkcja x^2")
x1=c(-1,0,1)
y1=c(2,3,2)
points(x1,y1, # dodanie trzech punktów do już narysowanego wykresu - są to pary z x1 i y1
       col="red",
       pch=10, #numer znacznika
       cex=5) #rozmiar znacznika
lines(x1,y1-1) #dodanie linii łączącej te same punkty, tylko przesunięte o 1 w dół
abline(a=2,b=0.5,col="orange") #dodanie prostej y=b*x+a (dziwne, wiem)

par(mfrow=c(1,1)) # powrót do jednego okna

6.2 Grafika w `ggplot2`

Pakiet ggplot2 daje dużo większe możliwości tworzenia różnych wykresów. Wymaga zainstalowania (jednorazowo, o ile nie jest) oraz załadowania poleceniem library(ggplot2). Najwygodniej dane do zilustrowania na wykresie umieszczamy w ramce danych.

Poniżej seria przykładów z przypomnieniem podstawowych elementów.

6.2.1 Podstawowy wykres

#Na początku ładujemy bibliotekę
library(ggplot2) 

#Tworzymy ramkę danych
wart_x=seq(0,10, by=0.25)
wart_f1=sin(wart_x)
wart_f2=cos(wart_x)
df = data.frame(wart_x=wart_x,fsin=wart_f1,fcos=wart_f2) # Tworzymy ramkę danych

#Dość podstawowy wykres sinusa i cosinusa, punktowy i liniowy
#Poprzez umieszczenie koloru w `aes` i dodaniu `scale_color_manual` dostajemy legendę
rys=ggplot(df, aes(x = wart_x)) +
  geom_point(aes(y = fsin, color = "sinus"), size = 1.5) +
  geom_line(aes(y = fcos, color = "cosinus"), linetype="dotdash")+
  scale_color_manual(values = c("sinus" = "red", "cosinus" = "blue"))
rys #rysunek przypiszemy sobie do zmiennej `rys`, żeby potem dodawać wybrane elementy

6.2.2 Tytuły i etykiety

rys+
  labs(title="Tytuł wykresu",
       subtitle="A to podtytuł",
       x="argumenty są na osi X",
       y="wartości są na osi Y",
       color="Nazwa legendy", #Można zadać też w `scale_color_manual` przez `name`
       caption="Podpis obrazka")

6.2.3 Kolory tła, legendy

rys+
  theme(
        panel.border = element_rect(colour = "black", fill=NA, linewidth = 3), #ramka dookoła wykresu, fill - jeżeli będzie inne, to zamaluje nam cały wykres
        panel.background = element_rect(fill ='lightyellow'), #kolor tła
        panel.grid.major = element_line(color = 'darkgreen'), #kolor głównych linii siatki
        panel.grid.minor=element_line(color="lightgreen"), # linie pomocnicze
        legend.title=element_text(size=20, #rozmiar tytułu legendy
                                  color="blue", # kolor
                                  hjust=0), #
        legend.text=element_text(size=5, hjust=0.5), #rozmiar etykiet w legendzie, ich wyrównanie do środka
        legend.position = "inside", #pozycja legendy "inside"
        legend.position.inside = c(0.9, 0.8), #dokładne umiejscowienie legendy
        legend.box.background = element_rect(colour = "violet",linewidth=3), # kolor ramki legendy
        legend.key = element_rect(fill = "cornflowerblue"), # kolor tła pod znacznikami legendy
        legend.background = element_rect(fill="floralwhite", #kolor tła pod legendą
                                         linewidth =0.8,
                                         linetype="solid",
                                         color="green") # kolor ramki panelu legendy 
  )

7 Zadania

7.1 Zadanie - logarytm naturalny

Ciąg \(y_n=\sum_{i=1}^n \frac{1}{i}-\ln(n)\) jest zbieżny, gdy \(n\rightarrow\infty\). Oblicz wartość \(y_n\) dla n=50, 100, 1000. Wskazana wektorowa wersja (tj. bez pętli, bez funkcji).

Uwaga: logarytm naturalny jest poza sumą

Spodziewane wyniki

## [1] 0.5871823

## [1] 0.5822073

## [1] 0.5777156

Narysować wykres punktowy pierwszych tysiąca elementów ciągu \(y_n\).

Wskazówka - tutaj warto utworzyć funkcję

Przykładowy wykres

7.2 Zadanie - aproksymacja liczby \(\pi\)

Wiemy, że \(\frac{\pi}{4}=\sum_{k=0}^{\infty}\frac{(-1)^k}{2k+1}\). Znajdź metodę, która pozwala znaleźć przybliżenie liczby \(\pi\) z dokładnością do 6 miejsc po przecinku.

Dopuszczalna wersja z pętlą lub ręcznego dobierania ilości elementów

Wskazówki

Polecenie R pi służy do wyświetlenia przybliżonej wartości \(\pi\).
Można zacząć od wygenerowania wektorów postaci (1,-1,1,-1,...) oraz (1,3,5,7,...) o ustalonej długości, a dopiero później zająć się ilorazem i sumą.
Dokładność można sprawdzić obliczając różnicę (a dokładniej, jej wartość bezwzględną) między uzyskanym wynikiem, a wartością przybliżaną.

7.3 Zadanie - próba prosta 50 elementowa

Wygeneruj próbę 50-cio elementową \(x=(x_1,…,x_{50})\) z rozkładu jednostajnego na przedziale \([0,1]\) używając funkcji runif().
Utwórz wektor \(y=(y_1,...,y_{50})\), taki że \(y_i=0,\) gdy \(x_i<0.5\) oraz \(y_i=1\) w przeciwnym przypadku.
Wyznacz liczbę, a następnie proporcję elementów wektora \(x\), które są większe od 0.5 (w stosunku do całkowitej liczby elementów). Powtórz eksperyment, generując kilkukrotnie próbkę \(x\) i testując różne wartości jej długości \(n\). Czy wyniki są zgodne z Twoją intuicją?

7.4 Zadanie - liczby zespolone w R

7.4.1 Informacje i przykłady

Jednostka zespolona w R oznaczana jest przez 1i. Zatem przykładową liczbę zespoloną możemy zapisać jako

z=3+5*1i
z

## [1] 3+5i

Zadziała także

z=3+5i
z

## [1] 3+5i

Jednak nie zadziała kod:

x=3
y=5
z=x+yi #nie zadziała
z=x+y*i #też nie zadziała

a konieczne jest bezpośrednie użycie jednostki 1i, tzn.:

x=3
y=5
z=x+y*1i #tak już zadziała

Części rzeczywiste i urojone, moduł oraz argument dostaniemy przez:

Re(z)
Im(z)
Arg(z) # z dużej litery koniecznie
abs(z)

Przydatna obserwacja: jak wiadomo, funkcja arcus sinus asin w wersji rzeczywistej ma dziedzinę równą przedziałowi \([-1,1]\), jednak jeśli rozważymy zbiór liczb zespolonych, to dziedziną staje się cała płaszczyzna zespolona. Przetestować działanie funkcji asin w następujących przypadkach (w niektórych otrzymamy błąd - dlaczego?):

asin(0.5)
asin(5)
asin(5+0i)
sin(asin(5+0i)) #sprawdzenie

7.4.2 Zadanie właściwe

Zdefiniować wybraną funkcję wymierną, może to być funkcja \[ w(z) = \frac{1+z+\frac{1}{2}z^2}{1-z+\frac{1}{2}z^2} \]
Jak wiadomo, okrąg na płaszczyźnie zespolonej można sparametryzować za pomocą \[\exp(i \cdot \theta) = \cos(\theta) + i \cdot \sin(\theta), \quad \theta \in [0, 2\pi].\] Na płaszczyźnie zespolonej narysować krzywą będącą obrazem okręgu jednostkowego poprzez zdefiniowaną wyżej funkcję wymierną, tj. zbiór \[ w(\exp(i \cdot \theta)), \quad \theta \in [0, 2\pi]. \]
Poeksperymentować z innymi funkcjami, np.: \[ w_2(z) = \frac{24\cdot z^3(z-1)}{55z^3-59z^2+37z-9} \] \[ w_3(z) = \frac{720\cdot z^4(z-1)}{1901z^4-2774z^3+2616z^2-1274z+251} \]

Wskazówki

podzielić przedział \([0,2\pi]\) np. na 1000 części za pomocą seq, dostaniemy wektor;
na otrzymanym wektorze wykonać złożenie funkcji exp oraz w;
zbudować ramkę mającą 2 kolumny - część rzeczywista i część urojona wyniku;
narysować krzywą, w pakiecie ggplot2 może się przydać geom_path zamiast geom_line.

Przykładowe wykresy

7.5 Zadanie - paproć Barnsleya

Paproć Barnsleya jest fraktalem wyglądem przypominającym liść paproci. Wszystkie potrzebne informacje można znaleźć na Wikipedii Paproć Barnsleya.

Zadanie - napisać kod, który wygeneruje odpowiednie punkty, a następnie narysuje paproć Barnsleya przy użyciu pakietu ggplot2. Można przystosować przykładowy kod Matlaba z Wikipedii (prostsza wersja) lub napisać własny (prawidłowa wersja).

Przykładowy wynik

7.6 Zadanie - zbiór Mandelbrota

Zbiór Mandelbrota (żuk Mandelbrota) - kolejny znany fraktal. Szczegóły - Wikipedia zbiór Mandelbrota.

Jest to podzbiór płaszczyzny zespolonej, który tworzą punkty \(p \in\mathbb{C}\), o tej własności, że ciąg dany rekurencyjnie przez \[ \left\{ \begin{array}{l} z_0=0\\ z_{n+1}=z_n^2+p \end{array} \right. \] nie dąży do nieskończoności.

Równoważny warunek: \[ \forall n\in\mathbb{N} \quad |z_n|<2, \] może być wykorzystany do praktycznego wyznaczenia przybliżenia tego zbioru poprzez zastąpienie warunku \(\forall n \in \mathbb{N}\) skończoną liczbą (sprawdzamy np. 50 początkowych wyrazów ciągu, czy spełniają warunek \(|z_n|<2\)).

Zadanie Napisać kod, który narysuje zbiór Mandelbrota przy użyciu pakietu ggplot2. Zadbać o szczegóły typu kolor, opis osi itp.

Wskazówki

Liczby zespolone w R: jednostka zespolona to 1i, liczba zespolona to 5+3i. W przypadku zmiennych x+y*1i
Utworzyć funkcję testującą, czy liczba zespolona należy czy też nie do zbioru Mandelbrota. Przyjąć np. \(n=40\) pierwszych wyrazów uwzględnianych.
Utworzyć ramkę o 3 kolumnach: \(\mathrm{Re} \ z,\ \mathrm{Im} \ z,\) wynik testu.
Może przydać się funkcja expand.grid().
Wynik testu przełożyć na kolor przez aes.
Przyjąć zakres osi \([-2,1] \times [-1,1]\)
Wykres dla dużej ilości punktów może dość długo się tworzyć.

Przykładowy wynik

7.7 Zadanie - schemat Eulera

Rozważmy problem początkowy dla równań różniczkowych zwyczajnych postaci: \[ \left\{ \begin{array}{l} x'=f(t,x)\\ x(t_0) = \eta \end{array} \right. \] na przedziale \([t_0,T]\). Zakładamy, że funkcja \(f \colon [t_0,T] \times \mathbb{R}^n \to \mathbb{R}^n\) jest ciągła i spełnia warunek Lipschitza ze względu na drugą zmienną. Takie założenia gwarantują istnienie i jednoznaczność rozwiązania równania różniczkowego. Oznaczmy to rozwiązanie przez \(\varphi(t)\)

Najprostsza metoda numeryczna przeznaczona do rozwiązania tego typu równań to schemat Eulera:

zaczynamy od wybrania liczby \(N \in \mathbb{N}\),
dzielimy przedział \([t_0,T\) na \(N\) równych części o długości \(h=\frac{T-t_0}{N}\),
dostajemy węzły \(t_k=t_0+h \cdot k\), \(k=0,1, \ldots, N\),
startujemy od warunku początkowego \(x_0=\eta\),
kolejne aproksymacje \(x_k \approx \varphi (t_k)\) rozwiązania dokładnego w węzłach \(t_k\) dostajemy zgodnie z formułą

\[ x_{k+1}=x_k + h \cdot f(t_k,x_k), \quad k=0,1,\dots, N-1. \]

Uwagi:

powyższy algorytm jest wprost implementowalny, uwagi wymaga odpowiednia organizacja danych.
rozwiązaniem (analitycznym) równania różniczkowego jest funkcja ciągła, różniczkowalna oznaczona tutaj przez \(\varphi(t)\).
rozwiązaniem numerycznym jest ciąg punktów \(\{x_k\}_{k=0}^N\) aproksymujących rozwiązanie dokładne w węzłach \(t_k\) (patrz przykładowy wykres do zadań)
węzły numerujemy od \(0\) do \(N\), czyli ich jest \(N+1\). Wektory itd. w R numerujemy od \(1\)!

7.7.1 Równanie Prothero - Robinsona

Rozwiązać za pomocą schematu Eulera szczególny przypadek równania Prothero - Robinsona postaci: \[ \left\{ \begin{array}{l} x'=-5 \cdot (x-\cos t) - \sin t\\ x(0) = 1, \quad t \in [0{,}10], \end{array} \right. \] którego rozwiązaniem jest funkcja \(\varphi(t) = \cos t\).

Na wspólnym rysunku zilustrować rozwiązanie dokładne (linia ciągła) oraz ciąg punktów rozwiązania numerycznego metodą Eulera (punkty postaci \((t_k,x_k)\))

Przykładowe rozwiązanie

7.7.2 Układ Lotki - Volterry, model drapieżnik - ofiara (zadanie dodatkowe, nieobowiązkowe)

Rozwiązać numerycznie za pomocą schematu Eulera szczególny przypadek układu Lotki-Volterry postaci: \[ \left\{ \begin{array}{l} x'=(1.1-0.4 \cdot y) \cdot x\\ y'=(0.1 \cdot x -0.4) \cdot y\\ x(0) = 20\\ y(0) = 5, \quad t \in [0,100]. \end{array} \right. \] Układ ten nie ma rozwiązania analitycznego. Rozwiązanie zilustrować na wykresie. Ze względu na to, że jest to układ, czyli rozwiązaniem jest ciąg punktów \((t_k,x_k,y_k)\), możemy narysować dwie serie: \((t_k,x_k)\) oraz \((t_k,x_k)\). Drugi sposób zilustrowania to portret fazowy \((x_k,y_k)\)

Przykładowe rozwiązanie

Punkty z rozwiązania z pierwszego są połączone łamaną (geom_line) - ze względu na liczbę \(N=20.000\) punktów, a dzięki temu wykres jest “ładniejszy”.

7.8 Zadanie - Zbiór Julii (zadanie dodatkowe 23.10.2025)

Podobnym (w pewnym sensie) do zadania ze zbiorem Mandelbrota jest wykonanie rysunku zbioru Julii - opis na Wikipedii.

Jest to zadanie dodatkowe, ale można spróbować się z nim zmierzyć.

8 Wybrane pakiety numeryczne

8.1 Numeryczne rozwiązywanie równań nieliniowych

Rozwiązywanie równań nieliniowych jest ważnym zagadnieniem analizy numerycznej. W R można użyć w tym celu na przykład pakietu nleqslv, za pomocą którego rozwiążemy układ równań nieliniowych (Dennis Schnabel example 6.5.1 page 149): \[ \begin{array}{l} x_1^2+x_2^2-2=0\\ \exp(x_1-1)+x_2^3-2=0. \end{array} \]

#ładujemy pakiet
library("nleqslv")

#definiujemy funkcję
dslnex <- function(x) {
  y <- numeric(2)
  y[1] <- x[1]^2 + x[2]^2 - 2
  y[2] <- exp(x[1]-1) + x[2]^3 - 2
  y
}

#definiujemy jej macierz Jacobiego (nie jest wymagana)
jacdsln <- function(x) {
  Df <- matrix(numeric(2*2),2,2) #pusta macierz
  Df[1,1] <- 2*x[1]
  Df[1,2] <- 2*x[2]
  Df[2,1] <- exp(x[1]-1)
  Df[2,2] <- 3*x[2]^2
  
  Df
}

#punkt startowy
xstart <- c(2,0.5)

#argumenty - pierwszy to punkt startowy, drugi to funkcja, której miejsce zerowe szukamy
#jac - macierz Jacobiego jest opcjonalna
rozw=nleqslv(xstart, dslnex, jac=jacdsln)

#Wypisanie rozwiązania, zwracana jest lista
rozw[1]

#Sprawdzenie - unlist robi z listy wektor
dslnex(unlist(rozw[1]))

Więcej informacji w pomocy.

8.1.1 Zadanie 1

W powyższym przykładzie rozwiązanie jest średnio dokładne, ponieważ residuum jest rzędu \(10^{-9}\). Znaleźć w dokumentacji argument odpowiadający za zwiększenie dokładności.

8.1.2 Zadanie 2

Narysować wykresy funkcji \(f(x) = \sin(x)\) oraz \(g(x) = \frac{x}{3}\) na jednym wykresie. Przecinają się one w 3 punktach.
Jeden punkt przecięcia to \((0,0)\). Znaleźć numerycznie 2 pozostałe punkty. Można użyć wykresu do oszacowania punktu startowego.

Podpowiedź: Jedno z rozwiązań to

## $x
## [1] 2.278863

8.2 Numeryczna optymalizacja

Dla przykładu użyjemy funkcji Rosenbrocka: \[ f(x) = \sum_{i=1}^{n-1}\left[100(x_{i+1}-x_i^2)^2+(1-x_i)^2\right], \quad x\in\mathbb{R}^n \] która ma minimum w punkcie \((1,1,\dots,1) \in \mathbb{R}^n\) równe 0.

#definiujemy funkcję
Rosenbrock=function(x){
  n=length(x)
  wynik= sum(100*(x[2:n]-x[1:(n-1)]^2)^2+(1-x[1:(n-1)])^2)
  wynik
}

#ładujemy pakiet
library(optimx)

## Warning: pakiet 'optimx' został zbudowany w wersji R 4.5.1

#wywołanie funkcji do minimalizacji - argumetny: punkt statowy i funkcja do optymalizacji
optimx(c(0,0,0,0),Rosenbrock)

## Warning in max(logpar): brak argumentów w max; zwracanie wartości -Inf

## Warning in min(logpar): brak argumentów w min; zwracanie wartości Inf

##                    p1        p2        p3        p4        value fevals gevals
## Nelder-Mead 0.9999981 0.9999901 0.9999764 0.9999401 2.210632e-08    395     NA
## BFGS        0.9999340 0.9998686 0.9997378 0.9994749 9.053041e-08    241     50
##             niter convcode  kkt1 kkt2 xtime
## Nelder-Mead    NA        0 FALSE TRUE  0.01
## BFGS           NA        0  TRUE TRUE  0.00

8.2.1 Zadanie 1

Na stronie Wikipedia - funkcje testowe znajdziemy szereg funkcji służących testowaniu procedur numerycznych służących minimalizacji funkcji. Wybrać co najmniej jedną z nich i przetestować działanie polecenia optimx

8.2.2 Zadanie 2

Sprawdzić w pomocy dodatkowe parametry funkcji optimx, w szczególności - jak zwiększyć dokładność otrzymanego wyniku.

8.3 Numeryczne rozwiązywanie równań różniczkowych zwyczajnych

Jako przykład weźmiemy szczególny przypadek równania Lotki-Volterry. Jest to model interakcji drapieżnik-ofiara, gdzie występują cykliczne zmiany liczby osobników obu populacji:

\[ \left\{ \begin{array}{l} x'=x-0.1\cdot x \cdot y\\ y'=0.75\cdot x \cdot y-1.5 \cdot y\\ x(0)=20\\ y(0)=40 \end{array} \right. \] Rozwiązanie numeryczne dostaniemy za pomocą pakietu deSolve (prosty pakiet, trochę mało dokładny). Pakiet ten wymaga dość specjalnego przygotowania funkcji

#ładowanie pakietu
library(deSolve)

#klasyczne równanie ma postać x'(t) = f(t,x), gdzie t jest skalarem, x jest wektorem.
#wymagane są opcjonalne parametry params funkcji, można je potem ustawić jako NULL
lotka_volterra=function(t,u,parms){ 
  #t nie występuje jawnie, nasz układ jest autonomiczny
  #u jest wektorem [x,y]
  x=u[1]
  y=u[2]
  #sam układ
  dx=1*x-0.1*x*y
  dy=0.075*x*y-1.5*y
  #funkcja musi zwrócić listę
  return(list(c(dx,dy)))
}

#wartości początkowe, przedział całkowania
x0=20
y0=40
t0=0
T=100

#wartości początkowe zbieramy w wektor
u0=c(ofiara=20,lowca=40)

#przedział [t0,T] dzielimy na równe części
siatka=seq(t0,T,by=0.2)

#argumenty funkcji
#y=u0 - warunek początkowy
#times=siatka - punkty, w których rozwiązanie numeryczne dostaniemy - nasza siatka.
#func=lotka_volterra nasz układ
#method='lsoda' wybrana metoda numeryczna, jest ich więcej w pomocy
#wynikiem jest macierz, zapiszemy ją w zmiennej `wynik`

wynik=ode(y=u0,times=siatka,func=lotka_volterra,method='lsoda')
head(wynik)

##      time    ofiara    lowca
## [1,]  0.0 20.000000 40.00000
## [2,]  0.2 11.199796 37.16840
## [3,]  0.4  6.876026 31.41014
## [4,]  0.6  4.764360 25.34549
## [5,]  0.8  3.703705 19.98968
## [6,]  1.0  3.174585 15.58526

Wynik to macierz z aproksymacjami rozwiązania w zadanych przez siatkę punktach. Jest to “duża” macierz, dlatego dla ilustracji narysujemy nasz wynik. Pamiętamy, że to układ, czyli mamy 2 składowe zależne od czasu.

#Zaczniemy od zamiany na ramkę danych
wynik=as.data.frame(wynik)
library(ggplot2)
ggplot(wynik,aes(x=time))+
  geom_line(aes(y=ofiara,color="ofiary"))+
  geom_line(aes(y=lowca,color="łowcy"))+
  scale_color_manual(values = c("ofiary" = "red", "łowcy" = "blue"))+
  labs(x="czas", color="populacja")

Na powyższym wykresie mamy spore przekłamanie - wynik numeryczny to seria punktów, a nie linia ciągła.

Wynik możemy też zilustrować w przestrzeni fazowej, tj. płaszczyźnie drapieżnik-ofiara

ggplot(wynik, aes(x=ofiara,y=lowca))+
  geom_path(color="blue")

Na powyższym obrazku widać jakieś “dziwne fluktuacje”, zwiększymy dokładność i zagęścimy siatkę

u0=c(ofiara=20,lowca=40)
#zagęszczenie siatki
siatka=seq(t0,T,by=0.02)
#dodanie tolerancji - zwiększenie dokładności
wynik=ode(y=u0,times=siatka,func=lotka_volterra,parms=NULL,method='lsoda', rtol=1e-14, atol=1e-14)
wynik=as.data.frame(wynik)

#poprawiony rysunek

ggplot(wynik, aes(x=ofiara,y=lowca))+
  geom_path(color="blue")

8.3.1 Zadanie 1

Rozwiązać numerycznie równanie różniczkowe zwyczajne: \[ \left\{ \begin{array}{l} x'=\frac{1}{1+t^2}-2x^2\\ x(0) = 0 \end{array} \right. \quad \quad t \in [0,10] \] którego rozwiązaniem dokładnym jest funkcja \[\varphi(t) = \frac{t}{1+t^2}.\]
Porównać na jednym wykresie rozwiązanie numeryczne (jako punkty) oraz rozwiązanie dokładne (jako linia ciągła).

Przykładowe rozwiązanie

Przetwarzanie i wizualizacja danych w R

Laboratorium 1 - struktury danych, funkcje i grafika w R

Michał Braś

23 październik, 2025

1 Wstęp i linki do laboratoriów z 1 stopnia

1.1 Wstęp

1.2 Linki do laboratoriów z 1 stopnia

2 O pakiecie R*

2.1 R i RStudio są darmowe - skąd ściągnąć?

3 RStudio*

3.1 Konsola

3.2 Tworzenie nowego skryptu i zapisanie go

3.3 Praca ze skryptami

3.4 Wybrane ustawienia RStudio

3.4.1 Sprawdzanie pisowni

3.4.2 Wygląd edytora

3.5 Instalowanie i ładowanie dodatkowych pakietów

3.5.1 Polecenie w R do instalowania i ładowania bibliotek

3.5.2 Alternatywny sposób instalowania pakietu

3.6 Zarządzanie środowiskiem pracy

3.7 Komentarze

4 Struktury danych - krótki przegląd

4.1 Wektory

4.1.1 Tworzenie wektorów

4.1.2 Zadanie

4.1.3 Wyodrębnianie podzbiorów danych

4.1.4 Testy logiczne:

4.1.5 Brakujące wartości w testach logicznych

4.1.6 Dodatkowe operatory logiczne

4.2 Macierze, listy

4.2.1 Macierze

4.2.2 Listy

4.3 Ramki danych

4.3.1 Tworzenie ramki danych

4.3.2 Dostęp do kolumn ramki danych

4.3.3 Pakiet dplyr

5 Funkcje i pętle

5.1 Funkcje

5.2 Instrukcje warunkowe i pętle

5.2.1 if…else

5.2.2 Pętla for

5.2.3 Pętla while

5.2.4 Pętla repeat

6 Grafika podstawowa i pakiet ggplot2

6.1 Podstawowa grafika przy wykorzystaniu plot

6.2 Grafika w ggplot2

6.2.1 Podstawowy wykres

6.2.2 Tytuły i etykiety

6.2.3 Kolory tła, legendy

7 Zadania

7.1 Zadanie - logarytm naturalny

7.2 Zadanie - aproksymacja liczby \(\pi\)

7.3 Zadanie - próba prosta 50 elementowa

7.4 Zadanie - liczby zespolone w R

7.4.1 Informacje i przykłady

7.4.2 Zadanie właściwe

7.5 Zadanie - paproć Barnsleya

7.6 Zadanie - zbiór Mandelbrota

7.7 Zadanie - schemat Eulera

7.7.1 Równanie Prothero - Robinsona

7.7.2 Układ Lotki - Volterry, model drapieżnik - ofiara (zadanie dodatkowe, nieobowiązkowe)

7.8 Zadanie - Zbiór Julii (zadanie dodatkowe 23.10.2025)

8 Wybrane pakiety numeryczne

8.1 Numeryczne rozwiązywanie równań nieliniowych

8.1.1 Zadanie 1

8.1.2 Zadanie 2

8.2 Numeryczna optymalizacja

8.2.1 Zadanie 1

8.2.2 Zadanie 2

8.3 Numeryczne rozwiązywanie równań różniczkowych zwyczajnych

8.3.1 Zadanie 1

4.3.3 Pakiet `dplyr`

5.2.1 `if`…`else`

5.2.2 Pętla `for`

5.2.3 Pętla `while`

5.2.4 Pętla `repeat`

6 Grafika podstawowa i pakiet `ggplot2`

6.1 Podstawowa grafika przy wykorzystaniu `plot`

6.2 Grafika w `ggplot2`