1 Dodatkowe czcionki

Prosty sposób na zmianę czcionek w wykresach (w tym ggplot2) możemy znaleźć na stronie

https://cran.rstudio.com/web/packages/showtext/vignettes/introduction.html

dodatkowo można też zajrzeć na stronie

https://r-graph-gallery.com/custom-fonts-in-R-and-ggplot2.html

1.1 Zadanie - czcionki

Zmienić czcionki na wybrane na wykresie

library(ggplot2)

x<-seq(-3.1,5.1, by =.01)
y1=dnorm(x)
y2=dnorm(x,sd=2)
y3=dnorm(x,sd=1/2)
y4=dnorm(x,mean=2)

xx=as.data.frame(cbind(x,y1, y2, y3, y4))

ggplot(xx)+
  geom_line(aes(x,y1,color="N(0,1)"),size=2)+
  geom_line(aes(x,y2,color="N(0,4)"),size=2)+
  geom_line(aes(x,y3,color="N(0,1/4)"),size=2)+
  geom_line(aes(x,y4,color="N(2,1)"),size=2)+
  labs(title="Wykresy rozkładów normalnych",
       subtitle="dla wybranych parametrów",
       caption="źródło inspiracji - dr K.N.")+
  scale_color_manual(name="opis",
                     breaks=c('N(0,1)', 'N(0,4)', 'N(0,1/4)', 'N(2,1)'),
                     values=c('N(0,1)'='blue', 'N(0,4)'='red', 'N(0,1/4)'='green', 'N(2,1)'='orange'))

2 Przydatne szablony

Na stronie

https://r-graph-gallery.com/

znajduje się zbiór ciekawych wykresów wraz z kodami w R

2.1 Zadanie - wykres typu Doughnut lub Piechart

Utworzyć wykres typu Doughnut albo Pie chart obrazujący podział mandatów w polskim parlamencie (można też wykorzystać inne dane).

2.2 Zadanie - grafy

Korzystając z

https://r-graph-gallery.com/network.html

utworzyć własny graf

3 Testowanie czasu wykonywania kody

Za pomocą polecenia system.time() możemy zmierzyć czas działania bloku instrukcji

3.1 Zadanie - czas działania kodu

Sprawdzić działanie poniższych instrukcji (wszystkie robią to samo, jednak w różnym czasie). Wyjaśnić różnice w sposobie implementacji. Uwaga: kod chwilę się wykonuje

system.time(
  {x=NULL
  for(i in 1:10^5) x=c(x,runif(1))
  }
)
system.time(
  {x=NULL
  for(i in 1:10^5) x[i]=runif(1)
  }
)
system.time(
  {x=numeric(10^5)
  for(i in 1:10^5) x[i]=runif(1)
  }
)
system.time(
  {x=NULL
  x=runif(10^5)
  }
)

3.2 Microbenchmark

Korzystając z pakietu microbenchmark możemy dostać wyniki z dokładnością do milisekundy. Uwaga: pakiet trzeba zainstalować:

#---------biblioteki
library(microbenchmark)
library(ggplot2)

#instrukcje testowe takie jak poprzednio
pomiary=microbenchmark(
  skladanie    = {x=NULL; for(i in 1:1000) x=c(x,runif(1))},
  skladanie2   = {x=NULL; for(i in 1:1000) x[i]=runif(1)},
  inicjacja    = {x=numeric(1000); for(i in 1:1000) x[i]=runif(1)},
  wektoryzacja = {x=NULL; x=runif(1000)}
)

#------ wynik
pomiary
## Unit: microseconds
##          expr    min      lq     mean  median     uq     max neval
##     skladanie 3905.5 4733.30 6810.377 5161.10 6845.1 18824.0   100
##    skladanie2 2701.8 2863.80 3578.712 2995.90 3818.4 16153.8   100
##     inicjacja 2489.5 2666.75 3423.268 2789.90 3577.9 13786.0   100
##  wektoryzacja   23.6   26.75   28.807   27.85   29.5    48.7   100
#------ na obrazku
autoplot(pomiary)

Uwaga: każdy blok instrukcji wykonywany jest 100 razy dla celów testowych, w tabelce i na wykresie mamy wyniki dla tych 100 wykonań.

3.3 Zadanie - microbenchmark

Zbudować i przetestować (porównać) za pomocą microbenchmarka kod wykonujący poniższe zadanie w dwóch wersjach: poprzez pętle i w wersji wektorowej.

Zadanie dla kodu:

  • przygotowanie (przed microbenchmarkiem): wylosować z rozkładu normalnego dwa wektory \(x\) i \(y\) o długości \(N=1000\) (a później \(N=10.000\) i może dłuższy).

  • właściwe zadanie: obliczenie długości różnicy \(x-y\)

3.4 Testowanie składowych bloków

Do testowania składowych bloku kodu można wykorzystać profiler:

Rprof("profiler.out", interval=0.01)
N=2000
df= as.data.frame(matrix(runif(2*N*N),2*N,N))
tmp =summary(lm(V1~.,data=df))
Rprof(NULL)
summaryRprof("profiler.out")$by.self
##                           self.time self.pct total.time total.pct
## "lm.fit"                       9.18    80.81       9.18     80.81
## "chol2inv"                     1.33    11.71       1.33     11.71
## "runif"                        0.20     1.76       0.20      1.76
## ".External2"                   0.14     1.23       0.28      2.46
## "[.data.frame"                 0.10     0.88       0.11      0.97
## "as.vector"                    0.09     0.79       0.09      0.79
## "matrix"                       0.08     0.70       0.28      2.46
## ".External"                    0.04     0.35       0.04      0.35
## "makepredictcall.default"      0.04     0.35       0.04      0.35
## "[["                           0.03     0.26       0.06      0.53
## "summary.lm"                   0.01     0.09       1.34     11.80
## "model.frame.default"          0.01     0.09       0.29      2.55
## "as.data.frame.matrix"         0.01     0.09       0.10      0.88
## "makepredictcall"              0.01     0.09       0.05      0.44
## "[[.data.frame"                0.01     0.09       0.03      0.26
## "FUN"                          0.01     0.09       0.03      0.26
## "<Anonymous>"                  0.01     0.09       0.02      0.18
## "deparse"                      0.01     0.09       0.02      0.18
## ".subset2"                     0.01     0.09       0.01      0.09
## "as.list.default"              0.01     0.09       0.01      0.09
## "dev.cur"                      0.01     0.09       0.01      0.09
## "is.na"                        0.01     0.09       0.01      0.09
## "pmatch"                       0.01     0.09       0.01      0.09

summaryRprof zwraca listę 4-elementową, gdzie by.self to informacja, która funkcja jak długo działała, uporządkowana zgodnie z czasem działania wyłącznie danej funkcji bez podfunkcji.

Graficznie: potrzebny dodatkowy pakiet profr

library(profr)
Rprof("profiler.out", interval=0.01)
tmp=table(outer(1:2000,1:2000,"*") %% 10)
tmp=sort(rnorm(10^7))
Rprof(NULL)
wynik=parse_rprof("profiler.out")
plot(wynik)

4 Trochę matematyki

Kilka przykładów “matematycznych” w R. No cóż, są lepsze narzędzia.

4.1 Wielomiany - pakiet polynom

Kilka przydatnych operacji na wielomianach:

library(polynom)
# 2 wielomiany
p1 = polynomial(c(2,0,1))
p2 = polynomial(c(2,2,1,1))
p1
## 2 + x^2
p2
## 2 + 2*x + x^2 + x^3
#działania na wielomianach
p1+p2
## 4 + 2*x + 2*x^2 + x^3
p1*p2
## 4 + 4*x + 4*x^2 + 4*x^3 + x^4 + x^5
p2/p1
## 1 + x
#całkowanie i różniczkowanie
integral(p1,c(0,1))
## [1] 2.333333
deriv(p2)
## 2 + 2*x + 3*x^2
library(polynom)
#definiowanie wielomianu poprzez zera lub punkty które ma zawierać
poly.calc(c(-1,1))
## -1 + x^2
poly.calc(c(0,2,4),c(3,2,3))
## 3 - x + 0.25*x^2
#pierwiastki wielomianu
library(polynom)
p1=poly.calc(c(-1,2,3))
p1
## 6 + x - 4*x^2 + x^3
solve(p1)
## [1] -1  2  3
#NWW i NWD dla dwóch wielomianów
library(polynom)
p1 = polynomial(c(2,0,1))
p2 = polynomial(c(2,2,1,1))

LCM(p1,p2)
## 2 + 2*x + x^2 + x^3
GCD(p1,p2)
## 2 + x^2
# i na koniec obliczenie wartości poprzez zamianę wielomianu na funkcję...
library(polynom)
p2 = polynomial(c(2,2,1,1))
as.function(p1)(2)
## [1] 6

4.2 Minima, maksima i pierwiastki funkcji

Do szukania ekstremów można wykorzystać funkcję optimize():

f=function(x) {(x-7)^2-x}
optimize(f,interval=c(01,10))
## $minimum
## [1] 7.5
## 
## $objective
## [1] -7.25
optimize(f,interval=c(01,10),maximum=TRUE)
## $maximum
## [1] 1.000046
## 
## $objective
## [1] 34.9994

Pierwiastki - funkcja uniroot

f=function(x) {(x-7)^2-x}
uniroot(f,interval=c(-1,10))
## $root
## [1] 4.807418
## 
## $f.root
## [1] -4.69523e-06
## 
## $iter
## [1] 8
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

4.3 Pochodne i całki

D(expression(3*(x-2)^2-15),"x")
## 3 * (2 * (x - 2))
a=D(expression(exp(x)*cos(x)),"x")
a
## exp(x) * cos(x) - exp(x) * sin(x)
typeof(a)
## [1] "language"
x=1:5
eval(a)
## [1]  -0.8186613  -9.7937820 -22.7190020   5.6322837 184.4161820
pochodna = deriv(~(x-2)^2+15,"x")
pochodna
## expression({
##     .expr1 <- x - 2
##     .value <- .expr1^2 + 15
##     .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
##     .grad[, "x"] <- 2 * .expr1
##     attr(.value, "gradient") <- .grad
##     .value
## })
x=1:3
eval(pochodna)
## [1] 16 15 16
## attr(,"gradient")
##       x
## [1,] -2
## [2,]  0
## [3,]  2
integrate(dnorm,0,Inf)
## 0.5 with absolute error < 4.7e-05
f=function(x){(sin(x))^2+cos(x)}
integrate(f, 0,10)
## 4.227743 with absolute error < 1.3e-12
#integrate(f, 0,Inf)