Open main menu

BasicsEdit

?Arithmetic
?Special

Linear AlgebraEdit

VectorsEdit

The inner productEdit

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 productEdit

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 AlgebraEdit

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 matrixEdit

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 calculationsEdit

> 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 transpositionEdit

  • 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 matrixEdit

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

Matrix inversionEdit

  • 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 equationEdit

> 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 eigenspaceEdit

  • 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

MiscEdit

  • 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.

AnalysisEdit

Logarithm and ExponentsEdit

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 equationsEdit

To solve  , where   are given numbers, use the command

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

So, for example, to calculate the roots of the equation   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  .

See also polynom and multipol packages

DerivativesEdit

Symbolic calculationsEdit

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 approximationEdit

  • numDeriv package

IntegrationEdit

R can perform one dimensional integration. For example we can integrate over the density of the normal distribution between   and  

> 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.

ProbabilityEdit

  • The number of combination of length k within n numbers :  
> 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

ArithmeticsEdit

The factorial functionEdit

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  . 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 divisionEdit

  • 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

GeometryEdit

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

Symbolic calculusEdit

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

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

See alsoEdit

The following command gives help on special mathematical functions related to the beta and gamma functions.

?Special

ReferencesEdit

  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