• R is a free software environment for statistical computing and graphics
• RStudio is an integrated development environment (IDE) for R.

• geting help

# help on log function
?log

• R comes with a number of packages, some of them are loaded by default

# install or update a package (only once!)
install.packages("igraph")

# load a package (when you need it)
library(igraph)

# list all packages where an update is available
old.packages()

# update all available packages
update.packages()

• arithmetic: sum (+), minus (-), product (*), division (/), integer division (%/%), modulus (%%), exponent (^)
• comparison: equal (==), different (!=), less than (<), greater than (>), less than or equal to (<=), greater than or equal to (>=)
• logic operators: conjunction (&), disjunction (|), negation (!), exclusive disjunction (xor) (short-circuit versions && and ||)

Explain the following mismatch between math and R:

\[ (\sqrt{2}) ^ 2 \stackrel{?}{=} 2\]

sqrt(2) ^ 2 == 2

## [1] FALSE

The computer uses finite binary arithmetic and the binary representation of \(\sqrt{2}\) has infinite figures, hence it is rounded.

Define the xor operator in terms of conjunction (&), disjunction (|), and negation (!).

x = TRUE
y = TRUE
# first solution
(x | y) & !(x & y)

## [1] FALSE

# second solution
(x & !y) | (y & !x)

## [1] FALSE

x = TRUE
y = FALSE
(x | y) & !(x & y)

## [1] TRUE

(x & !y) | (y & !x)

## [1] TRUE

• the value NA (not available) is used to represent missing values;
• the value NULL is the null object (not to be confused with NULL in databases);
• the value Inf stands for positive infinity;
• the value NaN (not a number) is the result of a computation that makes no sense.

NA & TRUE

## [1] NA

NA & FALSE

## [1] FALSE

NA | TRUE

## [1] TRUE

NA | FALSE

## [1] NA

!NA

## [1] NA

2^1024

## [1] Inf

1/0

## [1] Inf

0 / 0

## [1] NaN

Inf - Inf

## [1] NaN

Of course, you may use variables to store values. There are 3 equivalent ways to assign a value to a variable:

x = 42  # my favourite
x <- 42 # this is the politically correct one!
42 -> x # used rarely

# print x
x

## [1] 42

# print structure of x (with type)
str(x)

##  num 42

R has four main atomic types:

# double (double-precision number)
x = 108.801
typeof(x)

## [1] "double"

# integer (integer number)
x = 108L
typeof(x)

## [1] "integer"

# character (a string of characters)
x = "108L"
typeof(x)

## [1] "character"

# logical (a Boolean, either TRUE or FALSE)
x = TRUE
typeof(x)

## [1] "logical"

The main data structures used in R include:

• atomic vector
• list
• matrix
• data frame

Dim Homogeneous     Heterogeneous
1d  atomic vector   list
2d  matrix          data frame

A vector is a sequence of elements with the same type. Vector indexes start at 1 (not 0).

# create a vector with c() function
c(1, 3, 5, 7)

## [1] 1 3 5 7

# concatenate vectors
c(c(1, 3), c(5, 7))

## [1] 1 3 5 7

# element-wise sum
c(1, 2, 3, 4) + c(10, 20, 30, 40)

## [1] 11 22 33 44

# recyclying
10 + c(1, 2, 3, 4)

## [1] 11 12 13 14

c(10, 20) + c(1, 2, 3, 4)

## [1] 11 22 13 24

# element-wise product
c(1, 2, 3, 4) * c(10, 20, 30, 40)

## [1]  10  40  90 160

# recyclying
10 * c(1, 2, 3, 4)

## [1] 10 20 30 40

c(10, 20) * c(1, 2, 3, 4)

## [1] 10 40 30 80

# scalar product (the result is a 1x1 matrix)
c(1, 2, 3, 4) %*% c(10, 20, 30, 40)

##      [,1]
## [1,]  300

x = c(TRUE, FALSE, TRUE, FALSE)
(y = !x) # also prints result

## [1] FALSE  TRUE FALSE  TRUE

x & y

## [1] FALSE FALSE FALSE FALSE

x | y

## [1] TRUE TRUE TRUE TRUE

xor(x, y)

## [1] TRUE TRUE TRUE TRUE

You may refer to members of a vector in several ways:

primes = c(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
primes[5]

## [1] 11

primes[c(1, 5, 10)]

## [1]  2 11 29

primes[-1]

## [1]  3  5  7 11 13 17 19 23 29

primes[-c(1, 5, 10)]

## [1]  3  5  7 13 17 19 23

primes > 15

##  [1] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE

primes[primes > 15]

## [1] 17 19 23 29

# modify the vector
primes[primes > 15] = Inf
primes

##  [1]   2   3   5   7  11  13 Inf Inf Inf Inf

Create the vector containing the integers from 1 to 60. Extract:

1. all numbers divisible by 4;
2. all numbers divisible by 3 or by 5;
3. all numbers divisible by both 4 and 6.
4. all even numbers that are not divisible by 3;

Try to solve the exercise without loops.

x = 1:60

# divisible by 4
x[x %% 4 == 0]

##  [1]  4  8 12 16 20 24 28 32 36 40 44 48 52 56 60

# divisible by 3 or by 5
x[x %% 3 == 0 | x %% 5 == 0]

##  [1]  3  5  6  9 10 12 15 18 20 21 24 25 27 30 33 35 36 39 40 42 45 48 50 51 54
## [26] 55 57 60

# divisible by both 4 and 6
x[x %% 4 == 0 & x %% 6 == 0]

## [1] 12 24 36 48 60

# even but not divisible by 3
x[x %% 2 == 0 & x %% 3 != 0]

##  [1]  2  4  8 10 14 16 20 22 26 28 32 34 38 40 44 46 50 52 56 58

Write a logical condition that is TRUE is the number is prime (Hint: take advantage of the all function).

n = 109
n == 2L || all(n %% 2:(n-1) != 0)

## [1] TRUE

n = 111
n == 2L || all(n %% 2:(n-1) != 0)

## [1] FALSE

# more efficient
n == 2L || all(n %% 2:sqrt(n) != 0)

## [1] FALSE

All elements of an atomic vector must be the same type, so when you attempt to combine different types they will be casted to the most flexible type (coercion). Types from least to most flexible are: logical, integer, double, and character.

# cast logical to numeric (integer or double)
sum(c(TRUE, TRUE, FALSE, FALSE))

## [1] 2

# cast numeric to character
paste(108.801, "108.801")

## [1] "108.801 108.801"

# cast character to numeric (invalid)
sum("abc", "def")

Without running R first, predict the type of each object:

a = c(TRUE, FALSE, TRUE)

b = c(TRUE, 2L, FALSE)

c = c(TRUE, 2L, 3.5)

d = c(TRUE, 2L, 3.5, "hello")

Then predict the result of:

sum(a)
sum(b)
sum(c)
sum(d)

Finally, verify your answers using typeof().

Vector elements can have names:

x = c(a = 1, b = 2, c = 3)
# or
x = c(1, 2, 3)
names(x) = c("a", "b", "c")
x["a"]

## a 
## 1

x[c("a", "b")]

## a b 
## 1 2

A factor is a vector that can contain only predefined values, and is used to store categorical variables (for instance sex or religion).

Factors are built on top of integer vectors using the levels attribute, which defines the set of allowed values.

x = factor(c("male", "female", "female", "male", "male"))
x

## [1] male   female female male   male  
## Levels: female male

typeof(x)

## [1] "integer"

levels(x)

## [1] "female" "male"

# if you use values that are not levels 
# a warning is issued and a NA is generated
x[1] = "unknown"
x

## [1] <NA>   female female male   male  
## Levels: female male

A list is a sequence of elements that might have different types.

# create a list
l = list(thing = "hat", size = 8.25, female = TRUE)

# print the list
l

## $thing
## [1] "hat"
## 
## $size
## [1] 8.25
## 
## $female
## [1] TRUE

str(l)

## List of 3
##  $ thing : chr "hat"
##  $ size  : num 8.25
##  $ female: logi TRUE

# an element
l$thing

## [1] "hat"

l[[1]]

## [1] "hat"

# a sublist
l[c("thing", "size")]

## $thing
## [1] "hat"
## 
## $size
## [1] 8.25

l[c(1, 2)]

## $thing
## [1] "hat"
## 
## $size
## [1] 8.25

“If list x is a train carrying objects, then x[[5]] is the object in car 5; x[5] is car number 5.”

# a sublist containing the first element of the list
l[1]

## $thing
## [1] "hat"

typeof(l[1])

## [1] "list"

# the first element of the list
l[[1]]

## [1] "hat"

typeof(l[[1]])

## [1] "character"

List elements can have any atomic or complex type. Hence a list can contain other lists, making it a nested list.

l = list(1, list(1, 2, 3), list("a", 1, list("TRUE", "FALSE")))
str(l)

## List of 3
##  $ : num 1
##  $ :List of 3
##   ..$ : num 1
##   ..$ : num 2
##   ..$ : num 3
##  $ :List of 3
##   ..$ : chr "a"
##   ..$ : num 1
##   ..$ :List of 2
##   .. ..$ : chr "TRUE"
##   .. ..$ : chr "FALSE"

Consider the list:

l = list(1, list(1, 2, 3), list("a", 1, list("TRUE", "FALSE")))

Find:

• the list list(1, 2, 3)
• the element 1 of list list(1, 2, 3)
• the element TRUE of list list("TRUE", "FALSE")

l = list(1, list(1, 2, 3), list("a", 1, list("TRUE", "FALSE")))
l[[2]]

## [[1]]
## [1] 1
## 
## [[2]]
## [1] 2
## 
## [[3]]
## [1] 3

l[[2]][[1]]

## [1] 1

l[[3]][[3]][[1]]

## [1] "TRUE"

Write a list containing the information of the Porphyrian Tree. Then select the insensitive part of the tree.

substance = 
  list(immaterial = "spirit", 
       material = list(
                   body = list(
                     inanimate = "mineral", 
                     animate = list(
                       living = list(
                         insensitive = "plant", 
                         sensitive = list(
                           irrational = "beast", 
                           rational = 
                             list(human = c("Arendt", "Butler", "Barad"))))))))

str(substance)

## List of 2
##  $ immaterial: chr "spirit"
##  $ material  :List of 1
##   ..$ body:List of 2
##   .. ..$ inanimate: chr "mineral"
##   .. ..$ animate  :List of 1
##   .. .. ..$ living:List of 2
##   .. .. .. ..$ insensitive: chr "plant"
##   .. .. .. ..$ sensitive  :List of 2
##   .. .. .. .. ..$ irrational: chr "beast"
##   .. .. .. .. ..$ rational  :List of 1
##   .. .. .. .. .. ..$ human: chr [1:3] "Arendt" "Butler" "Barad"

substance$material$body$animate$living$insensitive

## [1] "plant"

Consider:

experiment = list(
  info = list(
    title = "Network experiment",
    year = 2026
  ),
  results = list(
    training = c(0.72, 0.78, 0.81),
    test = c(0.68, 0.73, 0.76)
  ),
  successful = TRUE
)

Extract:

1. "Network experiment";
2. the complete results list;
3. the test vector;
4. the second test value;
5. a sublist containing only info and successful.

Then change the year from 2026 to 2027.

# 1
experiment$info$title

## [1] "Network experiment"

# 2
experiment$results

## $training
## [1] 0.72 0.78 0.81
## 
## $test
## [1] 0.68 0.73 0.76

# alternatively
experiment[["results"]]

## $training
## [1] 0.72 0.78 0.81
## 
## $test
## [1] 0.68 0.73 0.76

# 3
experiment$results$test

## [1] 0.68 0.73 0.76

# 4
experiment$results$test[2]

## [1] 0.73

# 5
experiment[c("info", "successful")]

## $info
## $info$title
## [1] "Network experiment"
## 
## $info$year
## [1] 2026
## 
## 
## $successful
## [1] TRUE

# change year
experiment$info$year = 2027

A matrix is a 2-dimensional vector, that is a vector of vectors of the same type and length.

# by column (the default)
M = matrix(data = 1:9, ncol = 3)
M

##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

# by row
N = matrix(data = 1:9, nrow = 3, byrow = TRUE)
N

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9

nrow(M)

## [1] 3

ncol(M)

## [1] 3

dim(M)

## [1] 3 3

M

##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

# element in row 1 and column 2
M[1, 2]

## [1] 4

# first row
M[1, ]

## [1] 1 4 7

# first column
M[ ,1]

## [1] 1 2 3

# sub-matrix
M[1:2, 1:2]

##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5

M[-3, -3]

##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5

# diagonal
diag(M)

## [1] 1 5 9

P = matrix(data = runif(9), nrow = 3,)
P

##           [,1]       [,2]      [,3]
## [1,] 0.5948830 0.41407173 0.9760665
## [2,] 0.0543333 0.01581983 0.4619919
## [3,] 0.5793651 0.20779408 0.1268384

# add column
(cbind(P, c(0, 0, 0)))

##           [,1]       [,2]      [,3] [,4]
## [1,] 0.5948830 0.41407173 0.9760665    0
## [2,] 0.0543333 0.01581983 0.4619919    0
## [3,] 0.5793651 0.20779408 0.1268384    0

# add row
(rbind(P, c(0, 0, 0)))

##           [,1]       [,2]      [,3]
## [1,] 0.5948830 0.41407173 0.9760665
## [2,] 0.0543333 0.01581983 0.4619919
## [3,] 0.5793651 0.20779408 0.1268384
## [4,] 0.0000000 0.00000000 0.0000000

# modify matrix
P

##           [,1]       [,2]      [,3]
## [1,] 0.5948830 0.41407173 0.9760665
## [2,] 0.0543333 0.01581983 0.4619919
## [3,] 0.5793651 0.20779408 0.1268384

P = cbind(P, c(0, 0, 0))
P

##           [,1]       [,2]      [,3] [,4]
## [1,] 0.5948830 0.41407173 0.9760665    0
## [2,] 0.0543333 0.01581983 0.4619919    0
## [3,] 0.5793651 0.20779408 0.1268384    0

M

##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

N

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9

# element-wise sum
M + N

##      [,1] [,2] [,3]
## [1,]    2    6   10
## [2,]    6   10   14
## [3,]   10   14   18

# element-wise product
M * N

##      [,1] [,2] [,3]
## [1,]    1    8   21
## [2,]    8   25   48
## [3,]   21   48   81

# matrix product
M %*% N

##      [,1] [,2] [,3]
## [1,]   66   78   90
## [2,]   78   93  108
## [3,]   90  108  126

# matrix transpose
M

##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

t(M)

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9

# matrix inverse
C = matrix(c(1,0,1, 1,1,1, 1,1,0), nrow=3, byrow=TRUE)
D = solve(C)
D

##      [,1] [,2] [,3]
## [1,]    1   -1    1
## [2,]   -1    1    0
## [3,]    0    1   -1

D %*% C

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

C %*% D

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

# linear systems C x = b
C

##      [,1] [,2] [,3]
## [1,]    1    0    1
## [2,]    1    1    1
## [3,]    1    1    0

b = c(2, 1, 3)
# the system is:
# x1      + x3 = 2
# x1 + x2 + x3 = 1
# x1 + x2      = 3
x = solve(C,b)
x

## [1]  4 -1 -2

C %*% x

##      [,1]
## [1,]    2
## [2,]    1
## [3,]    3

# matrix spectrum
spectrum = eigen(C)
# columns are the eigenvectors
spectrum$vectors

##            [,1]       [,2]       [,3]
## [1,] -0.4151581 -0.4743098 -0.6026918
## [2,] -0.7480890 -0.2110877  0.7515444
## [3,] -0.5176936  0.8546767  0.2682231

# eigenvalues
spectrum$values

## [1]  2.2469796 -0.8019377  0.5549581

# check 
x = spectrum$vectors[, 1]
lambda = spectrum$values[1]
lambda * x

## [1] -0.9328517 -1.6809406 -1.1632470

C %*% x

##            [,1]
## [1,] -0.9328517
## [2,] -1.6809406
## [3,] -1.1632470

1. create a square matrix
2. verify that the trace of a matrix (the sum of the diagonal elements) is the sum of its eigenvalues
3. verify that the determinant of a matrix (function det) is the product of its eigenvalues (use function prod)

C = matrix(c(1,0,1, 1,1,1, 1,1,0), nrow=3, byrow=TRUE)

(v = eigen(C)$values)

## [1]  2.2469796 -0.8019377  0.5549581

# trace
sum(diag(C))

## [1] 2

sum(v)

## [1] 2

# determinant
det(C)

## [1] -1

prod(v)

## [1] -1

A data frame is a list of vectors (called columns). A data frame is like a database table:

• each column has a name and contain elements of the same type
• different columns have the same length and may have different types

name = c("John", "Samuel", "Uma", "Bruce", "Tim")
age = c(23, 31, 17, 41, 25)
married = c(TRUE, FALSE, FALSE, TRUE, TRUE)

pulp = data.frame(name, age, married)
pulp

##     name age married
## 1   John  23    TRUE
## 2 Samuel  31   FALSE
## 3    Uma  17   FALSE
## 4  Bruce  41    TRUE
## 5    Tim  25    TRUE

# first row
pulp[1, ]

##   name age married
## 1 John  23    TRUE

# first column
# matrix style
pulp[ ,1]

## [1] "John"   "Samuel" "Uma"    "Bruce"  "Tim"

pulp[, "name"]

## [1] "John"   "Samuel" "Uma"    "Bruce"  "Tim"

# list style (remember a data frame is a list)
pulp$name 

## [1] "John"   "Samuel" "Uma"    "Bruce"  "Tim"

pulp[[1]]

## [1] "John"   "Samuel" "Uma"    "Bruce"  "Tim"

# filtering
pulp[pulp$name == "Uma", ]

##   name age married
## 3  Uma  17   FALSE

pulp[pulp$age < 18, ]

##   name age married
## 3  Uma  17   FALSE

pulp[married == TRUE, "name"]

## [1] "John"  "Bruce" "Tim"

Consider:

students = data.frame(
  name = c("Ada", "Linus", "Grace", "Alan", "Barbara", "Dennis"),
  grade = c(28, 17, 30, 22, 19, 26),
  attended = c(TRUE, TRUE, FALSE, TRUE, TRUE, FALSE)
)

Find:

1. all students with grade at least 18;
2. the names of students with grade at least 18 who attended the course;
3. the mean grade of students who attended.

# 1
students[students$grade >= 18, ]

##      name grade attended
## 1     Ada    28     TRUE
## 3   Grace    30    FALSE
## 4    Alan    22     TRUE
## 5 Barbara    19     TRUE
## 6  Dennis    26    FALSE

# 2
students[students$grade >= 18 & students$attended == TRUE, "name"]

## [1] "Ada"     "Alan"    "Barbara"

# 3
mean(students$grade[students$attended == TRUE])

## [1] 21.5

1. Matrices are more efficient than data frames
2. In particular, if all columns have the same type, use a matrix, not a data frame

library(microbenchmark)

# A data frame
df = data.frame(
  a = lubridate::now() + runif(1000) * 86400, # timestamp
  b = lubridate::today() + runif(1000) * 30,  # date
  c = 1:1000, # integer
  d = runif(1000), # double
  e = sample(letters, 1000, replace = TRUE) # character
)

# A matrix
mat = matrix(runif(5000), ncol = 5)

# Get a column
microbenchmark(mat[, 3], df[, 3])

## Unit: microseconds
##      expr   min     lq    mean median     uq    max neval
##  mat[, 3] 2.160 2.4650 3.58374 2.7195 3.6015 23.478   100
##   df[, 3] 5.773 5.9255 8.05321 6.0405 6.3925 62.394   100

# Get a row
microbenchmark(mat[1, ], df[1, ])

## Unit: nanoseconds
##      expr   min      lq  mean  median    uq    max neval
##  mat[1, ]   377   442.5   724   641.5   763   8694   100
##   df[1, ] 45793 46993.0 52307 48166.0 52073 164400   100

Since a data frame is a list, and lists can contain other lists as elements, you can create nested data frames, that is data frames whose elements are data frames.

# a data frame
Venus = data.frame(
  x = c(17, 19), 
  y = c("Hello", "Venus")
)

# a data frame
Jupiter = data.frame(
  x = c(21, 23), 
  y = c("Hello", "Jupiter")
)

# a nested data frame
# I() treats the object ‘as is’
worlds = data.frame(
  x = I(list(Venus, Jupiter)), 
  y = c("Hello", "Worlds")
)

str(worlds)

## 'data.frame':    2 obs. of  2 variables:
##  $ x:List of 2
##   ..$ :'data.frame': 2 obs. of  2 variables:
##   .. ..$ x: num  17 19
##   .. ..$ y: chr  "Hello" "Venus"
##   ..$ :'data.frame': 2 obs. of  2 variables:
##   .. ..$ x: num  21 23
##   .. ..$ y: chr  "Hello" "Jupiter"
##   ..- attr(*, "class")= chr "AsIs"
##  $ y: chr  "Hello" "Worlds"

worlds$x[[1]]

##    x     y
## 1 17 Hello
## 2 19 Venus

worlds$x[[2]]

##    x       y
## 1 21   Hello
## 2 23 Jupiter

R is an Turing-complete (functional) programming language.

It includes conditional statements:

x = 49
if (x %% 7 == 0) x else -x

## [1] 49

And loops:

x = 108
i = 2
while (i <= x/2) {
 if (x %% i == 0) print(i)
 i = i + 1;
}  

## [1] 2
## [1] 3
## [1] 4
## [1] 6
## [1] 9
## [1] 12
## [1] 18
## [1] 27
## [1] 36
## [1] 54

for (i in 2:(x/2)) {
  if (x %% i == 0) print(i)
}

## [1] 2
## [1] 3
## [1] 4
## [1] 6
## [1] 9
## [1] 12
## [1] 18
## [1] 27
## [1] 36
## [1] 54

df = data.frame(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)

# we know the sequence and output lengths
# create a vector of a given size
output = vector("double", ncol(df))  
for (i in 1:ncol(df)) {               
  output[i] = mean(df[[i]])          
}

# we know the sequence length but we do NOT know the output length (slow solution)
means = c(0, 1, 2)
# a vector of doubles of length 0
output = double() 
for (i in 1:length(means)) {
  n = sample(1:100, 1)
  # dynamically increase the vector (slow)
  output = c(output, rnorm(n, means[i])) 
}

# create a list with length(means) elements (faster solution)
output = vector("list", length(means))
for (i in 1:length(means)) {
  n = sample(1:100, 1)
  output[[i]] = rnorm(n, means[i])
}
# unlist the list into a vector
output = unlist(output)  

# we do not know the sequence length
# iterate until a sequence of Heads of length difficulty is found
flips = 0
nheads = 0
difficulty = 10

while (nheads < difficulty) {
  if (sample(c("T", "H"), 1) == "H") {
    nheads = nheads + 1
  } else {
    nheads = 0
  }
  flips = flips + 1
}
flips

## [1] 235

Most of the times you can perform your task by applying functions, avoiding loops. This is typically faster.

x = 1:100

# compute the sum (bad)
s = 0
for (i in 1:length(x)) {               
  s = s + x[i]
}
s

## [1] 5050

# compute the sum (good)
sum(x)

## [1] 5050

# even faster
n = length(x)
n * (n+1) / 2

## [1] 5050

You may use built-in functions:

log

## function (x, base = exp(1))  .Primitive("log")

args(log)

## function (x, base = exp(1)) 
## NULL

log(x = 128, base = 2)

## [1] 7

log(base = 2, x = 128)

## [1] 7

log(128, 2)

## [1] 7

log(2, 128)

## [1] 0.1428571

log(128)

## [1] 4.85203

Or define your our functions:

euclidean = function(x=0, y=0) {sqrt(x^2 + y^2)}

euclidean(1, 1)

## [1] 1.414214

euclidean(1)

## [1] 1

euclidean()

## [1] 0

Or define your own binary operators using functions:

# xor
'%()%' = function(x, y) {(x | y) & !(x & y)}

TRUE %()% TRUE 

## [1] FALSE

TRUE %()% FALSE

## [1] TRUE

FALSE %()% TRUE

## [1] TRUE

FALSE %()% FALSE

## [1] FALSE

Functions may be recursive:

factorial = function(x) {
 if (x == 0) 1 else x * factorial(x-1)
}
factorial(5)

## [1] 120

You may write functionals, that are functions whose arguments are other functions:

# compute the sum of applications of f up to n
g = function(f, n) {
 sum = 0
 for (i in 1:n) sum = sum + f(i)
 return(sum)
}
 
g(factorial, 5)

## [1] 153

An application of functionals and iteration is the set of apply-like functionals:

df = data.frame(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)

# apply mean to each column of data frame, returns a list
lapply(df, mean)

## $a
## [1] -0.3489407
## 
## $b
## [1] -0.1119336
## 
## $c
## [1] 0.1162174
## 
## $d
## [1] 0.3462086

# apply mean to each column of data frame, returns an atomic vector
sapply(df, mean)

##          a          b          c          d 
## -0.3489407 -0.1119336  0.1162174  0.3462086

# apply to a vector
sapply(1:100, function(x) {x^2})

##   [1]     1     4     9    16    25    36    49    64    81   100   121   144
##  [13]   169   196   225   256   289   324   361   400   441   484   529   576
##  [25]   625   676   729   784   841   900   961  1024  1089  1156  1225  1296
##  [37]  1369  1444  1521  1600  1681  1764  1849  1936  2025  2116  2209  2304
##  [49]  2401  2500  2601  2704  2809  2916  3025  3136  3249  3364  3481  3600
##  [61]  3721  3844  3969  4096  4225  4356  4489  4624  4761  4900  5041  5184
##  [73]  5329  5476  5625  5776  5929  6084  6241  6400  6561  6724  6889  7056
##  [85]  7225  7396  7569  7744  7921  8100  8281  8464  8649  8836  9025  9216
##  [97]  9409  9604  9801 10000

mtx <- cbind(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)


# apply mean to each column of matrix, returns an atomic vector
apply(mtx, 2, mean)

##          a          b          c          d 
## -0.1787799 -0.3330621  0.2865619  0.4891630

# apply mean to each row of matrix, returns an atomic vector
apply(mtx, 1, mean)

##  [1] -0.06644800  0.66865979  0.26506432 -0.33323834  0.07575792 -0.08472629
##  [7] -0.83181827 -0.26466914  0.57091241  0.66021277

Write a function num_divisors(n) that returns the number of positive divisors of n.

For example:

num_divisors(6) = 4

because the divisors are 1, 2, 3, 6.

Then:

1. compute the number of divisors of every integer from 1 to 30 using sapply;
2. use the result to find all prime numbers between 1 and 30.

Try not to use an explicit loop.

num_divisors = function(n) {
  sum(n %% (1:n) == 0)
}

Test:

num_divisors(6)

## [1] 4

Apply it:

x = 1:30

d = sapply(x, num_divisors)
d

##  [1] 1 2 2 3 2 4 2 4 3 4 2 6 2 4 4 5 2 6 2 6 4 4 2 8 3 4 4 6 2 8

Prime numbers have exactly two positive divisors:

x[d == 2]

##  [1]  2  3  5  7 11 13 17 19 23 29

We repeatedly toss two fair coins.

We want to estimate:

\[ P(\text{both coins are heads}\mid\text{at least one coin is heads}). \]

1. Guess the theoretical probabilty of the event;
2. Write an R simulation using 100,000 experiments and compare the simulation with the exact theoretical probability.

The probability of the event is \(1/3\), since we are looking for the event HH over possible events are HH, HT, TH.

The simulation is the following:

n = 10^6

coin1 = sample(c("H", "T"), n, replace = TRUE)
coin2 = sample(c("H", "T"), n, replace = TRUE)

cond = (coin1 == "H" | coin2 == "H")
prob = (coin1[cond] == "H" & coin2[cond] == "H")

mean(prob)

## [1] 0.3340087

Write a iterative function that, given a square matrix \(A\) and an integer \(n \geq 0\), computes the power \(A^n\) (use diag function to build the identity matrix).

power = function(A, n) {
  I = diag(nrow(A))
  if (n == 0) return(I)
  if (n == 1) return(A)
  B = A
  for (i in 2:n) {
    B = B %*% A 
  }
  return(B)
}

A = matrix(c(1,0,0, 0,2,0, 0,0,3), nrow=3, byrow=TRUE)
power(A, 5)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0   32    0
## [3,]    0    0  243

The previous material computes \(A^n\) by repeatedly multiplying by \(A\).

Write instead a recursive function:

fast_power(A, n)

using the rules

\[ A^0 = I \]

and, when \(n\) is even,

\[ A^n = \left(A^{n/2}\right)^2, \]

while when \(n\) is odd,

\[ A^n = A A^{n-1}. \]

Test it using:

A = matrix(
  c(1, 0, 0,
    0, 2, 0,
    0, 0, 3),
  nrow = 3,
  byrow = TRUE
)

and compute \(A^{10}\).

Finally, benchmark the two functions power and fast_power using the microbenchmark package.

fast_power = function(A, n) {

  if (n == 0) {
    return(diag(nrow(A)))
  }

  if (n %% 2 == 0) {

    B = fast_power(A, n %/% 2)
    return(B %*% B)

  } else {

    return(A %*% fast_power(A, n - 1))
  }
}

A = matrix(
  c(1, 0, 0,
    0, 2, 0,
    0, 0, 3),
  nrow = 3,
  byrow = TRUE
)

fast_power(A, 10)

##      [,1] [,2]  [,3]
## [1,]    1    0     0
## [2,]    0 1024     0
## [3,]    0    0 59049

A = matrix(runif(100), nrow = 10)
microbenchmark(power(A, 100), fast_power(A, 100), times = 1000)

## Unit: microseconds
##                expr     min      lq     mean   median       uq       max neval
##       power(A, 100) 110.600 145.225 160.9471 149.2175 164.6385   645.693  1000
##  fast_power(A, 100)  18.492  22.295  45.6738  22.7830  23.8895 20370.427  1000

# a data matrix
M = matrix(c(
   c(1200, 1190, 1100, 1120, 890),
   c(6200, 6690, 6700, 7120, 7150),
   c(8900, 8790, 8760, 8800, 9010),
   c(3300, 3490, 3660, 4300, 4510),
   c(2190, 2000, 1890, 1740, 1500)), ncol = 5
)  

# give names to rows
rownames(M) = 2014:2018
# give names to columns
colnames(M) = LETTERS[1:5]
M

##         A    B    C    D    E
## 2014 1200 6200 8900 3300 2190
## 2015 1190 6690 8790 3490 2000
## 2016 1100 6700 8760 3660 1890
## 2017 1120 7120 8800 4300 1740
## 2018  890 7150 9010 4510 1500

# barplot
barplot(M[1,])

# stacked barplot
barplot(M, legend=TRUE)

#  juxtaposed barplot
barplot(M, beside=TRUE, legend=TRUE)

# histogram
x = rnorm(1000)
hist(x, probability=TRUE, main="Histogram of a normal sample")
## add distribution
rug(x)

# density plot
plot(density(x), main="Density of a normal sample")
rug(x)

# boxplot
# If range is positive, the whiskers extend to the most extreme 
# data point  which is no more than range times the interquartile 
# range from the box.  A value of zero causes the whiskers to extend 
# to the data extremes.
x = rnorm(1000)
boxplot(x, range = 1.5)

boxplot(x, range = 0)

# scatter plot
x = rnorm(100)
y = rnorm(100)
plot(x, y)

y = x + runif(100)
plot(x, y)

• R Markdown provides an unified authoring framework for data science, combining your code, its results, and your prose commentary
• R Markdown documents are fully reproducible and support many output formats, like HTML, PDF, and slideshows

R Markdown files are designed to be used in three ways:

1. for communicating to decision makers, who want to focus on the conclusions, not the code behind the analysis
2. for collaborating with other data scientists (including future you!), who are interested in both your conclusions, and how you reached them (the code)
3. as an environment in which to do data science, as a modern day lab notebook where you can capture not only what you did, but also what you were thinking

In an R Markdown (.Rmd) document, the section enclosed between two lines of --- at the beginning is called the YAML header (or YAML front matter).

It specifies metadata about the document and controls how the document is rendered.

Read about YAML parameters.

In an R Markdown document, knitr::opts_chunk$set() defines the default options for all R code chunks. Individual chunks can override these defaults.

Here are the most useful options:

Important distinctions:

• eval = FALSE: code is displayed but not executed.
• echo = FALSE: code is executed but not displayed.
• include = FALSE: code is executed, but neither code nor results are displayed.

Create an R Markdown document investigating the following simulated relationship:

x = rnorm(500)
y = x + runif(500)

Your report must contain:

1. a short explanation of the question being investigated;
2. creation of a data frame containing x and y;
3. the mean of both variables;
4. a histogram of x
5. a boxplot of y
6. a scatter plot of x against y;
7. a short written interpretation of what the plots show;
8. both code and results in the knitted output.

Here is an R Markdown solution.

• R Cookbook