R Programming/Mathematics

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 , 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 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 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. Probability 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 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 . 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 MaximaMaxima is open source http://maxima.sourceforge.net/, SAGESAGE is an open source package which includes R and Maxima : http://www.sagemath.org/, MathematicaMathematica is not open source http://www.wolfram.com/products/mathematica/index.html See also The following command gives help on special mathematical functions related to the beta and gamma functions. ?Special References 2 de:GNU R: Rechnen mit R pt:R (linguagem de programação)/Matemática