MATH AND SIMULATIONS IN R

Julius Schmid

R contains built-in functions for your favorite math operations and, of course, for statistical distributions.

For example, we can create a scatterplot for the built-in cars data set just by calling plot(cars):

plot(cars)

Since the cars data set has only two variables speed and dist, we do not need to specify which columns we use to create the scatterplot. There is no other possibility to create a scatterplot with the cars data at hand (apart from switching the axes).

Extended Example: Calculating a Probability

As our first example, we’ll work through calculating a probability using the prod() function. Suppose we have n independent events, and the i-th event has the probability p(i) of occurring. What is the probability of exactly one of these events occurring? Suppose first that n = 3 and our events are named A, B, and C. Then we break down the computation as follows:

P(exactly one event occurs) = P(A and not B and not C) + P(not A and B and not C) + P(not A and not B and C)

Here’s code to compute this, with our probabilities p(i) contained in the vector p:

exactlyone <- function(p) {
notp <- 1 - p
tot <- 0.0
for (i in 1:length(p))
tot <- tot + p[i] * prod(notp[-i])
return(tot)
}

How does it work? Well, the assignment notp <- 1 - p creates a vector of all the “not occur” probabilities 1 − p , using recycling. The expression notp[-i] computes the product of all the elements of notp, except the ith—exactly what we need.

Cumulative Sums and Products

As mentioned, the functions cumsum() and cumprod() return cumulative sums and products.

x <- c(2,7,13,4)
cumsum(x)

[1]  2  9 22 26

In x, the sum of the first element is 2, the sum of the first two elements is 9, the sum of the first three elements is 22, and the sum of the first four elements (all elements) is 26.

Now, apply the cumprod() function:

cumprod(x)

[1]   2  14 182 728

The function cumprod() works the same way as cumsum(), but with the product instead of the sum. The output is a four-dimensional vector. The first entry of the output vector is just the first entry of x (which is 2). The product of the first two numbers 14. the product of the first three numbers is 182, and the product of the first four numbers (all numbers) is 728.

Minima and Maxima

---

There is quite a difference between min() and pmin(). The former simply combines all its arguments into one long vector and returns the minimum value in that vector. In contrast, if pmin() is applied to two or more vectors, it returns a vector of the pair-wise minima, hence the name pmin.

Here’s an example:

z <- matrix(c(3,4,6,7,3,2), nrow = 2)
z

     [,1] [,2] [,3]
[1,]    3    6    3
[2,]    4    7    2

When applying the min() function, we determine the global minimum over all matrix entries:

min(z[,1],z[,2],z[,3])

[1] 2

The minimum value is 2 which belongs to the entry z[3,2].

Now, let us return the pairwise minima for each row which we get by applying the pmin() function:

pmin(z[,1],z[,2],z[,3])

[1] 3 2

The minimum in the first row is 3 (in z[1,1]), and the minimum value in the second row is 2 (belonging to z[3,2]).

In the first case, min() computed the smallest value in (3,4,6,7,3,2). But the call to pmin() computed the smaller of 3, 6, and 3, yielding 3; and then the smaller of 4, 7, and 2, which is 2. Thus, the call returned the vector (3,2).

The max() and pmax() functions act analogously to min() and pmin(). Function minimization/maximization can be done via nlm() and optim(). For example, let’s find the smallest value of f(x) = x^2 − sin(x).

Let us use the nlm() function to solve this where nlm is an abbreviation for nonlinear minimization. nlm() takes the function x^2 - sin(x) as an input together with 8 as a starting parameter for the optimization:

nlm(function(x) return(x^2-sin(x)),8)

$minimum
[1] -0.2324656

$estimate
[1] 0.4501831

$gradient
[1] 4.024558e-09

$code
[1] 1

$iterations
[1] 5

nlm() then returns the outputs minimum (y-value), estimate (argmin, x-value), a gradient (which is 0 for a minimum), a code value and the number of iterations in order to reach the minimum. In this case, it took the algorithm only 5 iteration steps to get to the result.

The minimum value was found to be approximately −0.23, occurring at x = 0.45. A Newton-Raphson method (a technique from numerical analysis for approximating roots) is used, running five iterations in this case. The second argument specifies the initial guess, which we set to be 8. (This second argument was picked pretty arbitrarily here, but in some problems, you may need to experiment to find a value that will lead to convergence.)

Recreate a similar case scenario shown above but this time using the function f(x) = 3x^2 − cos(x)

Again, we use the nlm() function in order to find the minimum. Information about the nlm() function can be found in the comments above. This time, our first input is the function 3x^2 - cos(x), and our second input is the starting value 3.

#Enter Answer Here
nlm(function(x) return(3*x^2-cos(x)),3)

$minimum
[1] -1

$estimate
[1] -5.132311e-07

$gradient
[1] -9.25926e-08

$code
[1] 1

$iterations
[1] 3

The minimum value was found to be − 1, occurring at x = 0. A Newton-Raphson method (a technique from numerical analysis for approximating roots) is used, running three iterations in this case. The gradient is 0 (nlm() returns almost 0 because of floating-point error) which is just what we would expect for an extreme point of a function.

Calculus

---

R also has some calculus capabilities, including symbolic differentiation and numerical integration, as you can see in the following example.

In order to calculate the first derivative of a function, we apply the D() function in R, where we need to state the function as well as the variable that we want to derive.

D(expression(exp(x^2)),"x") # derivative

exp(x^2) * (2 * x)

We find that the derivative for x of exp(x^2)) is given by exp(x^2) * (2x).

Find the derivative of exp(x^3)+x

We proceed the same way as we did above, applying the D() function, and setting the function and the variable that we want to derive:

#Enter Answer Here
D(expression(exp(x^3) + x),"x")

exp(x^3) * (3 * x^2) + 1

We find that the derivative for x of exp(x^3) + x is given by exp(x^3) * (3x^2) + 1.

We can also integrate functions, calling integrate() in R. Here, we need to sepcify the function, as well as the lower and upper bounds for the integral. For example, if we want to integrate x^2 between 0 and 1, we enter the following code:

integrate(function(x) x^2,0,1)

0.3333333 with absolute error < 3.7e-15

The integral of x^2 within the boundaries 0 and 1 is evaluated to 1/3.

Next, let us integrate 2sin(x) + 3 between 0 and pi, we enter the following code:

#Enter Answer Here
integrate(function(x) 2 * sin(x) + 3 ,0,pi)

13.42478 with absolute error < 1.5e-13

The integral of 2sin(x) + 3 within the boundaries 0 and pi is evaluated to be approx. 13.42.

Set Operations

---

R includes some handy set operations, including these: • union(x,y): Union of the sets x and y • intersect(x,y): Intersection of the sets x and y • setdiff(x,y): Set difference between x and y, consisting of all elements of x that are not in y • setequal(x,y): Test for equality between x and y • c %in% y: Membership, testing whether c is an element of the set y • choose(n,k): Number of possible subsets of size k chosen from a set of size n

Let us create two variables x and y to which we can apply the set functions. Note that x and y are not interpreted as vectors, but rather as sets when we apply the functions above to them.

x <- c(1,2,5)
y <- c(5,1,8,9)

First, we build the union of the sets x and y.

union(x,y)

[1] 1 2 5 8 9

We get a set that collects all values of both sets x and y, without duplicates.

Next, we want to see if x and y have values in common. Therefore, we use the intersect() function:

intersect(x,y)

[1] 1 5

The numbers 1 and 5 are contained in both sets x and y.

Now we want to evaluate the mathematical expression x/y, returning all values in x that are not contained in y:

setdiff(x,y)

[1] 2

Only the number 2 is part of x and not contained in y.

We do this the other way around, returning all numbers y that are not contained in x: (y/x)

setdiff(y,x)

[1] 8 9

We get two numbers 8 and 9, that are part of the set y and not part of x.

Create two vectors and x1 and y1 and simulate a similar case scenario.

x1 includes every second number up to 12, and y1 includes every third number up to 12.

#Enter Answer Here
x1 <- c(2,4,6,8,10,12)
y1 <- c(3,6,9,12)

First, we build the union of the sets x1 and y1.

union(x1,y1)

[1]  2  4  6  8 10 12  3  9

We get a set that collects all values of both sets x1 and y1, without duplicates. These are all values that are divisible by two or three, up to twelve.

Next, we want to see if x1 and y1 have values in common. Therefore, we use the intersect() function:

intersect(x1,y1)

[1]  6 12

The numbers 6 and 12 are contained in both sets x1 and y1. They are divisible by both 2 and 3.

Now we want to evaluate the mathematical expression x1/y1, returning all values in x1 that are not contained in y1:

setdiff(x1,y1)

[1]  2  4  8 10

The numbers 2, 4, 8, and 10 are divisible by 2, but not by 3.

We do this the other way around, returning all numbers y1 that are not contained in x1: (y1/x1)

setdiff(y1,x1)

[1] 3 9

The numbers 3 and 9 are divisible by 3, but not by 2. We conclude that half of the set y1 is intersected with x1 (6 and 12), the other half is not (3 and 9).