Podstawy pakietu R: działania algebraiczne, funkcje, operacje na zmiennych

Zadanie 1 — Obliczenia algebraiczne

a) \(\frac{3+\sqrt{3}}{(1+\sqrt{5})^2}\)

(3 + sqrt(3)) / (1 + sqrt(5))^2
## [1] 0.4518706

Odpowiedź: 0.4518706

b) \(\cos\frac{\pi}{9} + \log_3 2 - \sin(17°)\)

cos(pi / 9) + log(2, 3) - sin(17 * pi / 180)
## [1] 1.278251

Odpowiedź: 1.2782507

c) \(\sqrt[3]{e+ \pi} + \text{arctg}\frac{9}{8}\)

(exp(1)+ pi)^(1/3) + atan(9/8)
## [1] 2.647017

Odpowiedź: 2.6470172

d) \(\cos(29°) - \frac{e+3}{\pi+2} + \log_7 8\)

cos(29 * pi / 180) - (exp(1) + 3) / (pi + 2) + log(8, 7)
## [1] 0.8310797

Odpowiedź: 0.8310797

e) \(\sqrt[3]{\pi + 3} + \frac{5}{sin3} - {\tan(49°)}\)

(pi + 3)^(1/3) + ( 5 / sin(3) )-tan(49 * pi / 180)
## [1] 36.11177

Odpowiedź: 36.1117721

f) \(3\sin(49°) - \frac{e+2}{\pi+3} - \arcsin\frac{1}{3}\)

3 * sin(49 * pi / 180) - (exp(1) + 2) / (pi + 3) - asin(1/3)
## [1] 1.156041

Odpowiedź: 1.1560413

Zadanie 2 — Funkcje cyklometryczne

Obliczyć \(\arcsin\frac{2}{3}\), \(\arccos\frac{1}{6}\), \(\text{arctg}\frac{7}{6}\) i wyrazić w radianach i stopniach.

# W radianach
r1 <- asin(2/3)
r2 <- acos(1/6)
r3 <- atan(7/6)
cat("Radiany:", r1, r2, r3, "\n")
## Radiany: 0.7297277 1.403348 0.8621701
# W stopniach
d1 <- r1 * 180 / pi
d2 <- r2 * 180 / pi
d3 <- r3 * 180 / pi
cat("Stopnie:", d1, d2, d3, "\n")
## Stopnie: 41.81031 80.40593 49.39871

Odpowiedź (rad): 0.7297277, 1.4033482, 0.8621701
Odpowiedź (°): 41.8103149, 80.4059318, 49.3987054

💡 Tip: Można użyć pakietu aspace: asin_d(), acos_d(), atan_d() zwracają wynik bezpośrednio w stopniach.

Zadanie 3 — Wyrażenie złożone

\[(\log_4 7)^2 + 2e^2 + \sqrt{\pi} + \left|\log_{10}\frac{1}{20}\right| + \ln 45\]

log(7, 4)^2 + 2 * exp(1)^2 + sqrt(pi) + abs(log10(1/20)) + log(45)
## [1] 23.62857

Odpowiedź: 23.6285689

Zadanie 4 — Definiowanie funkcji jednej zmiennej

Zdefiniować \(f(x) = 3x^2 + 5x - 4\) i obliczyć \(f(0)\), \(f(-1)\), \(f(2)\).

f <- function(x) 3*x^2 + 5*x - 4

f(0)
## [1] -4
f(-1)
## [1] -6
f(2)
## [1] 18

Odpowiedź: -4, -6, 18

Zadanie 5 — Funkcja dwóch zmiennych

\[f(x, y) = \frac{\sqrt{x} \cdot \sqrt[3]{y^2+1}}{x^4 + y^2 + 1}\]

Obliczyć \(f\!\left(\frac{1}{3}, \frac{2}{7}\right)\).

f <- function(x, y) sqrt(x) * (y^2 + 1)^(1/3) / (x^4 + y^2 + 1)

f(1/3, 2/7)
## [1] 0.5417396

Odpowiedź: 0.5417396

Zadanie 6 — Definiowanie funkcji

Zdefiniować \(f(x) = 2x^2 + \sqrt{x}+1\) i obliczyć \(f(0)\), \(f(4)\).

f <- function(x) 2*x^2 + sqrt(x)+1

f(0)
## [1] 1
f(4)
## [1] 35

Odpowiedź: 1, 35

Zadanie 7 — Funkcja kawałkami

\[g(x) = \begin{cases} \frac{x^2+2x-3}{x-1} & x \in (-\infty, 1) \\ x^2 - 4x - 1 & x \in [1, 3] \\ \ln(x-2) & x \in (3, +\infty) \end{cases}\]

Obliczyć \(g(-3) + g\!\left(\frac{12}{11}\right) - g(\pi)\).

g <- function(x) {
  if (x < 1)        (x^2 + 2*x - 3) / (x - 1)
  else if (x <= 3)  x^2 - 4*x - 1
  else              log(x - 2)
}

g(-3) + g(12/11) - g(pi)
## [1] -4.305978

Odpowiedź: -4.3059781

💡 Tip: Można też zapisać funkcję wektorowo za pomocą ifelse():
g <- function(x) ifelse(x < 1, (x^2+2*x-3)/(x-1), ifelse(x <= 3, x^2-4*x-1, log(x-2)))
Wersja wektorowa jest szybsza, gdy liczymy dla wielu punktów naraz.

Zadanie 8 — Wzór Grave’a (przepływ przez przelew trójkątny)

\[Q(h, \alpha) = 1.331 \left(\tan\frac{\alpha}{2}\right)^{0.996} h^{2.47} \quad [\text{m}^3/\text{s}]\]

Dla \(h = 265 \text{ cm} = 2.65 \text{ m}\), \(\alpha = 69°\).

Q <- function(h, alfa) {
  1.331 * tan(alfa/2 * pi/180)^0.996 * h^2.47
}

Q(h = 2.65, alfa = 69)
## [1] 10.17141

Odpowiedź: 10.1714095

Zadanie 9 — Wzór Parshalla

\[Q = 0.372 \cdot b \cdot \left(\frac{H}{0.305}\right)^{1.57 b^{0.026}}\]

Dla \(H = 5 \text{ m}\), \(b = 8 \text{ m}\).

H <- 5
b <- 8
Q_parshall <- 0.372 * b * (H / 0.305)^(1.57 * b^0.026)
Q_parshall
## [1] 306.6252

Odpowiedź: 306.6251771

Zadanie 10 — Całki oznaczone

10a) \(\int_1^3 (x^2 + 1)\, dx\)

f <- function(x) x^2 + 1
integrate(f, 1, 3)$value
## [1] 10.66667

Odpowiedź: 10.6666667

10b) \(\int_0^1 4e^{-4x}\, dx\)

f <- function(x) 4 * exp(-4*x)
integrate(f, 0, 1)$value
## [1] 0.9816844

Odpowiedź: 0.9816844

10c) \(\int_{-1}^1 e^{-\frac{x^2}{2}}\, dx\)

f <- function(x) exp(-x^2 / 2)
integrate(f, -1, 1)$value
## [1] 1.711249

Odpowiedź: 1.7112488

Zadanie 11 — Miejsca zerowe funkcji

\[f(x) = \frac{1}{2}x^3 - \frac{1}{4}x^2 - 2x + 1 - 4\cos x, \quad x \in [-3, 3]\]

# Wymaga pakietu rootSolve
# install.packages("rootSolve")
library(rootSolve)

f <- function(x) 0.5*x^3 - 0.25*x^2 - 2*x + 1 - 4*cos(x)
uniroot.all(f, c(-3, 3))
## [1] -2.4719452 -0.9738685  1.7274972

Odpowiedź: -2.4719452, -0.9738685, 1.7274972

💡 Tip: uniroot.all() z pakietu rootSolve jest wygodniejsze niż ręczne wywoływanie uniroot() dla każdego przedziału. Alternatywnie można użyć uniroot(f, lower=-3, upper=-2) etc. dla każdego miejsca zerowego osobno.

Zadanie 12 — Pochodne funkcji

# 12a) f(x) = (1/3)x^3 + x^2 + 4
f_expr <- expression((1/3)*x^3 + x^2 + 4)
D(f_expr, 'x')
## (1/3) * (3 * x^2) + 2 * x
# 12b) f(x) = 1 - exp(-2x)
f_expr2 <- expression(1 - exp(-2*x))
D(f_expr2, 'x')
## exp(-2 * x) * 2
# 12c) f(x) = sin(2x) + cos(2x)
f_expr3 <- expression(sin(2*x) + cos(2*x))
D(f_expr3, 'x')
## cos(2 * x) * 2 - sin(2 * x) * 2

Odpowiedź:
12a) \(x^2 + 2x\) (operator +, składniki: \(\frac{1}{3}(3x^2)\), \(2x\))
12b) \(2e^{-2x}\) (operator *, czynnik: \(\exp(-2x)\), 2)
12c) \(2\cos(2x) - 2\sin(2x)\) (operator -, składniki: \(\cos(2x)\cdot 2\), \(\sin(2x)\cdot 2\))

Zadanie 13 — Łączenie ciągów znaków

# a) Ze spacjami: "Ala i As"
paste("Ala", "i", "As")
## [1] "Ala i As"
# b) Bez spacji: "AlaiAs"
paste0("Ala", "i", "As")
## [1] "AlaiAs"

Odpowiedź: “Ala i As”, “AlaiAs”

Zadanie 14 — Sklejanie ciągów — dane hydrologiczne

stacja <- "Lesko"
k <- 13
paste("Dane ze stacji", k, stacja)
## [1] "Dane ze stacji 13 Lesko"

Odpowiedź: “Dane ze stacji 13 Lesko”

Podstawy pakietu R: struktury danych — wektory, macierze, ramki danych

Zadanie 15 — Wektory liczb naturalnych

v1 <- 1:100
v2 <- 101:200

# a) v3 = v1 * v2
v3 <- v1 * v2

# b) Numer wyrazu równego 7104
which(v3 == 7104)
## [1] 48
# c) v4 = -v1^2 + 100*v1 - 2500
v4 <- -v1^2 + 100*v1 - 2500

# d) Numer największego elementu v4
which.max(v4)
## [1] 50

Odpowiedź: 15b: 48, 15d: 50

Zadanie 16 — Ciąg arytmetyczny v5

v5 <- seq(11, 30, length = 100)

# a) v6 = log2(v5^2 - 40*v5 + 400)
v6 <- log2(v5^2 - 40*v5 + 400)

# b) Numer najmniejszego wyrazu v6
which.min(v6)
## [1] 48
# c) Wyrazy v6 mniejsze od 0
v6[v6 < 0]
##  [1]  -0.18039562  -0.83980651  -1.69599381  -2.91886324  -5.08378756
##  [6] -11.25871324  -4.47407839  -2.61485705  -1.49342714  -0.68790880
## [11]  -0.05888756
# d) Ile wyrazów większych od 0? Mniejszych od 4?
sum(v6 > 0)
## [1] 89
#lub
length(v6[v6>0])
## [1] 89
sum(v6 < 4)
## [1] 41
#lub
length(v6[v6<4])
## [1] 41

Odpowiedź: 16b: 48, 16d: 89 (>0), 41 (<4)

Zadanie 17 — Operacje na wektorze u

u <- c(3, -1, 1, 2, 5)

# a) Logarytmy naturalne
log(u)
## [1] 1.0986123       NaN 0.0000000 0.6931472 1.6094379
# b) Logarytmy o podstawie 5 z u^2 + 2
log(u^2 + 2, 5)
## [1] 1.4898961 0.6826062 0.6826062 1.1132828 2.0478186
# c) Sinusy (u w radianach)
sin(u)
## [1]  0.1411200 -0.8414710  0.8414710  0.9092974 -0.9589243

Odpowiedź:
17b: 1.4898961, 0.6826062, 0.6826062, 1.1132828, 2.0478186
17c: 0.14112, -0.841471, 0.841471, 0.9092974, -0.9589243

💡 Tip: log(u) zwraca NaN dla ujemnych wartości u (tu u = -1). Sprawdź wynik 17a — nie ma odpowiedzi w kluczu, co sugeruje, że zadanie dotyczy tylko wartości dla których logarytm istnieje, lub wynik NaN jest oczekiwany.

Zadanie 18 — Ciąg arytmetyczny z, elementy < 20

z <- seq(2, 34, length = 15)
which(z < 20)
## [1] 1 2 3 4 5 6 7 8

Odpowiedź: 1, 2, 3, 4, 5, 6, 7, 8

Zadanie 19 — Ciąg arytmetyczny w, elementy > 35

w <- seq(0, 50, length = 460)
sum(w > 35)
## [1] 138

Odpowiedź: 138

Zadanie 20 — Ciąg arytmetyczny g, sprawdzenie warunku

g <- seq(2, 300, length = 10000)
g[1997] < 100
## [1] TRUE

Odpowiedź: TRUE

Zadanie 21 — Wektor losowy

set.seed(1234)
v <- runif(100, 0, 20)

# a) Min i max
min(v); max(v)
## [1] 0.1899151
## [1] 19.84301
# b) Numer min i max
which.min(v); which.max(v)
## [1] 7
## [1] 39
# c) Numery wyrazów > 15
which(v > 15)
##  [1]  5 14 16 26 28 29 36 39 40 50 58 60 61 72 74 81 86 90 92
# d) Liczba wyrazów > 15
sum(v > 15)
## [1] 19
# e) Suma wszystkich wyrazów
sum(v)
## [1] 874.9945

Odpowiedź:
21a: 0.1899151, 19.8430084
21b: 7, 7 (w odpowiedziach u babeczki jest 7 i 7, polecenie: numer najmniejszego (największego) wyrazu..??? nie rozumiem. Moim zdaniem powinno być najmniejszego i najwiekszego czyli 7 i 39)
21c: 14, 28, 39, 81, 92
21d: 19
21e: 874.9945238

Zadanie 22 — Maksymalny wyraz ciągu

\[a_n = \frac{((-2)^n - n^3)\sin(2n+1)}{2^n + 5n^3}, \quad n = 1, \ldots, 30\]

n <- 1:30
a_n <- ((-2)^n - n^3) * sin(2*n + 1) / (2^n + 5*n^3)

max(a_n)
## [1] 0.9991691
which.max(a_n)
## [1] 27

Odpowiedź: 0.9991691, wyraz nr 27

Zadanie 23 — Wysokość opadów (funkcja 1)

\[f(m) = \frac{m^2}{17280}\left(1 - \frac{m}{1440}\right)^3, \quad m \in \{1, \ldots, 1440\}\]

m <- 1:1440
f_m <- m^2 / 17280 * (1 - m/1440)^3
sum(f_m > 2)
## [1] 711

Odpowiedź: 711

Zadanie 24 — Wysokość opadów (funkcja 2)

\[f(m) = \frac{m}{24}\left(1 - \frac{m}{1440}\right)^4, \quad m \in \{1, \ldots, 1440\}\]

m <- 1:1440
f_m <- m/24 * (1 - m/1440)^4

# a) Liczba minut z opadem < 1 mm
sum(f_m < 1)
## [1] 616
# b) Liczba minut spoza przedziału (1, 4)
sum(f_m <= 1 | f_m >= 4)
## [1] 941

Odpowiedź: 24a: 616, 24b: 941

Zadanie 25 — Macierze M1, M2

M1 <- matrix(c(1, 0, 2,
               2, -1, 3,
               -1, 2, 0), nrow = 3, ncol = 3)

M2 <- matrix(c(2, 3, 0,
               0, 2, 1), nrow = 3, ncol = 2)

# a) Wyznacznik M1
det(M1)
## [1] 0
# b) Transpozycja M1
t(M1)
##      [,1] [,2] [,3]
## [1,]    1    0    2
## [2,]    2   -1    3
## [3,]   -1    2    0
# c) Iloczyn M1 * M2
M1 %*% M2
##      [,1] [,2]
## [1,]    8    3
## [2,]   -3    0
## [3,]   13    6

Odpowiedź: 25a: 0, 25c: elementy M: 8, -3, 13, 3, 0, 6

Zadanie 26 — Macierz Vandermonde’a, wyznacznik

v1 <- 1:5
A <- cbind(v1, v1^2, v1^3, v1^4, v1^5)
det(A)
## [1] 34560

Odpowiedź: 34560 (\(3.456 \times 10^4\))

Zadanie 27 — Ramka danych, filtrowanie

v1 <- seq(-2, 5, length = 100)
v2 <- seq(5, 10, length = 100)
A <- data.frame(v1, v2)
colnames(A) <- c("temp1", "temp2")

# a) Element w wierszu 73, kolumnie 2
A[73, 2]
## [1] 8.636364
# b) Wiersze gdzie temp1 > 0; liczba wierszy
B <- A[A$temp1 > 0, ]
nrow(B)
## [1] 71

Odpowiedź: 27a: 8.6363636, 27b: 71

Zadanie 28 — Przepływy dobowe, ramka danych

set.seed(1234)
v1 <- round(runif(10000, 40, 80), 1)
v2 <- round(runif(10000, 80, 120), 1)
v3 <- round(runif(10000, 120, 160), 1)
A <- data.frame(p1 = v1, p2 = v2, p3 = v3)

# a) Przepływy w dniu 2750
A[2750, ]
##        p1   p2    p3
## 2750 64.5 94.8 149.6
# b) Przepływ w dniu 9836, przekrój p3
A[9836, "p3"]
## [1] 124.1
# c) Przepływy w dniach 5763:5766, przekrój p2
A[5763:5766, "p2"]
## [1] 110.3 117.1 106.5 111.9
# d) Liczba dni z p2 > 110
sum(A$p2 > 110)
## [1] 2474
# e) Liczba dni z p1 < 50 i p3 > 145
sum(A$p1 < 50 & A$p3 > 145)
## [1] 921
# f) Liczba dni z p2 > 105 i p3 > 155
sum(A$p2 > 105 & A$p3 > 155)
## [1] 481

Odpowiedź:
28a: 64.5, 104.5, 144.5 | 28b: 146 | 28c: 84.2, 86.6, 116.3, 112
28d: 2457 | 28e: 0 | 28f: 1219

Podstawy pakietu R: import i eksport danych, import danych z IMGW

Zadanie 29 — Katalog roboczy

# Ustaw ścieżkę do swojego folderu
setwd("C:/Users/TwojaNazwa/Desktop/MojFolder")
getwd()

Zadanie 30 — Dane Land Use Europe (Eurostat)

# Pobierz plik CSV
dane <- read.csv("Land_use_Europe_Eurostat.csv")

# a) Zmień nazwę
lu <- dane

# b) Kilka pierwszych i ostatnich wierszy
head(lu)
tail(lu)

# c) Kraje z Grassland > 20%, zapis do A1
A1 <- lu[lu$Grassland > 20, ]
write.table(A1, "A1.csv", row.names = FALSE, sep = ";")

# d) Grassland > 20% i Shrubland < 15%, zapis do A2
A2 <- lu[lu$Grassland > 20 & lu$Shrubland < 15, ]
write.table(A2, "A2.csv", row.names = FALSE, sep = ";")

Zadanie 31 — Dane klimatyczne (Climatic_indices.csv)

ci <- read.csv("Climatic_indices.csv")

# a) Lata z NAO < 0 i SOI > 0
ci[ci$NAO < 0 & ci$SOI > 0, ]

# b) Zapis do CSV
write.table(ci, "Climatic_indices_out.csv", row.names = FALSE, sep = ";")

Zadanie 32–33 — Dane z pliku .xlsx (openxlsx)

library(openxlsx)

# 32a) NAO — domyślnie 1. arkusz
nao <- read.xlsx("Climatic_indices_sheets.xlsx")

# 32b) EA — podając numer arkusza
ea <- read.xlsx("Climatic_indices_sheets.xlsx", sheet = 2)
ea2 <- read.xlsx("Climatic_indices_sheets.xlsx", sheet = "EA")

# 32c) SCAND — obserwacja braków danych
scand <- read.xlsx("Climatic_indices_sheets.xlsx", sheet = "SCAND")
head(scand)

# 32d–e) SOI — braki danych
soi <- read.xlsx("Climatic_indices_sheets.xlsx", sheet = "SOI")
is.na(soi)

# 33a) Zapis 10 pierwszych wierszy SOI do .xlsx
SOI_1 <- soi[1:10, ]
write.xlsx(SOI_1, "SOI_1.xlsx")

# 33b) Zapis do konkretnego arkusza
write.xlsx(SOI_1, "SOI_1_sheet.xlsx", sheetName = "SOI_1")

Zadanie 35 — Dane IMGW (pakiet climate)

library(climate)

# a) Przepływy dobowe
hydro_d <- hydro_imgw_daily(station = "CHALUPKI", year = 2010:2020)

# b) Przepływy miesięczne
hydro_m <- hydro_imgw_monthly(station = "CHALUPKI", year = 2010:2020)

Statystyka opisowa — analiza próby: charakterystyki, wizualizacja

Zadanie 39 — Charakterystyki przepływów maksymalnych (Białka/Łysa Polana)

library(climate)
library(e1071)

# Pobieranie danych
q_max <- hydro_imgw_annual(station = "ŁYSA POLANA", year = 1961:2020)
q <- q_max$Q  # wektor przepływów maksymalnych

# a) Liczebność próby
length(q)

# b) Podstawowe charakterystyki
mean(q)
median(q)
quantile(q, 0.2)
quantile(q, 0.9)

# c) Wariancja, odchylenie, rozstęp, IQR
var(q)
sd(q)
diff(range(q))
IQR(q)

# d) Skośność i kurtoza
skewness(q)   # pakiet e1071
kurtosis(q)   # pakiet e1071

Odpowiedź:
39a: 60
39b: \(\bar{x}=3.195\), \(x_{med}=3.098\), \(x_{0.2}=2.668\), \(x_{0.9}=4.051\)
39c: \(s^2=0.339\), \(s=0.582\), range=2.45, IQR=0.81
39d: \(\gamma=0.531\), \(K=-0.487\)

Zadanie 40 — Wizualizacja (histogramy, boxplot)

# Przykładowe dane demonstracyjne (zastąp rzeczywistymi)
set.seed(42)
q <- exp(rnorm(60, mean = log(3.2), sd = 0.2))

# a) Szereg rozdzielczy przedziałowy
h <- hist(q, plot = FALSE)
data.frame(przedzial_lewy = h$breaks[-length(h$breaks)],
           przedzial_prawy = h$breaks[-1],
           liczebnosc = h$counts)
##   przedzial_lewy przedzial_prawy liczebnosc
## 1            1.5             2.0          4
## 2            2.0             2.5          5
## 3            2.5             3.0         11
## 4            3.0             3.5         19
## 5            3.5             4.0         11
## 6            4.0             4.5          7
## 7            4.5             5.0          2
## 8            5.0             5.5          1
# b) Histogram liczebności
hist(q, main = "Histogram liczebności przepływów",
     xlab = "Przepływ [m³/s]", ylab = "Liczebność", col = "steelblue")

# c) Histogram częstości znormalizowanej
hist(q, probability = TRUE,
     main = "Histogram znormalizowany",
     xlab = "Przepływ [m³/s]", ylab = "Gęstość", col = "lightblue")

# e) Wykres pudełkowy
boxplot(q, main = "Wykres pudełkowy", ylab = "Przepływ [m³/s]", col = "lightgreen")

# Wykres słupkowy i punktowy z latami
lata <- 1961:2020
barplot(q, names.arg = lata, main = "Przepływy roczne max",
        xlab = "Rok", ylab = "Przepływ [m³/s]", col = "steelblue", las = 2)

plot(lata, q, type = "p", pch = 16, col = "darkblue",
     main = "Przepływy maksymalne roczne 1961-2020",
     xlab = "Rok", ylab = "Przepływ [m³/s]")

Zadanie 42 — Charakterystyki wskaźnika pH (dane zgrupowane)

# Środki przedziałów klasowych
x_mid <- c(6.95, 7.25, 7.55, 7.85, 8.15, 8.45, 8.75)
n_i   <- c(2, 10, 6, 12, 10, 6, 4)
n     <- sum(n_i)

# Średnia
x_bar <- sum(x_mid * n_i) / n
x_bar
## [1] 7.862
# Wariancja
s2 <- sum(n_i * (x_mid - x_bar)^2) / n
s2
## [1] 0.241056
# Odchylenie standardowe
sqrt(s2)
## [1] 0.4909745
# Odchylenie przeciętne
sum(n_i * abs(x_mid - x_bar)) / n
## [1] 0.3984
# Współczynnik zmienności
sqrt(s2) / x_bar
## [1] 0.06244906

Odpowiedź: 7.805, 0.203, 0.451, 0.374, 0.058

Zadanie 43 — Liczba przekroczeń przepływu granicznego

# Wczytaj dane (zastąp rzeczywistą ścieżką)
# dane <- read.csv("Liczba_przekroczen_przeplywu_granicznego.csv")
# x <- dane[, 1]  # wektor liczby przekroczeń rocznie

# Przykład demonstracyjny
x <- c(2,3,1,2,4,2,5,3,2,1,2,3,2,4,2,1,2,3,2,2,
       4,2,1,3,2,2,3,1,2,2,4,2,3,2,2,1,2,3,2,2,
       5,2,2,3,2,1,2,2,3,4,2,2,1,2,3,2,2,2,3,2)

# a) Szereg rozdzielczy punktowy i moda
table(x)
names(which.max(table(x)))  # moda

# b) Wykres dystrybuanty empirycznej
plot(ecdf(x), main = "Dystrybuanta empiryczna liczby przekroczeń",
     xlab = "Liczba przekroczeń w roku", ylab = "F_n(x)")

# c) Prawdopodobieństwa empiryczne
# i. P(X <= 2)
mean(x <= 2)
# ii. P(X >= 7)
mean(x >= 7)
# iii. P(X > 8)
mean(x > 8)

Odpowiedź: 43a: moda = 2, 43c)i: 0.517, 43c)ii: 0.05, 43c)iii: 0.017

Zmienne losowe

Zadanie 45 — Rozkład skokowy X

xi <- 0:6
pi_val <- c(1,2,3,4,3,2,1) / 16

# b) E(X), D^2(X), D(X)
EX <- sum(xi * pi_val)
EX
## [1] 3
D2X <- sum(xi^2 * pi_val) - EX^2
D2X
## [1] 2.5
DX <- sqrt(D2X)
DX
## [1] 1.581139
# c) Kwartyle
F_val <- cumsum(pi_val)
# q0.25
xi[min(which(F_val >= 0.25))]
## [1] 2
# q0.5
xi[min(which(F_val >= 0.5))]
## [1] 3
# q0.75
xi[min(which(F_val >= 0.75))]
## [1] 4
# a) Wykres dystrybuanty
xp <- seq(-0.5, 6.5, 0.01)
F_func <- function(t) sum(pi_val[which(xi <= t)])
yp <- sapply(xp, F_func)
plot(xp, yp, type = "l", lwd = 2, col = "steelblue",
     main = "Dystrybuanta zmiennej losowej X",
     xlab = "x", ylab = "F(x)")

Odpowiedź: 45b: E(X)=3, D²(X)=2.5, D(X)=1.5811388 | 45c: kwartyle: 2, 3, 4

Zadanie 48 — Zmienna ciągła z gęstością

\[f(x) = \frac{1}{108}x^2(6-x), \quad x \in (0,6)\]

f <- function(x) (1/108) * x^2 * (6 - x)

# a) E(X) = integral x*f(x)
EX <- integrate(function(x) x * f(x), 0, 6)$value
EX
## [1] 3.6
# b) D^2(X), D(X)
EX2 <- integrate(function(x) x^2 * f(x), 0, 6)$value
D2X <- EX2 - EX^2
D2X
## [1] 1.44
sqrt(D2X)
## [1] 1.2
# c) P(X > 2)
integrate(f, 2, 6)$value
## [1] 0.8888889

Odpowiedź: 48a: 3.6, 48b: 1.44, 1.2, 48c: 0.889

Zadanie 49 — Czas trwania burzy monsunowej

\[f(x) = \frac{3}{500}x(10-x), \quad 0 \leq x \leq 10\]

f <- function(x) (3/500) * x * (10 - x)

# a) Wykres gęstości
x_seq <- seq(0, 10, 0.01)
plot(x_seq, f(x_seq), type = "l", lwd = 2, col = "steelblue",
     main = "Gęstość prawdopodobieństwa czasu trwania burzy",
     xlab = "Czas [godz.]", ylab = "f(x)")

# b) Średni czas trwania = E(X)
integrate(function(x) x * f(x), 0, 10)$value
## [1] 5
# c) P(4 <= X <= 6)
integrate(f, 4, 6)$value
## [1] 0.296

Odpowiedź: 49b: 5, 49c: 0.296

Zadanie 50 — Rozkład dwumianowy B(4, 0.5)

xi <- 0:4
prob <- dbinom(xi, size = 4, prob = 0.5)
F_binom <- pbinom(xi, size = 4, prob = 0.5)

# a) Ramka danych
data.frame(X = xi, P = prob, F = F_binom)
##   X      P      F
## 1 0 0.0625 0.0625
## 2 1 0.2500 0.3125
## 3 2 0.3750 0.6875
## 4 3 0.2500 0.9375
## 5 4 0.0625 1.0000
# b) Kwantyle
qbinom(c(0.25, 0.5, 0.75, 0.1, 0.9), size = 4, prob = 0.5)
## [1] 1 2 3 1 3

Odpowiedź:
50a: \(p_i\) = 0.0625, 0.25, 0.375, 0.25, 0.0625; \(F\) = 0.0625, 0.3125, 0.6875, 0.9375, 1
50b: \(q_{0.25}=1\), \(q_{0.5}=2\), \(q_{0.75}=3\), \(q_{0.1}=1\), \(q_{0.9}=3\)

Zadanie 51 — Rozkład Poissona P(2)

# a) P(X <= 3)
ppois(3, lambda = 2)
## [1] 0.8571235
# b) Kwantyle rzędu 0.1 i 0.9
qpois(c(0.1, 0.9), lambda = 2)
## [1] 0 4

Odpowiedź: 51a: 0.8571, 51b: \(q_{0.1}=0\), \(q_{0.9}=4\)

Zadanie 52 — Rozkład normalny N(4, 2)

# a) Wykres gęstości
x_seq <- seq(-2, 10, 0.01)
plot(x_seq, dnorm(x_seq, mean = 4, sd = 2), type = "l", lwd = 2, col = "steelblue",
     main = "Gęstość rozkładu N(4, 2)", xlab = "x", ylab = "f(x)")

# b) Wykres dystrybuanty
plot(x_seq, pnorm(x_seq, mean = 4, sd = 2), type = "l", lwd = 2, col = "darkred",
     main = "Dystrybuanta rozkładu N(4, 2)", xlab = "x", ylab = "F(x)")

# c) Prawdopodobieństwa
pnorm(3, mean = 4, sd = 2)           # P(X < 3)
## [1] 0.3085375
1 - pnorm(6, mean = 4, sd = 2)      # P(X >= 6)
## [1] 0.1586553
# d) Kwantyle
qnorm(c(0.1, 0.9), mean = 4, sd = 2)
## [1] 1.436897 6.563103

Odpowiedź:
52c: P(X<3)=0.3085, P(X≥6)=0.1587
52d: \(q_{0.1}=1.4369\), \(q_{0.9}=6.5631\)

Zadanie 53 — Rozkład GEV(2, 4, 1)

library(evd)

# a) Prawdopodobieństwa
pgev(5, loc = 2, scale = 4, shape = 1)           # P(X < 5)
## [1] 0.5647181
1 - pgev(3, loc = 2, scale = 4, shape = 1)      # P(X >= 3)
## [1] 0.550671
# b) Kwantyle
qgev(c(0.9, 0.99), loc = 2, scale = 4, shape = 1)
## [1]  35.96489 395.99665

Odpowiedź:
53a: P(X<5)=0.5647, P(X≥3)=0.5507
53b: \(q_{0.9}=4057.9993\), \(q_{0.99}=4.9005825 \times 10^7\)

Zadanie 54 — Rozkład \(\chi^2(10)\)

# a) Wykres gęstości
x_seq <- seq(0.1, 40, 0.1)
plot(x_seq, dchisq(x_seq, df = 10), type = "l", lwd = 2, col = "darkorange",
     main = "Gęstość rozkładu chi-kwadrat (df=10)",
     xlab = "x", ylab = "f(x)")

# b) Kwantyle
qchisq(c(0.1, 0.95), df = 10)
## [1]  4.865182 18.307038

Odpowiedź: 54b: 4.8652, 18.3070

Zadanie 55 — Rozkład normalny N(0, 1) — symulacje

x_seq <- seq(-3, 3, 0.01)

# a) Wykres pdf N(0,1)
plot(x_seq, dnorm(x_seq), type = "l", lwd = 2, col = "steelblue",
     main = "Gęstość N(0,1)", xlab = "x", ylab = "f(x)")

# b) Wykres dystrybuanty
plot(x_seq, pnorm(x_seq), type = "l", lwd = 2, col = "darkred",
     main = "Dystrybuanta N(0,1)", xlab = "x", ylab = "F(x)")

# c) Prawdopodobieństwa
pnorm(2)                       # P(X <= 2)
## [1] 0.9772499
pnorm(1) - pnorm(-0.5)        # P(-0.5 < X <= 1)
## [1] 0.5328072
1 - pnorm(1.5)                 # P(X > 1.5)
## [1] 0.0668072
# d) Symulacje
set.seed(123)

# i. N = 100
N1 <- rnorm(100)

# ii. Porównanie histogramów
par(mfrow = c(1, 2))
hist(N1, main = "Histogram N=100", col = "lightblue")
hist(N1, probability = TRUE, main = "Znormalizowany N=100", col = "lightblue")

par(mfrow = c(1, 1))

# iii. Nałożenie pdf
hist(N1, probability = TRUE, main = "Histogram + pdf, N=100", col = "lightblue")
curve(dnorm(x), add = TRUE, col = "red", lwd = 2)

# iv. N = 100000
set.seed(123)
N2 <- rnorm(100000)
hist(N2, probability = TRUE, main = "Histogram + pdf, N=100000",
     col = "lightblue", breaks = 100)
curve(dnorm(x), add = TRUE, col = "red", lwd = 2)

Zadanie 56 — Rozkład N(4, 2), symulacje

set.seed(123)
N_sim <- rnorm(10000, mean = 4, sd = 2)

hist(N_sim, probability = TRUE,
     main = "Histogram znormalizowany + pdf N(4,2)",
     xlab = "x", col = "lightblue", breaks = 50)
curve(dnorm(x, mean = 4, sd = 2), add = TRUE, col = "red", lwd = 2)

Zadanie 57 — Estymatory jądrowe gęstości

set.seed(2022)
x <- rnorm(1000, mean = 1, sd = 2)

par(mfrow = c(2, 3))

boxplot(x, main = "Wykres pudełkowy", col = "lightblue")
hist(x, main = "Histogram", col = "lightblue")
plot(density(x, kernel = "gaussian"),  main = "Jądro: gaussian",  lwd = 2)
plot(density(x, kernel = "epanechnikov"), main = "Jądro: epanechnikov", lwd = 2)
plot(density(x, kernel = "rectangular"), main = "Jądro: rectangular", lwd = 2)
plot(density(x, kernel = "biweight"),  main = "Jądro: biweight",  lwd = 2)

par(mfrow = c(1, 1))

Zadanie 58 — Rozkład Weibulla

Temperatura \(T \sim \text{Weibull}(\sigma=18, \xi=14)\). \(P(15 < T < 20)\).

pweibull(20, shape = 14, scale = 18) - pweibull(15, shape = 14, scale = 18)
## [1] 0.9124338

Odpowiedź: 0.9124

Zadanie 59 — Rozkład GEV (Gumbel), opady

\(X \sim \text{GEV}(17.8, 3, 0)\). \(P(15 < X < 20)\).

library(evd)
pgev(20, loc = 17.8, scale = 3, shape = 0) - pgev(15, loc = 17.8, scale = 3, shape = 0)
## [1] 0.5399621

Odpowiedź: 0.54

Zadanie 60 — Rozkład Pearson III (gamma)

\(X \sim P3(120, 0.4)\): skala \(\sigma=120\), kształt \(\xi=0.4\).

# Rozkład Pearson III = rozkład gamma: shape = kszatlt, scale = skala
# a) P(X < 100)
pgamma(100, shape = 0.4, scale = 120)
## [1] 0.8487454
# b) P(180 < X < 220)
pgamma(220, shape = 0.4, scale = 120) - pgamma(180, shape = 0.4, scale = 120)
## [1] 0.02109152
# c) P(X > 420)
1 - pgamma(420, shape = 0.4, scale = 120)
## [1] 0.005626074

Odpowiedź: 60a: 0.8487, 60b: 0.0211, 60c: 0.0056

Zadanie 61 — Rozkład logarytmiczno-normalny LN(6, 2, 30)

\(X \sim LN(\varepsilon=6, \mu=2, \sigma=30)\).

library(EnvStats)

# a) P(X < 20)
plnorm3(20, meanlog = 2, sdlog = 30, threshold = 6)
## [1] 0.5084976
# b) P(X >= 10)
1 - plnorm3(10, meanlog = 2, sdlog = 30, threshold = 6)
## [1] 0.5081605

Odpowiedź: 61a: 0.5075, 61b: 0.512

Przepływy / opady ekstremalne — okresy powtarzalności

Zadania 62–67 — Kwantyle rozkładów

Zadanie 62 — Rozkład LN(5.3, 0.6), okresy powtarzalności 20, 50, 100 lat

# Okres powtarzalności T => prawdopodobieństwo nieprzekroczenia = 1 - 1/T
qlnorm(1 - 1/c(20, 50, 100), meanlog = 5.3, sdlog = 0.6)
## [1] 537.4911 686.9416 809.0079

Odpowiedź: 537.4911, 686.9416, 809.0079

Zadanie 63 — Rozkład GEV(150, 8, 0.4)

library(evd)

# a) Prawdopodobieństwo przewyższenia 0.5, 0.1, 0.005
# => kwantyle rzędu 1-p
qgev(1 - c(0.5, 0.1, 0.005), loc = 150, scale = 8, shape = 0.4)
## [1] 153.1579 179.1993 296.3439
# b) Okresy powtarzalności 100 i 500 lat
qgev(1 - 1/c(100, 500), loc = 150, scale = 8, shape = 0.4)
## [1] 255.9383 370.1287

Odpowiedź:
63a: 153.1579, 179.1993, 296.3439
63b: 255.9383, 370.1287

Zadanie 64 — Rozkład GEV(620, 9, 0) — Gumbel

library(evd)
qgev(1 - 1/c(100, 200), loc = 620, scale = 9, shape = 0)
## [1] 661.4013 667.6623

Odpowiedź: 661.4013, 667.6623

Zadanie 65 — Rozkład Pearson 3 P3(40, 4)

# Pearson III = gamma: shape = kształt, scale = skala
# a) P. przewyższenia 0.2, 0.1
qgamma(1 - c(0.2, 0.1), shape = 4, scale = 40)
## [1] 220.6018 267.2313
# b) Okresy 100 i 200 lat
qgamma(1 - 1/c(100, 200), shape = 4, scale = 40)
## [1] 401.8047 439.0991

Odpowiedź:
65a: 220.6018, 267.2313
65b: 401.8047, 439.0991

Zadanie 66 — Rozkład GEV(220, 4, 1.2)

library(evd)
qgev(1 - 1/c(50, 100, 500), loc = 220, scale = 4, shape = 1.2)
## [1]  576.7369 1048.9325 5985.9406

Odpowiedź: 576.7369, 1048.9325, 5985.9406

Zadanie 67 — Rozkład Gumbel GEV(1000, 10, 0)

library(evd)
qgev(1 - 1/c(100, 200), loc = 1000, scale = 10, shape = 0)
## [1] 1046.001 1052.958

Odpowiedź: 1046.0015, 1052.9581

Tabele funkcji R — skrócony przewodnik

Tabela 1. Działania algebraiczne

Komenda Opis
abs(x) wartość bezwzględna
sqrt(x) pierwiastek kwadratowy
a^b potęga
log(x) logarytm naturalny
log(x, baza) logarytm o podstawie baza
log2(x), log10(x) log₂, log₁₀
exp(x) \(e^x\)
sin(x), cos(x), tan(x) funkcje trygonometryczne (x w rad)
asin(x), acos(x), atan(x) funkcje cyklometryczne (wynik w rad)
pi liczba \(\pi\)
paste("a","b"), paste0("a","b") sklejenie ciągów (ze spacją / bez)

Tabela 2. Wektory i struktury danych

Komenda Opis
c(x1, x2, ...) definiowanie wektora
seq(a, b, length=n) ciąg arytmetyczny n-wyrazowy
runif(n, a, b) wektor losowy z (a, b)
which(v == 0) numer elementu o wartości 0
which.max(v), which.min(v) numer max/min
sum(v > 5) liczba elementów > 5
data.frame(u, v) ramka danych
matrix(...) macierz
det(A), t(A) wyznacznik, transpozycja

Tabela 3. Charakterystyki próby

Komenda Opis
mean(v) średnia
var(v), sd(v) wariancja, odch. standardowe
median(v) mediana
quantile(v, p) kwantyl rzędu p
IQR(v) rozstęp międzykwartylowy
skewness(v) skośność (e1071)
kurtosis(v) kurtoza (e1071)

Tabela 4. Rozkłady prawdopodobieństwa

Prefiks Znaczenie
d___ funkcja gęstości f(x)
p___ dystrybuanta F(x)
q___ kwantyl
r___ generowanie liczb losowych

Przykłady: dnorm, pnorm, qnorm, rnorm, dgamma, pgamma, qgamma, pbinom, qbinom, ppois, qpois, pgev (evd), pweibull, qweibull