R Programming/Mathematics

R Basics
Data Management
Graphics
- Grammar of graphics
Publication quality output
Descriptive Statistics
Mathematics
Statistical Core Methods
Regression Models
Time Series
Factor Analysis
Classification
- Ordination
- Clustering
Network Analysis
High Performance Computing
- Profiling R code
- Parallel computing with R
Appendix
- Sources
- Index

Basics

?Arithmetic
?Special

Linear Algebra

Vectors

The inner product

The inner product is also called the dot product or the scalar product. It is the sum of the item-by-item product.

> u <- rep(3,3)
> v <- 1:3
> u%*%v # the inner product
     [,1]
[1,]   18

The outer product

The outer product is also called the cross product or the vector product. It is a matrix resulting from the product of the elements of the two vectors.

> v <- rep(3,3)
> u <- 1:3
> u%o%v # The outer product
     [,1] [,2] [,3]
[1,]    3    3    3
[2,]    6    6    6
[3,]    9    9    9

Matrix Algebra

If you want to create a new matrix, one way is to use the matrix() function. You have to enter a vector of data, the number of rows and/or columns and finally you can specify if you want R to read your vector by row or by column (the default option) with byrow. You can also combine vectors using cbind() or rbind(). The dimension of a matrix can be obtained using the dim() function or alternatively nrow() and ncol().

> matrix(data = NA, nrow = 5, ncol = 5, byrow = T)
> matrix(data = 1:15, nrow = 5, ncol = 5, byrow = T)
> v1 <- 1:5
> v2 <- 5:1
> cbind(v1,v2)
> rbind(v1,v2)
> dim(X)
> nrow(X)
> ncol(X)

Some special matrix

The identity matrix has ones on the diagonal and zeros outside the diagonal.

eye() (matlab)
diag(1,nrow=10,ncol=10)
diag(rep(1,10))

J matrix is full of ones

ones() (matlab)

A matrix full of zeros

zeros() (matlab)

> library(matlab)
> eye(3)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
> ones(3)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1
> zeros(3) 
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

Diagonal matrix

> diag(3)

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

Upper triangular

> round(upper.tri(matrix(1, n, n))) 

for n=3
     [,1] [,2] [,3]
[1,]    0    1    1
[2,]    0    0    1
[3,]    0    0    0

If you also need the diagonal of one's 

> round(upper.tri(matrix(1, 3, 3), diag = TRUE))

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

Lower triangular

Same as upper triangular but using lower.tri instead

create an Hilbert matrix using hilbert() (fUtilities).

Matrix calculations

compute a matrix multiplication X%*%Y.

> b <- matrix(nrow = 2, ncol = 2, c(1, 2, 3, 4))
> a <- matrix(nrow = 2, ncol = 2, c(1, 0, 0, -1))
> a
     [,1] [,2]
[1,]    1    0
[2,]    0   -1
> b
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> a%*%b
     [,1] [,2]
[1,]    1    3
[2,]   -2   -4
> b%*%a
     [,1] [,2]
[1,]    1   -3
[2,]    2   -4

compute the Kronecker product using %x% or kron() (fUtilities).

> M <- matrix(rep(2,4),nrow = 2) 
> M
     [,1] [,2]
[1,]    2    2
[2,]    2    2
> I <- eye(2) 
> I
     [,1] [,2]
[1,]    1    0
[2,]    0    1
> I %x% M 
     [,1] [,2] [,3] [,4]
[1,]    2    2    0    0
[2,]    2    2    0    0
[3,]    0    0    2    2
[4,]    0    0    2    2
> library(fUtilities)
> kron(I,M)
     [,1] [,2] [,3] [,4]
[1,]    2    2    0    0
[2,]    2    2    0    0
[3,]    0    0    2    2
[4,]    0    0    2    2

Matrix transposition

Transpose the matrix

> t(M)
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    1    2
[3,]    0    0    1

The trace and determinant of a matrix

compute the trace of a matrix using tr() (fUtilities)
returns the rank of a matrix using rk() (fBasics:)

Matrix inversion

Invert a matrix using solve() or inv() (fUtilities). We can also compute the generalized inverse using ginv() in the MASS package.

> M <- cbind(c(1,0,1),c(0,1,2),c(0,0,1))
> solve(M)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]   -1   -2    1
> solve(M)%*%M
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

Solving a linear equation

> m=matrix(nrow=2,ncol=2,c(1,-.8,1,.2))
> m
     [,1] [,2]
[1,]  1.0  1.0
[2,] -0.8  0.2
> 
> l=matrix(c(1.0+25.0/18,25.0/18.0))
> l
         [,1]
[1,] 2.388889
[2,] 1.388889
> 
> k=solve(m,l)
> k
           [,1]
[1,] -0.9111111
[2,]  3.3000000
> 
> m%*%k          #checking the answer
         [,1]
[1,] 2.388889
[2,] 1.388889
>

Eigenvalue, eigenvector and eigenspace

Eigenvalues and eigenvectors

> eigen(M)
$values
[1] 1 1 1

$vectors
     [,1]          [,2]          [,3]
[1,]    0  2.220446e-16  0.000000e+00
[2,]    0  0.000000e+00  1.110223e-16
[3,]    1 -1.000000e+00 -1.000000e+00

Misc

compute the norm of a matrix using norm() (fUtilities).
check if a matrix is positive definite isPositiveDefinite() (fUtilities).
make a matrix positive definite makePositiveDefinite() (fUtilities).

computes row statistics and column statistics (fUtilities).
extract the upper and the lower part of a matrix triang() and Triang() (fUtilities).
See also the matrix, matlab, matrixcalc, matrixStats packages.

Analysis

Logarithm and Exponents

We have the power function 10^3 or 10**3 , the logarithm and the exponential log(2.71), log10(10),exp(1).

> 10^3 # exponent
[1] 1000
> 10**3 # exponent
[1] 1000
> exp(1) # exponential
[1] 2.718282
> log(2.71) # natural logarithm
[1] 0.9969486
> log10(1000) # base 10 logarithm
[1] 3
> log(1000,base = 10) # base 10 logarithm
[1] 3

Polynomial equations

To solve $ax^{k}+bx^{k-1}+\cdots +n=0$ , where $a,b,\dots ,n$ are given numbers, use the command

> polyroot(c(n,...,b,a))

So, for example, to calculate the roots of the equation $2x^{2}-5x-3=0$ one would do as follows:

> polyroot(c(-3,-5,2))
 [1] -0.5+0i  3.0-0i

and the solution can be read to be $x_{1}=-0.5\land x_{2}=3$ .

See also polynom and multipol packages

Derivatives

Symbolic calculations

R can give the derivative of an expression. You need to convert your function as an expression using the expression() function. Otherwise you get an error message.

Here are some examples :

> D(expression(x^n),"x")
x^(n - 1) * n
> D(expression(exp(a*x)),"x")
exp(a * x) * a
> D(expression(1/x),"x")
-(1/x^2)
> D(expression(x^3),"x")
3 * x^2
> D(expression(pnorm(x)),"x")
dnorm(x)
> D(expression(dnorm(x)),"x")
-(x * dnorm(x))

Numerical approximation

numDeriv package

Integration

R can perform one dimensional integration. For example we can integrate over the density of the normal distribution between $-\infty$ and $+\infty$

> integrate(dnorm,-Inf,Inf)
1 with absolute error < 9.4e-05
> integrate(dnorm,-1.96,1.96)
0.9500042 with absolute error < 1.0e-11
> integrate(dnorm,-1.64,1.64)
0.8989948 with absolute error < 6.8e-14
# we can also store the result in an object
> ci90 <- integrate(dnorm,-1.64,1.64)
> ci90$value
[1] 0.8989948
> integrate(dnorm,-1.64,1.64)$value
[1] 0.8989948

see the adapt package for multivariate integration.

> library(adapt)
> ?adapt
> ir2pi <- 1/sqrt(2*pi)
> fred <- function(z) { ir2pi^length(z) * exp(-0.5 * sum(z * z))}
> 
> adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred)
       value       relerr       minpts       lenwrk        ifail 
    1.039222 0.0007911264          231           73            0 
> adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred, eps = 1e-4)
       value       relerr       minpts       lenwrk        ifail 
    1.000237 1.653498e-05          655          143            0 
> adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred, eps = 1e-6)
      value      relerr      minpts      lenwrk       ifail 
   1.000039 3.22439e-07        1719         283           0

See also integrate.gh() in the ecoreg package.

Probability

The number of combination of length k within n numbers : $C_{n}^{k}$

> choose(100, 5)
[1] 75287520

Union and intersection

> union(1:10, 5:7)
[1]  1  2  3  4  5  6  7  8  9 10
> intersect(1:10, 5:7)
[1] 5 6 7

Arithmetics

The factorial function

factorial returns the factorial of an integer. This can also be computed using the prod() (product) applied to the vector of integers between 1 and the number of interest.

> factorial(3)
[1] 6
> prod(1:3)
[1] 6

Note that by convention $0!=1$ . factorial() returns 1 in 0. This is not the case with the prod() functions.

> factorial(0)
[1] 1
> prod(0)
[1] 0

Factorial numbers can be very large and cannot be computed for high values.

> factorial(170)
[1] 7.257416e+306
> factorial(171)
[1] Inf
Message d'avis :
In factorial(171) : value out of range in 'gammafn'

The modulo function and euclidian division

Modulo and integer division (i.e. euclidean division)

> 5%%2
[1] 1
>5%/%2
[1] 2

Note: R is affected by the problem with non integer numbers and euclidian divisions.

> .5%/%.1 # we get 4 instead of 5
[1] 4
> .5%%.1 # we get .1 instead of 0
[1] 0.1

Geometry

pi the constant
cos(), sin(), tan() the trigonometric functions.

Symbolic calculus

rSymPy (rsympy) provides sympy (link) functions in R.

If you want to do more symbolic calculus, see Maxima^[1], SAGE^[2], Mathematica^[3]

References

↑ Maxima is open source http://maxima.sourceforge.net/
↑ SAGE is an open source package which includes R and Maxima : http://www.sagemath.org/
↑ Mathematica is not open source http://www.wolfram.com/products/mathematica/index.html

[1] Maxima is open source http://maxima.sourceforge.net/

[2] SAGE is an open source package which includes R and Maxima : http://www.sagemath.org/

[3] Mathematica is not open source http://www.wolfram.com/products/mathematica/index.html

[1]

[2]

[3]