# R Programming/Mathematics

## 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 <- matrix(1:3,1,3)
> u <- matrix(1:3,3,1)
> u%*%v
[,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    4    6
[3,]    3    6    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
>
[,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).

## 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 $a x^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$.

### 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 $-\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)
> 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


## ProbabilityEdit

• 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


## 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 $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 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]

?Special