Przeglądaj wersję html pliku:

Matlab R Reference (ang)


MATLAB R / R Reference
March 3, 2009 David Hiebeler Dept. of Mathematics and Statistics University of Maine Orono, ME 04469-5752 http://www.math.umaine.edu/faculty/hiebeler

I wrote the first version of this reference during the Spring 2007 semester, as I learned R while teaching my course “MAT400, Modeling & Simulation” at the University of Maine. The course covers population and epidemiological modeling, including deterministic and stochastic models in discrete and continuous time, along with spatial models. Half of the class meetings are in a regular classroom, and half are in a computer lab where students work through modeling & simulation exercises. When I taught earlier versions of the course, it was based on Matlab only. In Spring 2007, some biology graduate students in the class who had learned R in statistics courses asked if they could use R in my class as well, and I said yes. My colleague Bill Halteman was a great help as I frantically learned R to stay ahead of the class. As I went, every time I learned how to do something in R for the course, I added it to this reference, so that I wouldn’t forget it later. Some items took a huge amount of time searching for a simple way to do what I wanted, but at the end of the semester, I was pleasantly surprised that almost everything I do in Matlab had an equivalent in R. I was also inspired to do this after seeing the “R for Octave Users” reference written by Robin Hankin. I’ve continued to add to the document, with many additions based on topics that came up while teaching courses on Advanced Linear Algebra and Numerical Analysis. This reference is organized into general categories. There is also a Matlab index and an R index at the end, which should make it easy to look up a command you know in one of the languages and learn how to do it in the other (or if you’re trying to read code in whichever language is unfamiliar to you, allow you to translate back to the one you are more familiar with). The index entries refer to the item numbers in the first column of the reference document, rather than page numbers. Any corrections, suggested improvements, or even just notification that the reference has been useful will be appreciated. I hope all the time I spent on this will prove useful for others in addition to myself and my students. Note that sometimes I don’t necessarily do things in what you may consider the “best” way in a particular language; I often tried to do things in a similar way in both languages. But if you believe you have a “better” way (either simpler, or more computationally efficient) to do something, feel free to let me know. Acknowledgements: Thanks to Alan Cobo-Lewis and Isaac Michaud for correcting some errors; and Stephen Eglen, David Khabie-Zeitoune, Lee Pang, Manas A. Pathak, and Corey Yanofsky for contributions.

Permission is granted to make and distribute verbatim copies of this manual provided this permission notice is preserved on all copies. Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one. Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions, except that this permission notice may be stated in a translation approved by the Free Software Foundation. Copyright c 2007–2009 David Hiebeler

1

D. Hiebeler, Matlab / R Reference

2

Contents
1 Online help 2 Entering/building/indexing matrices 2.1 Cell arrays and lists . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Structs and data frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Computations 3.1 Basic computations . . . . . . . . . . 3.2 Complex numbers . . . . . . . . . . 3.3 Matrix/vector computations . . . . . 3.4 Root-finding . . . . . . . . . . . . . . 3.5 Function optimization/minimization 3.6 Numerical integration / quadrature . 3.7 Curve fitting . . . . . . . . . . . . . 4 Conditionals, control structure, loops 5 Functions, ODEs 6 Probability and random values 7 Graphics 7.1 Various types of plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Printing/saving graphics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Animating cellular automata / lattice simulations . . . . . . . . . . . . . . . . . . . . . . . 8 Working with files 9 Miscellaneous 9.1 Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2 Strings and Misc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Spatial Modeling Index of MATLAB commands and concepts Index of R commands and concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 4 6 6 7 7 7 8 13 14 14 15 16 19 21 25 25 32 33 34 35 35 36 39 40 44

D. Hiebeler, Matlab / R Reference

3

1
No. 1 2 3

Online help
Description Show help for a function (e.g. sqrt) Show help for a built-in keyword (e.g. for) General list of many help topics Matlab help sqrt, or helpwin sqrt to see it in a separate window help for help R help(sqrt) or ?sqrt help(’for’) or ?’for’ library() to see available libraries, or library(help=’base’) for very long list of stuff in base package which you can see help for help.start()

4

Explore main documentation in browser Search documentation for keyword or partial keyword (e.g. functions which refer to “binomial”)

5

doc or helpbrowser (previously it was helpdesk, which is now being phased out) lookfor binomial

help.search(’binomial’)

D. Hiebeler, Matlab / R Reference

4

2
No. 6

Entering/building/indexing matrices
Description Enter a row vector v 1 2 3 4 = Matlab v=[1 2 3 4] R v=c(1,2,3,4) or alternatively v=scan() then enter “1 2 3 4” and press Enter twice (the blank line terminates input) c(1,2,3,4) (R does not distinguish between row and column vectors.) [1 2 3 ; 4 5 6] To enter values by row: matrix(c(1,2,3,4,5,6), nrow=2, byrow=TRUE) To enter values by column: matrix(c(1,4,2,5,3,6), nrow=2) v[3] A[2,3] A[5]

7

 1  2  Enter a column vector    3  4 1 2 4 5 3 6



[1; 2; 3; 4]

8

Enter a matrix

9 10 11

12 13 14 15

16 17

18 19 20 21

Access an element of vector v Access an element of matrix A Access an element of matrix A using a single index: indices count down the first column, then down the second column, etc. Build the vector [2 3 4 5 6 7] Build the vector [7 6 5 4 3 2] Build the vector [2 5 8 11 14] Build a vector containing n equally-spaced values between a and b inclusive Build a vector of length k containing all zeros Build a vector of length k containing the value j in all positions Build an m×n matrix of zeros Build an m × n matrix containing j in all positions n × n identity matrix In Build diagonal matrix A using elements of vector v as diagonal entries Extract diagonal elements of matrix A “Glue” two matrices a1 and a2 (with the same number of rows) side-by-side “Stack” two matrices a1 and a2 (with the same number of columns) on top of each other

v(3) A(2,3) A(5)

2:7 7:-1:2 2:3:14 linspace(a,b,n)

2:7 7:2 seq(2,14,3) seq(a,b,length.out=n) seq(a,b,len=n) rep(0,k) rep(j,k)

or

just

zeros(k,1) (for a column vector) or zeros(1,k) (for a row vector) j*ones(k,1) (for a column vector) or j*ones(1,k) (for a row vector) zeros(m,n) j*ones(m,n) eye(n) diag(v)

22 23

v=diag(A) [a1 a2]

matrix(0,nrow=m,ncol=n) or just matrix(0,m,n) matrix(j,nrow=m,ncol=n) or just matrix(j,m,n) diag(n) diag(v,nrow=length(v)) (Note: if you are sure the length of vector v is 2 or more, you can simply say diag(v).) v=diag(A) cbind(a1,a2)

24

[a1; a2]

rbind(a1,a2)

D. Hiebeler, Matlab / R Reference No. 25 26 Description Reverse the order of elements in vector v Column 2 of matrix A Matlab v(end:-1:1) A(:,2) R rev(v)

5

27

Row 7 of matrix A

A(7,:)

28 29 30

31

32

33

34 35

36 37 38

All elements of A as a vector, column-by-column Rows 2–4, columns 6–10 of A (this is a 3 × 5 matrix) A 3 × 2 matrix consisting of rows 7, 7, and 6 and columns 2 and 1 of A (in that order) Given a single index ind into an m × n matrix A, compute the row r and column c of that position (also works if ind is a vector) Given the row r and column c of an element of an m × n matrix A, compute the single index ind which can be used to access that element of A (also works if r and c are vectors) Given equal-sized vectors r and c (each of length k), set elements in rows (given by r) and columns (given by c) of matrix A equal to 12. That is, k elements of A will be modified. Truncate vector v, keeping only the first 10 elements Reshape matrix A, making it an m × n matrix with elements taken columnwise from the original A (which must have mn elements) Extract the lower-triangular portion of matrix A Extract the upper-triangular portion of matrix A Enter n × n Hilbert matrix H where Hij = 1/(i + j − 1) Enter an n-dimensional array, e.g. a 3 × 4 × 2 array with the values 1 through 24

A(:) (gives a column vector) A(2:4,6:10) A([7 7 6], [2 1])

A[,2] Note: that gives the result as a vector. To make the result a m×1 matrix instead, do A[,2,drop=FALSE] A[7,] Note: that gives the result as a vector. To make the result a 1×n matrix instead, do A[7,,drop=FALSE] c(A) A[2:4,6:10] A[c(7,7,6),c(2,1)]

[r,c] = ind2sub(size(A), ind)

r = ((ind-1) %% m) + 1 c = floor((ind-1) / m) + 1

ind = sub2ind(size(A), r, c)

ind = (c-1)*m + r

inds = sub2ind(size(A),r,c); A(inds) = 12;

inds = cbind(r,c) A[inds] = 12

v = v(1:10) A = reshape(A,m,n)

v = v[1:10], or length(v) = 10 also works dim(A) = c(m,n)

L = tril(A) U = triu(A) hilb(n)

L = A; L[upper.tri(L)]=0 U = A; U[lower.tri(U)]=0 Hilbert(n), but this is part of the Matrix package which you’ll need to install (see item 295 for how to install/load packages). array(1:24, c(3,4,2)) (Note that a matrix is 2-D, i.e. rows and columns, while an array is more generally N -D)

39

reshape(1:24, 3, 4, 2) reshape(1:24, [3 4 2])

or

D. Hiebeler, Matlab / R Reference

6

2.1
No. 40

Cell arrays and lists
Matlab v = cell(1,n) In general, cell(m,n) makes an m × n cell array. Then you can do e.g.: v{1} = 12 v{2} = ’hi there’ v{3} = rand(3) R v = vector(’list’,n) can do e.g.: Then you

Description Build a vector v of length n, capable of containing different data types in different elements (called a cell array in Matlab, and a list in R)

v[[1]] = 12 v[[2]] = ’hi there’ v[[3]] = matrix(runif(9),3)

41

Extract the ith element of a cell/list vector v

w = v{i} If you use regular indexing, i.e. w = v(i), then w will be a 1 × 1 cell matrix containing the contents of the ith element of v. (Matlab does not have names associated with elements of cell arrays.)

w = v[[i]] If you use regular indexing, i.e. w = v[i], then w will be a list of length 1 containing the contents of the ith element of v. names(v)[3] = ’myrandmatrix’ Use names(v) to see all names, and names(v)=NULL to clear all names.

42

Set the name of the ith element in a list.

2.2
No. 43

Structs and data frames

Description Matlab R Create a matrix-like object avals=2*ones(1,6); v=c(1,5,3,2,3,7); d=data.frame( with different named columns yvals=6:-1:1; v=[1 5 3 2 3 7]; cbind(a=2, yy=6:1), v) (a struct in Matlab, or a d=struct(’a’,avals, data frame in R) ’yy’, yyvals, ’fac’, v); Note that I (surprisingly) don’t use R for statistics, and therefore have very little experience with data frames (and also very little with Matlab structs). I will try to add more to this section later on.

D. Hiebeler, Matlab / R Reference

7

3
3.1
No. 44 45 46 47

Computations
Basic computations
Matlab a+b, a-b, a*b, a/b sqrt(a) a^b abs(a) R a+b, a-b, a*b, a/b sqrt(a) a^b abs(a) Description a + b, a − b, ab, a/b √ a ab |a| (note: for complex arguments, this computes the modulus) ea ln(a) log2 (a), log10 (a) sin(a), cos(a), tan(a) sin−1 (a), cos−1 (a), tan−1 (a) sinh(a), cosh(a), tanh(a) sinh−1 (a), cosh−1 (a), −1 tanh (a) n MOD k (modulo arithmetic) Round to nearest integer

48 49 50 51 52 53 54 55 56

exp(a) log(a) log2(a), log10(a) sin(a), cos(a), tan(a) asin(a), acos(a), atan(a) sinh(a), cosh(a), tanh(a) asinh(a), acosh(a), atanh(a) mod(n,k) round(x)

exp(a) log(a) log2(a), log10(a) sin(a), cos(a), tan(a) asin(a), acos(a), atan(a) sinh(a), cosh(a), tanh(a) asinh(a), acosh(a), atanh(a) n %% k round(x) (Note: R uses IEC 60559 standard, rounding 5 to the even digit — so e.g. round(0.5) gives 0, not 1.) floor(x) ceiling(x) sign(x) (Does not work with complex values) 2*pnorm(x*sqrt(2))-1 2*pnorm(x*sqrt(2),lower=FALSE)

57 58 59 60 61

Round down to next lowest integer Round up to next largest integer Sign of x (+1, 0, or -1) Error function erf(x) √ 2 x (2/ π) 0 e−t dt =

floor(x) ceil(x) sign(x) (Note: for complex values, this computes x/abs(x).) erf(x) erfc(x)

Complementary error function cerf(x) = √ 2 ∞ (2/ π) x e−t dt = 1-erf(x)

62 63

Inverse error function erfinv(x) qnorm((1+x)/2)/sqrt(2) Inverse complementary error erfcinv(x) qnorm(x/2,lower=FALSE)/sqrt(2) function Note: the various functions above (logarithm, exponential, trig, abs, and rounding functions) all work with vectors and matrices, applying the function to each element, as well as with scalars.

3.2
No. 64 65 66 67 68 69

Complex numbers
Matlab 1+2i abs(z) angle(z) conj(z) real(z) imag(z) R 1+2i abs(z) or Mod(z) Arg(z) Conj(z) Re(z) Im(z)

Description Enter a complex number Modulus (magnitude) Argument (angle) Complex conjugate Real part of z Imaginary part of z

D. Hiebeler, Matlab / R Reference

8

3.3
No. 70 71 72

Matrix/vector computations
Matlab A * B A .* B A’ (This is actually the complex conjugate (i.e. Hermitian) transpose; use A.’ for the non-conjugate transpose if you like; they are equivalent for real matrices.) A\b Warning: if there is no solution, Matlab gives you a least-squares “best fit.” If there are many solutions, Matlab just gives you one of them. rref(A) inv(A) A/B A ./ B A\B A^2 A^k R A %*% B A * B t(A) for transpose, or Conj(t(A)) for conjugate (Hermitian) transpose

Description Matrix multiplication AB Element-by-element multiplication of A and B Transpose of a matrix, AT

73

Solve Ax = b

solve(A,b) Warning: this only works with square invertible matrices.

74 75 76 77 78 79 80

Reduced echelon form of A Compute inverse of A Compute AB −1 Element-by-element division of A and B Compute A−1 B Square the matrix A Raise matrix A to the k th power Raise each element of A to the k th power Rank of matrix A Set w to be a vector of eigenvalues of A, and V a matrix containing the corresponding eigenvectors Permuted LU factorization of a matrix

R does not have a function to do this solve(A) A %*% solve(B) A / B solve(A,B) A %*% A (No easy way to do this in R other than repeated multiplication A %*% A %*% A...) A^k qr(A)$rank tmp=eigen(A); w=tmp$values; V=tmp$vectors

81 82 83

A.^k rank(A) [V,D]=eig(A) and then w=diag(D) since Matlab returns the eigenvalues on the diagonal of D [L,U,P]=lu(A) then the matrices satisfy P A = LU . Note that this works even with non-square matrices

84

tmp=expand(lu(Matrix(A))); L=tmp$L; U=tmp$U; P=tmp$P then the matrices satisfy A = P LU , i.e. P −1 A = LU . Note that the lu and expand functions are part of the Matrix package (see item 295 for how to install/load packages). Also note that this doesn’t seem to work correctly with non-square matrices. L, U, and P will be of class Matrix rather than class matrix; to make them the latter, instead do L=as.matrix(tmp$L), U=as.matrix(tmp$U), and P=as.matrix(tmp$P) above.

D. Hiebeler, Matlab / R Reference No. 85 Description Singular-value decomposition: given m × n matrix A with rank r, find m × r matrix P with orthonormal columns, diagonal r × r matrix S, and r × n matrix QT with orthonormal rows so that P SQT = A Schur decomposition of square matrix, A = QT QH = QT Q−1 where Q is unitary (i.e. QH Q = I) and T is upper triangular; QH = QT is the Hermitian (conjugate) transpose Cholesky factorization of a square, symmetric, positive definite matrix A = RT R, where R is upper-triangular QR factorization of matrix A, where Q is orthogonal (satisfying QQT = I) and R is upper-triangular Vector norms Matlab [P,S,Q]=svd(A,’econ’)

9 R tmp=svd(A); U=tmp$u; V=tmp$v; S=diag(tmp$d)

86

[Q,T]=schur(A)

87

R = chol(A)

tmp=Schur(Matrix(A)); T=tmp@T; Q=tmp@Q Note that Schur is part of the Matrix package (see item 295 for how to install/load packages). T and Q will be of class Matrix rather than class matrix; to make them the latter, instead do T=as.matrix(tmp@T) and Q=as.matrix(tmp@Q) above. R = chol(A) Note that chol is part of the Matrix package (see item 295 for how to install/load packages). z=qr(A); Q=qr.Q(z); R=qr.R(z); E=diag(n)[,z$pivot] (where n is the number of columns in A) gives permuted QR factorization satisfying AE = QR R does not have a norm function for vectors; only one for matrices. But the following will work: norm(matrix(v),’1’) for 1-norm v 1 , norm(matrix(v),’i’) for infinity-norm v ∞, and sum(abs(v)^p)^(1/p) for p-norm 1/p v p = ( |vi |p ) norm(A,’1’) for 1-norm A 1, max(svd(A)$d) for 2-norm A 2 , norm(A,’i’) for infinity-norm A ∞ , and norm(A,’f’) for Frobenius norm 1/2 T i (A A)ii 1/rcond(A,’1’)

88

[Q,R]=qr(A) satisfying QR = A, or [Q,R,E]=qr(A) to do permuted QR factorization satisfying AE = QR

89

norm(v,1) for 1-norm v 1, norm(v,2) for Euclidean norm v 2 , norm(v,inf) for infinity-norm v ∞ , and norm(v,p) for p-norm 1/p v p = ( |vi |p )

90

Matrix norms

91

Condition number cond(A) = A 1 A−1 1 of A, using 1norm Condition number cond(A) = A 2 A−1 2 of A, using 2norm Condition number cond(A) = A ∞ A−1 ∞ of A, using infinity-norm

92

norm(A,1) for 1-norm A 1, norm(A) for 2-norm A 2, norm(A,inf) for infinity-norm A ∞ , and norm(A,’fro’) for 1/2 T Frobenius norm i (A A)ii cond(A,1) (Note: Matlab also has a function rcond(A) which computes reciprocal condition estimator using the 1-norm) cond(A,2)

93

cond(A,inf)

kappa(A, exact=TRUE) (leave out the “exact=TRUE” for an estimate) 1/rcond(A,’I’)

D. Hiebeler, Matlab / R Reference No. 94 95 96 97 Description Compute mean of all elements in vector or matrix Compute means of columns of a matrix Compute means of rows of a matrix Compute standard deviation of all elements in vector or matrix Compute standard deviations of columns of a matrix Compute standard deviations of rows of a matrix Compute variance of all elements in vector or matrix Compute variance of columns of a matrix Compute variance of rows of a matrix Compute covariance for two vectors of observations Compute covariance matrix, giving covariances between columns of matrix A Given matrices A and B, build covariance matrix C where cij is the covariance between column i of A and column j of B Compute Pearson’s linear correlation coefficient between elements of vectors v and w Compute Kendall’s tau correlation statistic for vectors v and w Compute Spearman’s rho correlation statistic for vectors v and w Compute pairwise Pearson’s correlation coefficient between columns of matrix A Compute matrix C of pairwise Pearson’s correlation coefficients between each pair of columns of matrices A and B, i.e. so cij is the correlation between column i of A and column j of B Matlab mean(v) for vectors, mean(A(:)) for matrices mean(A) mean(A,2) std(v) for vectors, std(A(:)) for matrices. This normalizes by n − 1. Use std(v,1) to normalize by n. std(A). This normalizes by n − 1. Use std(A,1) to normalize by n std(A,0,2) to normalize by n − 1, std(A,1,2) to normalize by n var(v) for vectors, var(A(:)) for matrices. This normalizes by n − 1. Use var(v,1) to normalize by n. var(A). This normalizes by n − 1. Use var(A,1) to normalize by n var(A,0,2) to normalize by n − 1, var(A,1,2) to normalize by n cov(v,w) computes the 2 × 2 covariance matrix; the off-diagonal elements give the desired covariance cov(A) R mean(v) or mean(A) colMeans(A) rowMeans(A)

10

sd(v) for vectors, sd(c(A)) for matrices. This normalizes by n − 1. sd(A). This normalizes by n − 1. apply(A,1,sd). This normalizes by n − 1. var(v) for vectors, var(c(A)) for matrices. This normalizes by n − 1. apply(A,2,var). This normalizes by n − 1. apply(A,1,var). This normalizes by n − 1. cov(v,w)

98 99 100

101 102 103

104

var(A) or cov(A)

105

I don’t know of a direct way to do this in Matlab. But one way is [Y,X]=meshgrid(std(B),std(A)); X.*Y.*corr(A,B) corr(v,w) Note: v and w must be column vectors. To make it work regardless of whether they are row or column vectors, do corr(v(:),w(:)) corr(v,w,’type’,’kendall’)

cov(A,B)

106

cor(v,w)

107

cor(v,w,method=’kendall’)

108

corr(v,w,’type’,’spearman’)

cor(v,w,method=’spearman’)

109

corr(A) The ’type’ argument may also be used as in the previous two items corr(A,B) The ’type’ argument may also be used as just above

cor(A) The method argument may also be used as in the previous two items cor(A,B) The method argument may also be used as just above

110

D. Hiebeler, Matlab / R Reference No. 111 112 113 114 Description Compute sum of all elements in vector or matrix Compute sums of columns of matrix Compute sums of rows of matrix Compute matrix exponential ∞ eA = k=0 Ak /k! Compute cumulative sum of values in vector Compute cumulative sums of columns of matrix Compute cumulative sums of rows of matrix Compute cumulative sum of all elements of matrix (column-by-column) Cumulative product of elements in vector v Cumulative minimum or maximum of elements in vector v Compute differences between consecutive elements of vector v. Result is a vector w 1 element shorter than v, where element i of w is element i + 1 of v minus element i of v Make a vector y the same size as vector x, which equals 4 everywhere that x is greater than 5, and equals 3 everywhere else (done via a vectorized computation). Compute minimum of values in vector v Matlab sum(v) for vectors, sum(A(:)) for matrices sum(A) sum(A,2) expm(A) R sum(v) or sum(A) colSums(A) rowSums(A)

11

115 116 117 118

cumsum(v) cumsum(A) cumsum(A,2) cumsum(A(:))

expm(Matrix(A)), but this is part of the Matrix package which you’ll need to install (see item 295 for how to install/load packages). cumsum(v) apply(A,2,cumsum) t(apply(A,1,cumsum)) cumsum(A)

119 120

cumprod(v) (Can also be used in the various ways cumsum can) I don’t know of an easy way to do this in Matlab diff(v)

cumprod(v) (Can also be used in the various ways cumsum can) cummin(v) or cummax(v)

121

diff(v)

122

z = [3 4]; y = z((x > 5)+1)

y = ifelse(x > 5, 4, 3)

123

min(v)

min(v)

D. Hiebeler, Matlab / R Reference No. 124 125 126 127 Description Compute minimum of all values in matrix A Compute minimum value of each column of matrix A Compute minimum value of each row of matrix A Given matrices A and B, compute a matrix where each element is the minimum of the corresponding elements of A and B Given matrix A and scalar c, compute a matrix where each element is the minimum of c and the corresponding element of A Find minimum among all values in matrices A and B Find index of the first time min(v) appears in v, and store that index in ind Notes: Matlab min(A(:)) min(A) (returns a row vector) min(A, [ ], 2) (returns a column vector) min(A,B) R min(A)

12

apply(A,2,min) (returns a vector) apply(A,1,min) (returns a vector) pmin(A,B)

128

min(A,c)

pmin(A,c)

129 130

min([A(:)

; B(:)])

min(A,B) ind = which.min(v)

[y,ind] = min(v)

• Matlab and R both have a max function (and R has pmax and which.max as well) which behaves in the same ways as min but to compute maxima rather than minima. • Functions like exp, sin, sqrt etc. will operate on arrays in both Matlab and R, doing the computations for each element of the matrix. No. 131 132 133 134 135 136 137 138 Description Number of rows in A Number of columns in A Dimensions of A, listed in a vector Number of elements in vector v Total number of elements in matrix A Max. dimension of A Sort values in vector v Sort values in v, putting sorted values in s, and indices in idx, in the sense that s[k] = x[idx[k]] To count how many values in the vector x are between 4 and 7 (inclusive on the upper end) Given vector v, return list of indices of elements of v which are greater than 5 Matlab size(A,1) size(A,2) size(A) length(v) numel(A) length(A) sort(v) [s,idx]=sort(v) R nrow(A) ncol(A) dim(A) length(v) length(A) max(dim(A)) sort(v) tmp=sort(v,index.return=TRUE); s=tmp$x; idx=tmp$ix

139

sum((x > 4) & (x <= 7))

sum((x > 4) & (x <= 7))

140

find(v > 5)

which(v > 5)

D. Hiebeler, Matlab / R Reference No. 141 Description Given matrix A, return list of indices of elements of A which are greater than 5, using single-indexing Given matrix A, generate vectors r and c giving rows and columns of elements of A which are greater than 5 Given vector x (of presumably discrete values), build a vector v listing unique values in x, and corresponding vector c indicating how many times those values appear in x Given vector x (of presumably continuous values), divide the range of values into k equally-sized bins, and build a vector m containing the midpoints of the bins and a corresponding vector c containing the counts of values in the bins Convolution / polynomial multiplication (given vectors x and y containing polynomial coefficients, their convolution is a vector containing coefficients of the product of the two polynomials) Matlab find(A > 5) R which(A > 5)

13

142

[r,c] = find(A > 5)

w = which(A > 5, arr.ind=TRUE); r=w[,1]; c=w[,2]

143

v = unique(x); c = hist(x,v);

w=table(x); c=as.numeric(w); v=as.numeric(names(w))

144

[c,m] = hist(x,k)

w=hist(x,seq(min(x),max(x), length.out=k+1), plot=FALSE); m=w$mids; c=w$counts

145

conv(x,y)

convolve(x,rev(y),type=’open’) Note: the accuracy of this is not as good as Matlab; e.g. doing v=c(1,-1); for (i in 2:20) v=convolve(v,c(-i,1), type=’open’) to generate the 20th -degree Wilkinson polynomial 20 W (x) = i=1 (x−i) gives a coefficient of ≈ −780.19 for x19 , rather than the correct value -210.

3.4
No. 146

Root-finding
Matlab roots(v) R polyroot(rev(v)) (This function really wants the vector to have the constant coefficient first in v; rev reverses their order to achieve this.) Define function f(x), then do uniroot(f, c(a,b)) to find a root between a and b, assuming the sign of f (x) differs at x = a and x = b. Default forward error tolerance (i.e. error in x) is fourth root of machine epsilon, (ǫmach )0.25 . To specify e.g. a tolerance of 2−52 , do uniroot(f, c(a,b), tol=2^-52).

147

Description Find roots of polynomial whose coefficients are stored in vector v (coefficients in v are highest-order first) Find zero (root) of a function f (x) of one variable

Define function f(x), then do fzero(f,x0) to search for a root near x0, or fzero(f,[a b]) to find a root between a and b, assuming the sign of f (x) differs at x = a and x = b. Default forward error tolerance (i.e. error in x) is machine epsilon ǫmach .

D. Hiebeler, Matlab / R Reference

14

3.5
No. 148

Function optimization/minimization
Matlab Define function f(x), then do m = fminbnd(f, a, b) Define function f(x,p1,p2), then use an “anonymous function”: % first define values for p1 % and p2, and then do: m=fminbnd(@(x) f(x,p1,p2),a,b) First write function f(v) which accepts a vector argument v containing values of x, y, and z, and returns the scalar value f (x, y, z), then do: fminsearch(@f,[1 2.2 3.4]) R Define function f(x), then do m = optimize(f,c(a,b))$minimum Define function f(x,p1,p2), then: # first define values for p1 # and p2, and then do: m = optimize(f, c(a,b), p1=p1, p2=p2)$minimum First write function f(v) which accepts a vector argument v containing values of x, y, and z, and returns the scalar value f (x, y, z), then do: optim(c(1,2.2,3.4),f)$par First write function f(v,p1,p2) which accepts a vector argument v containing values of x, y, and z, along with the extra parameters, and returns the scalar value f (x, y, z, p1 , p2 ), then do: optim(c(1,2.2,3.4), f, p1=p1, p2=p2)$par

149

150

Description Find value m which minimizes a function f (x) of one variable within the interval from a to b Find value m which minimizes a function f (x, p1 , p2 ) with given extra parameters (but minimization is only occuring over the first argument), in the interval from a to b. Find values of x, y, z which minimize function f (x, y, z), using a starting guess of x = 1, y = 2.2, and z = 3.4.

151

Find values of x, y, z which minimize function f (x, y, z, p1 , p2 ), using a starting guess of x = 1, y = 2.2, and z = 3.4, where the function takes some extra parameters (useful e.g. for doing things like nonlinear least-squares optimization where you pass in some data vectors as extra parameters).

First write function f(v,p1,p2) which accepts a vector argument v containing values of x, y, and z, along with the extra parameters, and returns the scalar value f (x, y, z, p1 , p2 ), then do: fminsearch(@f,[1 2.2 3.4], ... [ ], p1, p2) Or use an anonymous function: fminsearch(@(x) f(x,p1,p2), ... [1 2.2 3.4])

3.6
No. 152

Numerical integration / quadrature
Matlab quad(f,a,b) uses son’s quadrature, absolute tolerance specify absolute quad(f,a,b,tol) adaptive Simpwith a default of 10−6 . To tolerance, use R integrate(f,a,b) uses adaptive quadrature with default absolute and relative error tolerances being the fourth root of machine epsilon, (ǫmach )0.25 ≈ 1.22 × 10−4 . Tolerances can be specified by using integrate(f,a,b, rel.tol=tol1, abs.tol=tol2). Note that the function f must be written to work even when given a vector of x values as its argument.

Description Numerically integrate function f (x) over interval from a to b

D. Hiebeler, Matlab / R Reference

15

3.7
No. 153

Curve fitting
Matlab p = polyfit(x,y,1) The return vector p has the coefficients in descending order, i.e. p(1) is c1 , and p(2) is c0 . R p = coef(lm(y ~ x)) The return vector p has the coefficients in ascending order, i.e. p[1] is c0 , and p[2] is c1 . p = coef(lm(y ~ x + I(x^2))) The return vector p has the coefficients in ascending order, i.e. p[1] is c0 , p[2] is c1 , and p[3] is c2 . There isn’t a simple function built into the standard R distribution to do this, but see the polyreg function in the mda package (see item 295 for how to install/load packages).

Description Fit the line y = c1 x + c0 to data in vectors x and y.

154

Fit the quadratic polynomial y = c2 x2 + c1 x + c0 to data in vectors x and y.

p = polyfit(x,y,2) The return vector p has the coefficients in descending order, i.e. p(1) is c2 , p(2) is c1 , and p(3) is c0 .

155

Fit nth degree polynomial y = cn xn + cn−1 xn−1 + . . . + c1 x + c0 to data in vectors x and y.

p = polyfit(x,y,n) The return vector p has the coefficients in descending order, p(1) is cn , p(2) is cn−1 , etc. (I don’t know a simple way do this in Matlab, other than to write a function which computes the sum of squared residuals and use fminsearch on that function. There is likely an easy way to do it in the Statistics Toolbox.) pp=csape(x,y,’variational’); yy=ppval(pp,xx) but note that csape is in Matlab’s Spline Toolbox

156

Fit the quadratic polynomial with zero intercept, y = c2 x2 + c1 x to data in vectors x and y.

p=coef(lm(y ~ -1 + x + I(x^2))) The return vector p has the coefficients in ascending order, i.e. p[1] is c1 , and p[2] is c2 . tmp=spline(x,y,method=’natural’, xout=xx); yy=tmp$y

157

158

Fit natural cubic spline (S ′′ (x) = 0 at both endpoints) to points (xi , yi ) whose coordinates are in vectors x and y; evaluate at points whose x coordinates are in vector xx, storing corresponding y’s in yy Fit cubic spline using Forsythe, Malcolm and Moler method (third derivatives at endpoints match third derivatives of exact cubics through the four points at each end) to points (xi , yi ) whose coordinates are in vectors x and y; evaluate at points whose x coordinates are in vector xx, storing corresponding y’s in yy

I’m not aware of a function to do this in Matlab

tmp=spline(x,y,xout=xx); yy=tmp$y

D. Hiebeler, Matlab / R Reference No. 159 Description Fit cubic spline such that first derivatives at endpoints match first derivatives of exact cubics through the four points at each end) to points (xi , yi ) whose coordinates are in vectors x and y; evaluate at points whose x coordinates are in vector xx, storing corresponding y’s in yy Fit cubic spline with periodic boundaries, i.e. so that first and second derivatives match at the left and right ends (the first and last y values of the provided data should also agree), to points (xi , yi ) whose coordinates are in vectors x and y; evaluate at points whose x coordinates are in vector xx, storing corresponding y’s in yy Fit cubic spline with “nota-knot” conditions (the first two piecewise cubics coincide, as do the last two), to points (xi , yi ) whose coordinates are in vectors x and y; evaluate at points whose x coordinates are in vector xx, storing corresponding y’s in yy Matlab pp=csape(x,y); yy=ppval(pp,xx) but csape is in Matlab’s Spline Toolbox

16 R I’m not aware of a function to do this in R

160

pp=csape(x,y,’periodic’); yy=ppval(pp,xx) but csape is in Matlab’s Spline Toolbox

tmp=spline(x,y,method= ’periodic’, xout=xx); yy=tmp$y

161

yy=spline(x,y,xx)

I’m not aware of a function to do this in R

4
No. 162

Conditionals, control structure, loops
Description “for” loops over values in a vector v (the vector v is often constructed via a:b) Matlab for i=v command1 command2 end R If only one command inside the loop: for (i in v) command or for (i in v) command If multiple commands inside the loop: for (i in v) { command1 command2 }

D. Hiebeler, Matlab / R Reference No. 163 Description “if” statements with no else clause Matlab if cond command1 command2 end

17 R If only one command inside the clause: if (cond) command or if (cond) command If multiple commands: if (cond) { command1 command2 }

164

“if/else” statement if cond command1 command2 else command3 command4 end Note: Matlab also has an “elseif” statement, e.g.: if cond1 command1 elseif cond2 command2 elseif cond3 command3 else command4 end

If one command in clauses: if (cond) command1 else command2 or if (cond) cmd1 else cmd2 If multiple commands: if (cond) { command1 command2 } else { command3 command4 } Warning: the “else” must be on the same line as command1 or the “}” (when typed interactively at the command prompt), otherwise R thinks the “if” statement was finished and gives an error. R does not have an “elseif” statement.

Logical comparisons which can be used on scalars in “if” statements, or which operate element-byelement on vectors/matrices: Matlab x<a x>a x <= a x >= a x == a x ~= a R x x x x x x <a >a <= a >= a == a != a Description True if x is less than a True if x is greater than a True if x is less than or equal to a True if x is greater than or equal to a True if x is equal to a True if x is not equal to a

D. Hiebeler, Matlab / R Reference Scalar logical operators: Description a AND b a OR b a XOR b NOT a Matlab a && b a || b xor(a,b) ~a R a && b a || b xor(a,b) !a

18

The && and || operators are short-circuiting, i.e. && stops as soon as any of its terms are FALSE, and || stops as soon as any of its terms are TRUE.

Matrix logical operators (they operate element-by-element): Description a AND b a OR b a XOR b NOT a No. 165 Description To test whether a scalar value x is between 4 and 7 (inclusive on the upper end) To count how many values in the vector x are between 4 and 7 (inclusive on the upper end) Test whether all values in a logical/boolean vector are TRUE Test whether any values in a logical/boolean vector are TRUE Matlab a & b a | b xor(a,b) ~a R a & b a | b xor(a,b) !a R if ((x > 4) && (x <= 7))

Matlab if ((x > 4) && (x <= 7))

166

sum((x > 4) & (x <= 7))

sum((x > 4) & (x <= 7))

167

all(v)

all(v)

168

any(v)

any(v)

No. 169

Description “while” statements to do iteration (useful when you don’t know ahead of time how many iterations you’ll need). E.g. to add uniform random numbers between 0 and 1 (and their squares) until their sum is greater than 20:

Matlab mysum = 0; mysumsqr = 0; while (mysum < 20) r = rand; mysum = mysum + r; mysumsqr = mysumsqr + r^2; end

R mysum = 0 mysumsqr = 0 while (mysum < 20) { r = runif(1) mysum = mysum + r mysumsqr = mysumsqr + r^2 } (As with “if” statements and “for” loops, the curly brackets are not necessary if there’s only one statement inside the “while” loop.)

D. Hiebeler, Matlab / R Reference No. 170 Description “Switch” statements for integers Matlab switch (x) case 10 disp(’ten’) case {12,13} disp(’dozen (bakers?)’) otherwise disp(’unrecognized’) end

19 R R doesn’t have a switch statement capable of doing this. It has a function which is fairly limited for integers, but can which do string matching. See ?switch for more. But a basic example of what it can do for integers is below, showing that you can use it to return different expressions based on whether a value is 1, 2, . . .. mystr = switch(x, ’one’, ’two’, ’three’) print(mystr) Note that switch returns NULL if x is larger than 3 in the above case. Also, continuous values of x will be truncated to integers.

5
No. 171

Functions, ODEs
Description Implement add(x,y) a function Matlab Put the following in add.m: function retval=add(x,y) retval = x+y; Then you can do e.g. add(2,3) R Enter the following, or put it in a file and source that file: add = function(x,y) { return(x+y) } Then you can do e.g. add(2,3). Note, the curly brackets aren’t needed if your function only has one line. Write function as follows: f = function(x,y,z) { a = x*y+z; b=2*sin(x-z) return(list(a,b)) } Then call the function by doing: tmp=f(2,8,12); u=tmp[[1]]; v=tmp[[2]]. The above is most general, and will work even when u and v are different types of data. If they are both scalars, the function could simply return them packed in a vector, i.e. return(c(a,b)). If they are vectors of the same size, the function could return them packed together into the columns of a matrix, i.e. return(cbind(a,b)).

172

Implement a function f(x,y,z) which returns multiple values, and store those return values in variables u and v

Write function as follows: function [a,b] = f(x,y,z) a = x*y+z; b=2*sin(x-z); Then call the function by doing: [u,v] = f(2,8,12)

D. Hiebeler, Matlab / R Reference No. 173 Description Numerically solve ODE dx/dt = 5x from t = 3 to t = 12 with initial condition x(3) = 7 Matlab First implement function function retval=f(t,x) retval = 5*x; Then do ode45(@f,[3,12],7) to plot solution, or [t,x]=ode45(@f,[3,12],7) to get back vector t containing time values and vector x containing corresponding function values. If you want function values at specific times, e.g. 3, 3.1, 3.2, . . . , 11.9, 12, you can do [t,x]=ode45(@f,3:0.1:12,7). Note: in older versions of Matlab, use ’f’ instead of @f. First implement function function retval=myfunc(t,x) w = x(1); z = x(2); retval = zeros(2,1); retval(1) = 5*w; retval(2) = 3*w + 7*z; Then do ode45(@myfunc,[3,12],[7; 8.2]) to plot solution, or [t,x]=ode45(@myfunc,[3,12],[7; 8.2]) to get back vector t containing time values and matrix x, whose first column containing corresponding w(t) values and second column contains z(t) values. If you want function values at specific times, e.g. 3, 3.1, 3.2, . . . , 11.9, 12, you can do [t,x]=ode45(@myfunc,3:0.1:12,[7; 8.2]). Note: in older versions of Matlab, use ’f’ instead of @f. First implement function function retval=func2(t,x,r,K) retval = r*x*(1-x/K) Then do ode45(@func2,[0 20], 2.5, [ ], 1.3, 50). The empty matrix is necessary between the initial condition and the beginning of your extra parameters. R First implement function

20

f = function(t,x,parms) { return(list(5*x)) } Then do y=lsoda(7, seq(3,12, 0.1), f,NA) to obtain solution values at times 3, 3.1, 3.2, . . . , 11.9, 12. The first column of y, namely y[,1] contains the time values; the second column y[,2] contains the corresponding function values. Note: lsoda is part of the deSolve package (see item 295 for how to install/load packages). First implement function myfunc = function(t,x,parms) { w = x[1]; z = x[2]; return(list(c(5*w, 3*w+7*z))) } Then do y=lsoda(c(7,8.2), seq(3,12, 0.1), myfunc,NA) to obtain solution values at times 3, 3.1, 3.2, . . . , 11.9, 12. The first column of y, namely y[,1] contains the time values; the second column y[,2] contains the corresponding values of w(t); and the third column contains z(t). Note: lsoda is part of the deSolve package (see item 295 for how to install/load packages).

174

Numerically solve system of ODEs dw/dt = 5w, dz/dt = 3w + 7z from t = 3 to t = 12 with initial conditions w(3) = 7, z(3) = 8.2

175

Pass parameters such as r = 1.3 and K = 50 to an ODE function from the command line, solving dx/dt = rx(1 − x/K) from t = 0 to t = 20 with initial condition x(0) = 2.5.

First implement function func2=function(t,x,parms) { r=parms[1]; K=parms[2] return(list(r*x*(1-x/K))) } Then do y=lsoda(2.5,seq(0,20,0.1) func2,c(1.3,50)) Note: lsoda is part of the deSolve package (see item 295 for how to install/load packages).

D. Hiebeler, Matlab / R Reference

21

6
No. 176

Probability and random values
Description Generate a continuous uniform random value between 0 and 1 Generate vector of n uniform random vals between 0 and 1 Generate m×n matrix of uniform random values between 0 and 1 Generate m×n matrix of continuous uniform random values between a and b Generate a random integer between 1 and k Matlab rand R runif(1)

177 178

rand(n,1) or rand(1,n) rand(m,n)

runif(n) matrix(runif(m*n),m,n) matrix(runif(m*n),m) matrix(runif(m*n,a,b),m) or just

179

180

a+rand(m,n)*(b-a) or if you have the Statistics toolbox then unifrnd(a,b,m,n) floor(k*rand) + 1

181

182

Generate m×n matrix of discrete uniform random integers between 1 and k Generate m × n matrix where each entry is 1 with probability p, otherwise is 0

floor(k*rand(m,n))+1 or if you have the Statistics toolbox then unidrnd(k,m,n) (rand(m,n)<p)*1 Note: multiplying by 1 turns the logical (true/false) result back into numeric values. You could also do double(rand(m,n)<p)

floor(k*runif(1)) + 1 Note: sample(k)[1] would also work, but I believe in general will be less efficient, because that actually generates many random numbers and then just uses one of them. floor(k*matrix(runif(m*n),m))+1

183

184

Generate m × n matrix where each entry is a with probability p, otherwise is b Generate a random integer between a and b inclusive Flip a coin which comes up heads with probability p, and perform some action if it does come up heads Generate a random permutation of the integers 1, 2, . . . , n Generate a random selection of k unique integers between 1 and n (i.e. sampling without replacement)

b + (a-b)*(rand(m,n)<p)

(matrix(runif(m,n),m)<p)*1 (Note: multiplying by 1 turns the logical (true/false) result back into numeric values; using as.numeric() to do it would lose the shape of the matrix.) b + (a-b)*(matrix( runif(m,n),m)<p) floor((b-a+1)*runif(1))+a

floor((b-a+1)*rand)+a or if you have the Statistics toolbox then unidrnd(b-a+1)+a-1 if (rand < p) ...some commands... end randperm(n) [s,idx]=sort(rand(n,1)); ri=idx(1:k) or another way is ri=randperm(n); ri=ri(1:k). Or if you have the Statistics Toolbox, then randsample(n,k)

185

if (runif(1) < p) { ...some commands... } sample(n) ri=sample(n,k)

186 187

D. Hiebeler, Matlab / R Reference No. 188 Description Choose k values (with replacement) from the vector v, storing result in w Choose k values (without replacement) from the vector v, storing result in w Set the random-number generator back to a known state (useful to do at the beginning of a stochastic simulation when debugging, so you’ll get the same sequence of random numbers each time) Matlab L=length(v); w=v(floor(L*rand(k,1))+1) Or, if you have the Statistics Toolbox, w=randsample(v,k,replace=true) L=length(v); ri=randperm(L); ri=ri(1:k); w=v(ri) Or, if you have the Statistics Toolbox, w=randsample(v,k,replace=false) rand(’state’, 12)

22 R w=sample(v,k,replace=TRUE)

189

w=sample(v,k,replace=FALSE)

190

set.seed(12)

No. 191

192

193

194

195

196

197

Note that the “*rnd,” “*pdf,” and “*cdf” functions described below are all part of the Matlab Statistics Toolbox, and not part of the core Matlab distribution. Description Matlab R Generate a random value binornd(n,p) rbinom(1,n,p) from the Binomial(n, p) distribution Generate a random value poissrnd(lambda) rpois(1,lambda) from the Poisson distribution with parameter λ Generate a random value exprnd(mu) or -mu*log(rand) will rexp(1, 1/mu) from the Exponential distri- work even without the Statistics bution with mean µ Toolbox. Generate a random value unidrnd(k) or floor(rand*k)+1 sample(k,1) from the discrete uniform dis- will work even without the Statistics tribution on integers 1 . . . k Toolbox. Generate n iid random values unidrnd(k,n,1) or sample(k,n,replace=TRUE) from the discrete uniform dis- floor(rand(n,1)*k)+1 will work tribution on integers 1 . . . k even without the Statistics Toolbox. Generate a random value unifrnd(a,b) or (b-a)*rand + a runif(1,a,b) from the continuous uniform will work even without the Statistics distribution on the interval Toolbox. (a, b) Generate a random value normrnd(mu,sigma) or rnorm(1,mu,sigma) from the normal distribution mu + sigma*randn will work with mean mu and standard even without the Statistics Toolbox. deviation σ Notes: • The Matlab “*rnd” functions above can all take additional r,c arguments to build an r × c matrix of iid random values. E.g. poissrnd(3.5,4,7) for a 4 × 7 matrix of iid values from the Poisson distribution with mean λ = 3.5. The unidrnd(n,k,1) command above is an example of this, to generate a k × 1 column vector. • The first parameter of the R “r*” functions above specifies how many values are desired. E.g. to generate 28 iid random values from a Poisson distribution with mean 3.5, use rpois(28,3.5). To get a 4 × 7 matrix of such values, use matrix(rpois(28,3.5),4).

D. Hiebeler, Matlab / R Reference No. 198 Description Compute probability that a random variable from the Binomial(n, p) distribution has value x (i.e. the density, or pdf). Compute probability that a random variable from the Poisson(λ) distribution has value x. Matlab binopdf(x,n,p) or nchoosek(n,x)*p^x*(1-p)^(n-x) will work even without the Statistics Toolbox, as long as n and x are non-negative integers and 0 ≤ p ≤ 1. poisspdf(x,lambda) or exp(-lambda)*lambda^x / factorial(x) will work even without the Statistics Toolbox, as long as x is a non-negative integer and lambda ≥ 0. exppdf(x,mu) or (x>=0)*exp(-x/mu)/mu will work even without the Statistics Toolbox, as long as mu is positive. normpdf(x,mu,sigma) or exp(-(x-mu)^2/(2*sigma^2))/ (sqrt(2*pi)*sigma) will work even without the Statistics Toolbox. R dbinom(x,n,p)

23

199

dpois(x,lambda)

200

Compute probability density dexp(x,1/mu) function at x for a random variable from the exponential distribution with mean µ. 201 Compute probability density dnorm(x,mu,sigma) function at x for a random variable from the Normal distribution with mean µ and standard deviation σ. 202 Compute probability density unifpdf(x,a,b) or dunif(x,a,b) function at x for a random ((x>=a)&&(x<=b))/(b-a) will variable from the continuous work even without the Statistics uniform distribution on inter- Toolbox. val (a, b). 203 Compute probability that a unidpdf(x,n) or ((x==floor(x)) ((x==round(x)) && (x >= 1) && random variable from the dis- && (x>=1)&&(x<=n))/n will work (x <= n))/n crete uniform distribution on even without the Statistics Toolbox, integers 1 . . . n has value x. as long as n is a positive integer. Note: one or more of the parameters in the above “*pdf” (Matlab) or “d*” (R) functions can be vectors, but they must be the same size. Scalars are promoted to arrays of the appropriate size.

D. Hiebeler, Matlab / R Reference No. 204 The corresponding CDF functions are below: Description Matlab Compute probability that a binocdf(x,n,p). Without the random variable from the Statistics Toolbox, as long Binomial(n, p) distribution is as n is a non-negative inless than or equal to x (i.e. teger, this will work: r = the cumulative distribution 0:floor(x); sum(factorial(n)./ function, or cdf). (factorial(r).*factorial(n-r)) .*p.^r.*(1-p).^(n-r)). (Unfortunately, Matlab ’s nchoosek function won’t take a vector argument for k.) Compute probability that a poisscdf(x,lambda). Withrandom variable from the out the Statistics Toolbox, as Poisson(λ) distribution is less long as lambda ≥ 0, this than or equal to x. will work: r = 0:floor(x); sum(exp(-lambda)*lambda.^r ./factorial(r)) Compute cumulative distri- expcdf(x,mu) or bution function at x for a (x>=0)*(1-exp(-x/mu)) will random variable from the ex- work even without the Statistics ponential distribution with Toolbox, as long as mu is positive. mean µ. Compute cumulative distri- normcdf(x,mu,sigma) or 1/2 bution function at x for a ran- erf(-(x-mu)/(sigma*sqrt(2)))/2 dom variable from the Nor- will work even without the Statismal distribution with mean µ tics Toolbox, as long as sigma is and standard deviation σ. positive. Compute cumulative distri- unifcdf(x,a,b) or bution function at x for a ran- (x>a)*(min(x,b)-a)/(b-a) will dom variable from the contin- work even without the Statistics uous uniform distribution on Toolbox, as long as b > a. interval (a, b). Compute probability that a unidcdf(x,n) or random variable from the dis- (x>=1)*min(floor(x),n)/n will crete uniform distribution on work even without the Statistics integers 1 . . . n is less than or Toolbox, as long as n is a positive equal to x. integer. R pbinom(x,n,p)

24

205

ppois(x,lambda)

206

pexp(x,1/mu)

207

pnorm(x,mu,sigma)

208

punif(x,a,b)

209

(x>=1)*min(floor(x),n)/n

D. Hiebeler, Matlab / R Reference

25

7
7.1
No. 210

Graphics
Various types of plotting
Matlab figure R windows() (when running R in Windows), quartz() (in Mac OS-X), or x11() (in Linux) dev.set(n) (returns the actual device selected; will be different from n if there is no figure device with number n) dev.list() dev.off() to close the currently active figure device, dev.off(n) to close a specified one, and graphics.off() to close all figure devices. plot(x,y) plot(x,y,type=’l’) (Note: that’s a lower-case ’L’, not the number 1) plot(x,y,type=str1, pch=arg2,col=str3, lty=arg4) See tables below for possible values of the 4 parameters plot(..., log=’x’), plot(..., log=’y’), and plot(..., log=’xy’) plot with logarithmic scales for x, y, and both axes, respectively Can’t do this in R; but barplot(y) makes a bar graph where you specify the heights, barplot(y,w) also specifies the widths of the bars, and hist can make plots like this too. Description Create a new figure window

211

Select figure number n

figure(n) (will create the figure if it doesn’t exist)

212 213

List open figure windows Close figure window(s)

214 215 216

Plot points using open circles Plot points using solid lines Plotting: color, point markers, linestyle

get(0,’children’) (The 0 handle refers to the root graphics object.) close to close the current figure window, close(n) to close a specified figure, and close all to close all figures plot(x,y,’o’) plot(x,y) plot(x,y,str) where str is a string specifying color, point marker, and/or linestyle (see table below) (e.g. ’gs--’ for green squares with dashed line)

217

Plotting axes

with

logarithmic

semilogx, semilogy, and loglog functions take arguments like plot, and plot with logarithmic scales for x, y, and both axes, respectively bar(x,y) Or just bar(y) if you only want to specify heights. Note: if A is a matrix, bar(A) interprets each column as a separate set of observations, and each row as a different observation within a set. So a 20 × 2 matrix is plotted as 2 sets of 20 observations, while a 2 × 20 matrix is plotted as 20 sets of 2 observations. hist(x) v=unique(x); c=hist(x,v); bar(v,c)

218

Make bar graph where the x coordinates of the bars are in x, and their heights are in y

219 220

Make histogram of values in x Given vector x containing integer values, make a bar graph where the x coordinates of bars are the values, and heights are the counts of how many times the values appear in x

hist(x) hist(x,(min(x)-.5):(max(x)+.5))

D. Hiebeler, Matlab / R Reference No. 221 Description Given vector x containing continuous values, lump the data into k bins and make a histogram / bar graph of the binned data Make a plot containing errorbars of height s above and below (x, y) points Make a plot containing errorbars of height a above and b below (x, y) points Other types of 2-D plots Matlab [c,m] = hist(x,k); bar(m,c) or for slightly different plot style use hist(x,k)

26 R hist(x,seq(min(x), max(x), length.out=k+1))

222

errorbar(x,y,s)

223

errorbar(x,y,b,a)

224

225

Make a 3-D plot of some data points with given x, y, z coordinates in the vectors x, y, and z. Surface plot of data in matrix A

stem(x,y) and stairs(x,y) for other types of 2-D plots. polar(theta,r) to use polar coordinates for plotting. plot3(x,y,z) This works much like plot, as far as plotting symbols, linetypes, and colors.

errbar(x,y,y+s,y-s) Note: errbar is part of the Hmisc package (see item 295 for how to install/load packages). errbar(x,y,y+a,y-b) Note: errbar is part of the Hmisc package (see item 295 for how to install/load packages). pie(v)

cloud(z~x*y) You can also use arguments pch and col as with plot. To make a 3-D plot with lines, do cloud(z~x*y,type=’l’, panel.cloud=panel.3dwire) persp(A) You can include shading in the image via e.g. persp(A,shade=0.5). There are two viewing angles you can also specify, among other parameters, e.g. persp(A, shade=0.5, theta=50, phi=35). x = seq(0,10,100) y = seq(2,8,90) f = function(x,y) return(sin(x+y)*sqrt(y)) z = outer(x,y,f) persp(x,y,z) contour(x,y,z) Or do s=expand.grid(x=x,y=y), and then wireframe(z~x*y,s) or wireframe(z~x*y,s,shade=TRUE) (Note: wireframe is part of the lattice package; see item 295 for how to load packages). If you have vectors x, y, and z all the same length, you can also do symbols(x,y,z).

226

surf(A) You can then click on the small curved arrow in the figure window (or choose “Rotate 3D” from the “Tools” menu), and then click and drag the mouse in the figure to rotate it in three dimensions.

227

Surface plot √ sin(x + y) y of x between 90 values of y 8

of f (x, y) = for 100 values 0 and 10, and between 2 and

x = linspace(0,10,100); y = linspace(2,8,90); [X,Y] = meshgrid(x,y); Z = sin(X+Y).*sqrt(Y); surf(X,Y,Z) shading flat mesh(X,Y,Z), surfc(X,Y,Z), surfl(X,Y,Z), contour(X,Y,Z), pcolor(X,Y,Z), waterfall(X,Y,Z). Also see the slice command.

228

Other ways of plotting the data from the previous command

D. Hiebeler, Matlab / R Reference No. 229 Adding various labels or making adjustments to plots Description Matlab Set axis ranges in a figure axis([x1 x2 y1 y2]) window

27 R You have to do this when you make the plot, e.g. plot(x,y,xlim=c(x1,x2), ylim=c(y1,y2)) title(main=’somestring’) adds a main title, title(sub=’somestring’) adds a subtitle. You can also include main= and sub= arguments in a plot command. title(xlab=’somestring’, ylab=’anotherstr’). You can also include xlab= and ylab= arguments in a plot command. plot(x,y,xlab= expression(phi^2 + mu[’i,j’])) or plot(x,y,xlab=expression( paste(’fecundity ’, phi))) See also help(plotmath) and p. 98 of the R Graphics book by Paul Murrell for more.

230

Add title to plot

title(’somestring’)

231

Add axis labels to plot

xlabel(’somestring’) ylabel(’somestring’)

and

232

Include Greek letters or symbols in plot axis labels

233

For on-screen graphics, do par(ps=16) followed by e.g. a plot command. For PostScript or PDF plots, add a pointsize=16 argument, e.g. pdf(’myfile.pdf’, width=8, height=8, pointsize=16) (see items 245 and 246) 234 Add grid lines to plot grid on (and grid off to turn off) grid() Note that if you’ll be printing the plot, the default style for grid-lines is to use gray dotted lines, which are almost invisible on some printers. You may want to do e.g. grid(lty=’dashed’, col=’black’) to use black dashed lines which are easier to see. 235 Add figure legend to top-left legend(’first’, ’second’, legend(’topleft’, corner of plot ’Location’, ’NorthWest’) legend=c(’first’, ’second’), col=c(’red’, ’blue’), pch=c(’*’,’o’)) Matlab note: sometimes you build a graph piece-by-piece, and then want to manually add a legend which doesn’t correspond with the order you put things in the plot. You can manually construct a legend by plotting “invisible” things, then building the legend using them. E.g. to make a legend with black stars and solid lines, and red circles and dashed lines: h1=plot(0,0,’k*-’); set(h1,’Visible’, ’off’); h2=plot(0,0,’k*-’); set(h2,’Visible’, ’off’); legend([h1 h2], ’blah, ’whoa’). Just be sure to choose coordinates for your “invisible” points within the current figure’s axis ranges.

Change font size to 16 in plot labels

You can use basic TeX commands, e.g. plot(x,y); xlabel(’\phi^2 + \mu_{i,j}’) or xlabel(’fecundity \phi’) See also help tex.m and parts of doc text props for more about building labels using general LaTeX commands For the legends and numerical axis labels, use set(gca, ’FontSize’, 16), and for text labels on axes do e.g. xlabel(’my x var’, ’FontSize’, 16)

D. Hiebeler, Matlab / R Reference No. 236 Description Adding more things to a figure Matlab hold on means everything plotted from now on in that figure window is added to what’s already there. hold off turns it off. clf clears the figure and turns off hold.

28 R points(...) and lines(...) work like plot, but add to what’s already in the figure rather than clearing the figure first. points and lines are basically identical, just with different default plotting styles. Note: axes are not recalculated/redrawn when adding more things to a figure. matplot(x,y) where x and y are 2-D matrices. Each column of x is plotted against the corresponding column of y. If x has only one column, it will be re-used. curve(sin(2*x), 7, 18, 200) makes the plot, by sampling the value of the function at 200 values between 7 and 18 (if you don’t specify the number of points, 101 is the default). You could do this manually yourself via commands like tmpx=seq(7,18,200); plot(tmpx, sin(2*tmpx)). image(A) (it rotates the matrix 90 degrees counterclockwise: it draws row 1 of A as the left column of the image, and column 1 of A as the bottom row of the image, so the row number is the x coord and column number is the y coord). It also rescales colors. If you are using a colormap with k entries, but the value k does not appear in A, use image(A,zlim=c(1,k)) to avoid rescaling of colors. Or e.g. image(A,zlim=c(0,k-1)) if you want values 0 through k−1 to be plotted using the k colors. Use filled.contour(A) rather than image(A), although it “blurs” the data via interpolation, or use levelplot(A) from the lattice package (see item 295 for how to load packages). To use a colormap with the latter, do e.g. levelplot(A,col.regions= terrain.colors(100)). image(A,col=terrain.colors(100)). The parameter 100 specifies the length of the colormap. Other colormaps are heat.colors(), topo.colors(), and cm.colors().

237

Plot multiple data sets at once

238

Plot sin(2x) for x between 7 and 18

plot(x,y) where x and y are 2-D matrices. Each column of x is plotted against the corresponding column of y. If x has only one column, it will be re-used. fplot(’sin(2*x)’, [7 18])

239

Plot color image of integer values in matrix A

image(A) to use array values as raw indices into colormap, or imagesc(A) to automatically scale values first (these both draw row 1 of the matrix at the top of the image); or pcolor(A) (draws row 1 of the matrix at the bottom of the image). After using pcolor, try the commands shading flat or shading interp.

240

Add colorbar legend to image plot

colorbar, pcolor.

after using image or

241

Set colormap in image

colormap(hot). Instead of hot, you can also use gray, flag, jet (the default), cool, bone, copper, pink, hsv, prism. By default, the length of the new colormap is the same as the currently-installed one; use e.g. colormap(hot(256)) to specify the number of entries.

D. Hiebeler, Matlab / R Reference No. 242 Description Build your own colormap using Red/Green/Blue triplets Matlab Use an n × 3 matrix; each row gives R,G,B intensities between 0 and 1. Can use as argument with colormap. E.g. for 2 colors: mycmap = [0.5 0.8 0.2 ; 0.2 0.2 0.7]

29 R Use a vector of hexadecimal strings, each beginning with ’#’ and giving R,G,B intensities between 00 and FF. E.g. c(’#80CC33’,’#3333B3’); can use as argument to col= parameter to image. You can build such a vector of strings from vectors of Red, Green, and Blue intensities (each between 0 and 1) as follows (for a 2-color example): r=c(0.5,0.2); g=c(0.8,0.2); b=c(0.2,0.7); mycolors=rgb(r,g,b).

Matlab plotting specifications, Symbol Color Symbol b blue . g green o r red x c cyan + m magenta * y yellow s k black d w white v ^ < > p h

for use with plot, fplot, semilogx, semilogy, loglog, etc: Marker Symbol Linestyle point (.) solid line circle (◦) : dotted line cross (×) -. dash-dot line plus sign (+) -dashed line asterisk (∗) square ( ) diamond (♦) triangle (down) (▽) triangle (up) (△) triangle (left) (⊳) triangle (right) (⊲) pentragram star hexagram star

R plotting specifications for col (color), pch (plotting character), and type arguments, for use with plot, matplot, points, and lines: col Description pch Description type Description ’blue’ Blue ’a’ a (similarly for other p points characters, but see ’.’ below for an exception ’green’ Green 19 solid circle l lines ’red’ Red 20 bullet (smaller circle) b both ’cyan’ Cyan 21 open circle c lines part only of “b” ’magenta’ Magenta 22 square o lines, points overplotted ’yellow’ Yellow 23 diamond h histogram-like lines ’black’ Black 24 triangle point-up s steps ’#RRGGBB’ hexadecimal specifica25 triangle point-down S another kind of steps tion of Red, Green, Blue (Other names) See colors() for list of ’.’ rectangle of size 0.01 n no plotting available color names. inch, 1 pixel, or 1 point (1/72 inch) depending on device (See table on next page for more)

D. Hiebeler, Matlab / R Reference

30

R plotting specifications for lty (line-type) argument, for use with plot, matplot, points, and lines: lty Description 0 blank 1 solid 2 dashed 3 dotted 4 dotdash 5 longdash 6 twodash

24

25

A A

b b

.

# #

18

19

20

21

22

23

12

13

14

15

16

17

6

7

8

9

10

11

0

1

2

3

4

5

R plotting characters, i.e. values for pch argument (from the book R Graphics, by Paul Murrell, Chapman & Hall / CRC, 2006)

D. Hiebeler, Matlab / R Reference No. 243 Description Divide up a figure window into smaller sub-figures Matlab subplot(m,n,k) divides the current figure window into an m × n array of subplots, and draws in subplot number k as numbered in “reading order,” i.e. left-to-right, top-tobottom. E.g. subplot(2,3,4) selects the first sub-figure in the second row of a 2 × 3 array of sub-figures. You can do more complex things, e.g. subplot(5,5,[1 2 6 7]) selects the first two subplots in the first row, and first two subplots in the second row, i.e. gives you a bigger subplot within a 5 × 5 array of subplots. (If you that command followed by e.g. subplot(5,5,3) you’ll see what’s meant by that.)

31 R There are several ways to do this, e.g. using layout or split.screen, although they aren’t quite as friendly as Matlab E.g. if you let A = ’s.  1 1 2  1 1 3 , then layout(A) will 4 5 6 divide the figure into 6 sub-figures: you can imagine the figure divide into a 3 × 3 matrix of smaller blocks; subfigure 1 will take up the upper-left 2 × 2 portion, and sub-figures 2–6 will take up smaller portions, according to the positions of those numbers in the matrix A. Consecutive plotting commands will draw into successive subfigures; there doesn’t seem to be a way to explicitly specify which sub-figure to draw into next. To use split.screen, you can do e.g. split.screen(c(2,1)) to split into a 2 × 1 matrix of subfigures (numbered 1 and 2). Then split.screen(c(1,3),2) splits subfigure 2 into a 1 × 3 matrix of smaller sub-figures (numbered 3, 4, and 5). screen(4) will then select sub-figure number 4, and subsequent plotting commands will draw into it. A third way to accomplish this is via the commands par(mfrow=) or par(mfcol=) to split the figure window, and par(mfg=) to select which sub-figure to draw into. Note that the above methods are all incompatible with each other. R automatically updates graphics windows even before functions/scripts finish executing, so it’s not necessary to explictly request it. But note that some graphics functions (particularly those in the lattice package) don’t display their results when called from scripts or functions; e.g. rather than levelplot(...) you need to do print(levelplot(...)). Such functions will automatically display their plots when called interactively from the command prompt.

244

Force graphics windows to update

drawnow (Matlab normally only updates figure windows when a script/function finishes and returns control to the Matlab prompt, or under a couple of other circumstances. This forces it to update figure windows to reflect any recent plotting commands.)

D. Hiebeler, Matlab / R Reference

32

7.2
No. 245

Printing/saving graphics
Matlab print -dpdf fname saves the contents of currently active figure window R First do pdf(’fname.pdf’). Then, do various plotting commands to make your image, as if you were plotting in a window. Finally, do dev.off() to close/save the PDF file. To print the contents of the active figure window, do dev.copy(device=pdf, file=’fname.pdf’); dev.off(). (But this will not work if you’ve turned off the display list via dev.control(displaylist= ’inhibit’).) postscript(’fname.eps’), followed by your plotting commands, followed by dev.off() to close/save the file. Note: you may want to use postscript(’fname.eps’, horizontal=FALSE) to save your figure in portrait mode rather than the default landscape mode. To print the contents of the active figure window, do dev.copy(device=postscript, file=’fname.eps’); dev.off(). (But this will not work if you’ve turned off the display list via dev.control(displaylist= ’inhibit’).) You can also include the horizontal=FALSE argument with dev.copy(). jpeg(’fname.jpg’,quality=90), followed by your plotting commands, followed by dev.off() to close/save the file.

Description To print/save to a PDF file named fname.pdf

246

To print/save to a PostScript file fname.ps or fname.eps

print -dps fname for black & white PostScript; print -dpsc fname for color PostScript; print -deps fname for black & white Encapsulated PostScript; print -depsc fname for color Encapsulated PostScript. The first two save to fname.ps, while the latter two save to fname.eps.

247

To print/save to a JPEG file fname.jpg with jpeg quality = 90 (higher quality looks better but makes the file larger)

print -djpeg90 fname

D. Hiebeler, Matlab / R Reference

33

7.3
No. 248

Animating cellular automata / lattice simulations
Matlab Repeatedly use either pcolor or image to display the data. Don’t forget to call drawnow as well, otherwise the figure window will not be updated with each image. R If you simply call image repeatedly, there is a great deal of flickering/flashing. To avoid this, after drawing the image for the first time using e.g. image(A), from then on only use image(A,add=TRUE), which avoids redrawing the entire image (and the associated flicker). However, this will soon consume a great deal of memory, as all drawn images are saved in the image buffer. There are two solutions to that problem: (1) every k time steps, leave off the “add=TRUE” argument to flush the image buffer (and get occasional flickering), where you choose k to balance the flickering vs. memory-usage tradeoff; or (2) after drawing the first image, do dev.control(displaylist= ’inhibit’) to prohibit retaining the data. However, the latter solution means that after the simulation is done, the figure window will not be redrawn if it is resized, or temporarily obscured by another window. (A call to dev.control(displaylist= ’enable’) and then one final image(A) at the end of the simulation will re-enable re-drawing after resizing or obscuring, without consuming extra memory.)

Description To display images of cellular automata or other lattice simulations while running in real time

D. Hiebeler, Matlab / R Reference

34

8
No. 249 250 251 252 253

Working with files
Description Create a folder (also known as a “directory”) Set/change working directory See list of files in current working directory Run commands in file ‘foo.m’ or ‘foo.R’ respectively Read data from text file “data.txt” into matrix A Matlab mkdir dirname cd dirname dir foo A=load(’data.txt’) or A=importdata(’data.txt’) Note that both routines will ignore comments (anything on a line following a “%” character) R dir.create(’dirname’) setwd(’dirname’) dir() source(’foo.R’) A=as.matrix(read.table( ’data.txt’)) This will ignore comments (anything on a line following a “#” character). To ignore comments indicated by “%”, do A=as.matrix(read.table( ’data.txt’, comment.char=’%’)) write(A, file=’data.txt’, ncolumn=dim(A)[2])

254

Write data from matrix A into text file “data.txt”

save data.txt A -ascii

D. Hiebeler, Matlab / R Reference

35

9
9.1
No. 255 256

Miscellaneous
Variables
Matlab x = 5 assignin(’base’, ’y’, 7) R x <- 5 or x = 5 y <<- 7 Description Assigning to variables From within a function, assign a value to variable y in the base environment (i.e. the command prompt environment) From within a function, access the value of variable y in the base environment (i.e. the command prompt environment) Short list of defined variables Long list of defined variables See detailed info about the variable ab See detailed info about all variables with “ab” in their name Open graphical data editor, to edit the value of variable A (useful for editing values in a matrix, though it works for non-matrix variables as well) Clear one variable Clear two variables Clear all variables See what type of object x is (Variable names)

257

evalin(’base’, ’y’)

y (In R, if there isn’t a local variable y within the function, it will look for one in the base environment.)

258 259 260 261

who whos whos ab whos *ab*

ls() ls.str() str(ab) ls.str(pattern=’ab’)

262

openvar(A), or double-click on the variable in the Workspace pane (if it’s being displayed) of your Matlabdesktop clear x clear x y clear all class(x) Variable names must begin with a letter, but after that they may contain any combination of letters, digits, and the underscore character. Names are case-sensitive. ans contains the result of the last command which did not assign its value to a variable. E.g. after 2+5; x=3, then ans will contain 7.

fix(A)

263 264 265 266 267

268

Result of last command

rm(x) rm(x,y) rm(list=ls()) class(x) Variable names may contain letters, digits, the period, and the underscore character. They cannot begin with a digit or underscore, or with a period followed by a digit. Names are casesensitive. .Last.value contains the result of the last command, whether or not its value was assigned to a variable. E.g. after 2+5; x=3, then .Last.value will contain 3.

D. Hiebeler, Matlab / R Reference

36

9.2
No. 269

Strings and Misc.
Matlab If you want to break up a Matlab command over more than one line, end all but the last line with three periods: “...”. E.g.: x = 3 + ... 4 of format short g format long g are help format and see R In R, you can spread commands out over multiple lines, and nothing extra is necessary. R will continue reading input until the command is complete. E.g.: x = 3 + 4 options(digits=6) tells R you’d like to use 6 digits of precision in values it displays (it is only a suggestion, not strictly followed) q() or quit() # this is a comment print(’hi there’)

Description Line continuation

270

Controlling output

formatting

handy;

271 272 273

Exit the program Comments Print a string

274

Print a string containing single quotes Give prompt and read input from user Concatenate strings Concatenate strings stored in a vector

275 276 277

quit or exit % this is a comment disp(’hi there’) or to omit trailing newline use fprintf(’hi there’) disp(’It’’s nice’) or to omit trailing newline fprintf(’It’’s nice’) x = input(’Enter data:’) [’two hal’ ’ves’] v={’two ’, ’halves’}; strcat(v{:}) But note that this drops trailing spaces on strings. To avoid that, instead do strcat([v{:}]) text1=’hi there’; text2=text(2:6) x = ’a’, ’aa’, ’bc’, ’c’; y = ’da’, ’a’, ’bc’, ’a’, ’bc’, ’aa’; [tf, loc]=ismember(x,y) Then loc contains the locations of last occurrences of elements of x in the set y, and 0 for unmatched elements. num2str(x)

print(’It\’s nice’) print("It’s nice")

or

print(’Enter data:’) x = scan() paste(’two hal’, ’ves’, sep=’’) v=c(’two ’, ’halves’); paste(v, collapse=’’)

278 279

Extract substring of a string Determine whether elements of a vector are in a set, and give positions of corresponding elements in the set.

text1=’hi there’; text2=substr(text1,2,6) x = c(’a’, ’aa’, ’bc’, ’c’); y = c(’da’, ’a’, ’bc’, ’a’, ’bc’, ’aa’); loc=match(x,y) Then loc contains the locations of first occurences of elements of x in the set y, and NA for unmatched elements. as.character(x)

280

Convert number to string

D. Hiebeler, Matlab / R Reference No. 281 Description Use sprintf to create a formatted string. Use %d for integers (“d” stands for “decimal”, i.e. base 10), %f for floating-point numbers, %e for scientific-notation floating point, %g to automatically choose %e or %f based on the value. You can specify field-widths/precisions, e.g. %5d for integers with padding to 5 spaces, or %.7f for floating-point with 7 digits of precision. There are many other options too; see the docs. Machine epsilon ǫmach , i.e. difference between 1 and the next largest double-precision floating-point number Pause for x seconds Wait for user to press any key Matlab x=2; y=3.5; s=sprintf(’x is %d, y=%g’, ... x, y) R

37

x=2; y=3.5 s=sprintf(’x is %d, y is %g’, x, y)

282

eps (See help eps for various other things eps can give.)

.Machine$double.eps

283 284

pause(x) pause

285 286

287 288 289

Measure CPU time used to do some commands Measure elapsed (“wallclock”) time used to do some commands Print an error message an interrupt execution Print a warning message Putting multiple statements on one line

t1=cputime; ...commands... ; cputime-t1 tic; ...commands... ; toc or t1=clock; ...commands... ; etime(clock,t1) error(’Problem!’) warning(’Smaller problem!’) Separate statements by commas or semicolons. A semicolon at the end of a statement suppresses display of the results (also useful even with just a single statement on a line), while a comma does not. eval(s) which sqrt shows you where the file defining the sqrt function is (but note that many basic functions are “built in,” so the Matlab function file is really just a stub containing documentation). This is useful if a command is doing something strange, e.g. sqrt isn’t working. If you’ve accidentally defined a variable called sqrt, then which sqrt will tell you, so you can clear sqrt to erase it so that you can go back to using the function sqrt.

Sys.sleep(x) Don’t know of a way to do this in R, but scan(quiet=TRUE) will wait until the user presses the Enter key t1=proc.time(); ...commands... ; (proc.time()-t1)[1] t1=proc.time(); ...commands... ; (proc.time()-t1)[3] stop(’Problem!’) warning(’Smaller problem!’) Separate statements by semicolons.

290 291

Evaluate contents of a string s as command(s). Show where a command is

eval(parse(text=s)) R does not execute commands directly from files, so there is no equivalent command.

D. Hiebeler, Matlab / R Reference No. 292 Description Query/set the search path. Matlab path displays the current search path (the list of places Matlab searches for commands you enter). To add a directory ~/foo to the beginning of the search path, do addpath ~/foo -begin or to add it to the end of the path, do addpath ~/foo -end (Note: you should generally add the full path of a directory, i.e. in Linux or Mac OS-X something like ~/foo as above or of the form /usr/local/lib/foo, while under Windows it would be something like C:/foo) If a file startup.m exists in the startup directory for Matlab, its contents are executed. (See the Matlab docs for how to change the startup directory.)

38 R R does not use a search path to look for files.

293

Startup sequence

294

Shutdown sequence

Upon typing quit or exit, Matlab will run the script finish.m if present somewhere in the search path. Matlab does not have packages. It has toolboxes, which you can purchase and install. “Contributed” code (written by end users) can simply be downloaded and put in a directory which you then add to Matlab’s path (see item 292 for how to add things to Matlab’s path).

295

Install and load a package.

If a file .Rprofile exists in the current directory or the user’s home directory (in that order), its contents are sourced; saved data from the file .RData (if it exists) are then loaded. If a function .First() has been defined, it is then called (so the obvious place to define this function is in your .Rprofile file). Upon typing q() or quit(), R will call the function .Last() if it has been defined (one obvious place to define it would be in the .Rprofile file) To install e.g. the deSolve package, you can use the command install.packages(’deSolve’). You then need to load the package in order to use it, via the command library(’deSolve’). When running R again later you’ll need to load the package again to use it, but you should not need to re-install it. Note that the lattice package is typically included with binary distributions of R, so it only needs to be loaded, not installed.

D. Hiebeler, Matlab / R Reference

39

10
No. 296

Spatial Modeling
Matlab A = (A | (rand(L) < p))*1; R A = (A | (matrix(runif(L^2),L) < p))*1

297

298

299

Description Take an L×L matrix A of 0s and 1s, and “seed” fraction p of the 0s (turn them into 1s), not changing entries which are already 1. Take an L × L matrix A of 0s and 1s, and “kill” fraction p of the 1s (turn them into 0s), not changing the rest of the entries Do “wraparound” on a coordinate newx that you’ve already calculated. You can replace newx with x+dx if you want to do wraparound on an offset x coordinate. Randomly initialize a portion of an array: set fraction p of sites in rows iy1 through iy2 and columns ix1 through ix2 equal to 1 (and set the rest of the sites in that block equal to zero). Note: this assume iy1 < iy2 and ix1 < ix2.

A = (A & (rand(L) < 1-p))*1;

A = (A & (matrix(runif(L^2),L) < 1-p))*1

mod(newx-1,L)+1 Note: for portability with other languages such as C which handle MOD of negative values differently, you may want to get in the habit of instead doing mod(newx-1+L,L)+1 dx=ix2-ix1+1; dy=iy2-iy1+1; A(iy1:iy2,ix1:ix2) = ... (rand(dy,dx) < p0)*1;

((newx-1) %% L) + 1 Note: for portability with other languages such as C which handle MOD of negative values differently, you may want to get in the habit of instead doing ((newx-1+L)%%L) + 1 dx=ix2-ix1+1; dy=iy2-iy1+1; A[iy1:iy2,ix1:ix2] = (matrix(runif(dy*dx),dy) < p0)*1

INDEX OF MATLAB COMMANDS AND CONCEPTS

40

Index of MATLAB commands and concepts
’, 72 ,, 289 .*, 71 ..., 269 ./, 77 .^, 81 /, 76 :, 12–14 ;, 289 =, 255 [, 6–8 %, 272 &, 165, 166 ^, 46, 79, 80 \, 73, 78 { 41 abs, 47, 65 acos, 52 acosh, 54 addpath, 292 all, 167 angle, 66 ans, 268 any, 168 asin, 52 asinh, 54 assignin, 256 atan, 52 atanh, 54 average, see mean axis, 229 bar, 218, 220, 221 binocdf, 204 binopdf, 198 binornd, 191 boolean tests scalar, 165 vector, 166–168 cd, 250 ceil, 58 cell, 40 cell arrays, 40 extracting elements of, 41 cellular automata animation, 248 chol, 87 class, 266 clear, 263–265 clf, 236 clock, 286 close, 213 colon, see : colorbar, 240 colormap building your own, 242 colormap, 241, 242 column vector, 7 comments, 272 complex numbers, 64–69 cond, 91–93 conj, 67 contour, 228 conv, 145 corr, 105–110 cos, 51 cosh, 53 cov, 103, 104 cputime, 285 csape, 157, 159, 160 cubic splines, 158, 159 natural, 157 not-a-knot, 161 periodic, 160 cumprod, 119 cumsum, 115–118 cumulative distribution functions binomial, 204 continuous uniform on interval (a, b), 208 discrete uniform from 1..n, 209 exponential, 206 normal, 207 Poisson, 205 diag, 21, 22 diff, 121 differential equations, see ode45 dir, 251 disp, 273, 274 doc, 4 drawnow, 244, 248 echelon form, see matrix eig, 83 element-by-element matrix operations, see matrix else, 164 elseif, 164 eps, 282 erf, 60

INDEX OF MATLAB COMMANDS AND CONCEPTS erfc, 61 erfcinv, 63 erfinv, 62 error, 287 errorbar, 222, 223 etime, 286 eval, 290 evalin, 257 exit, 271 exp, 48 expcdf, 206 expm, 114 exppdf, 200 exprnd, 193 eye, 20 figure, 210, 211 file reading data from, 254 running commands in, 252 text reading data from, 253 saving data to, 254 find, 140–142 finish.m, 294 floor, 57 fminbnd, 148, 149 fminsearch, 150, 151 font size in plots, 233 for, 162 format, 270 fplot, 238 fprintf, 273, 274 function multi-variable minimization, 150 minimization over first parameter only, 149 minimization over only some parameters, 151 single-variable minimization, 148 user-written, 171 returning multiple values, 172 fzero, 147 gca, 233 get, 212 Greek letters in plot labels, 232 grid, 234 help, 1–3 helpbrowser, 4 helpdesk, 4 hilb, 38 hist, 143, 144, 219, 220 hold, 236 identity, see matrix if, 163–165 imag, 69 image, 239, 248 imagesc, 239 importdata, 253 ind2sub, 31 indexing matrix, 10 with a single index, 11 vector, 9 input, 275 inv, 75 inverse, see matrix ismember, 279 legend, 235 length, 134, 136 linspace, 15 load, 253, 254 log, 49 log10, 50 log2, 50 loglog, 217 lookfor, 5 lu, 84

41

matrix, 8 boolean operations on, 141, 142 changing shape of, 35 Cholesky factorization, 87 condition number, 91–93 containing all indentical entries, 19 containing all zeros, 18 converting row, column to single index, 32 converting single-index to row, column, 31 cumulative sums of all elements of, 118 cumulative sums of columns, 116 cumulative sums of rows, 117 diagonal, 21 echelon form, 74 eigenvalues and eigenvectors of, 83 equation solving, 73 exponential of, 114 extracting a column of, 26 extracting a rectangular piece of, 29 extracting a row of, 27 extracting specified rows and columns of, 30 “gluing” together, 23, 24

INDEX OF MATLAB COMMANDS AND CONCEPTS identity, 20 inverse, 75 lower-triangular portion of, 36 LU factorization, 84 minimum of values of, 124 minimum value of each column of, 125 minimum value of each row of, 126 modifying elements given lists of rows and columns, 33 multiplication, 70 element-by-element, 71 N -dimensional, 39 norm, 90 powers of, 80 QR factorization, 88 rank, 82 re-shaping its elements into a vector, 28 Schur decomposition, 86 singular value decomposition, 85 size of, 131–133, 135, 136 sum of all elements, 111 of columns of, 112 of rows of, 113 transpose, 72 upper-triangular portion of, 37 max, see min mean, 94–96 mesh, 228 meshgrid, 105 min, 123–126, 128–130 mind, 127 mkdir, 249 mod, 55, 298 modulo arithmetic, 55, 298 multiple statements on one line, 289 norm, 89, 90 normcdf, 207 normpdf, 201 normrnd, 197 num2str, 280 numel, 135 ode45, 173–175 ones, 17, 19 openvar, 262 optimization, 148–151 path, 292 pause, 283, 284 pcolor, 228, 239, 248 perform some commands with probability p, 185 permutation of integers 1..n, 186

42

plot, 214–216, 237 Greek letters in axis labels, 232 plot3, 225 poisscdf, 205 poisspdf, 199 poissrnd, 192 polar, 224 polyfit, 153–155 polynomial least-squares fitted, 154–156 multiplication, 145 roots of, 146 ppval, 157, 159, 160 print, 245–247 probability density functions binomial, 198 continuous uniform on interval (a, b), 202 discrete uniform from 1..n, 203 exponential, 200 normal, 201 Poisson, 199 qr, 88 quad, 152 quit, 271 rand, 176–184, 190 random values Bernoulli, 182 binomial, 191 continuous uniform distribution on interval (a, b), 179, 196 continuous uniform distribution on interval (0,1), 176–178 discrete uniform distribution from a..b, 184 discrete uniform distribution from 1..k, 181, 194, 195 discrete uniform distribution, 180 exponential, 193 k unique values sampled from integers 1..n, 187 normal, 197 Poisson, 192 setting the seed, 190 randperm, 186, 187 randsample, 187–189 rank, 82 rcond, 91 real, 68 reshape, 35, 39 roots of general single-variable function, 147 polynomial, 146

INDEX OF MATLAB COMMANDS AND CONCEPTS roots, 146 round, 56 row vector, 6 rref, 74 sampling values from a vector, 188, 189 save, 254 schur, 86 semilogx, 217 semilogy, 217 set, 233 sign, 59 sin, 51 sinh, 53 size, 131–133 slice, 228 sort, 137, 138, 187 spline, 161 splines, see cubic splines sprintf, 281 sqrt, 45 stairs, 224 standard deviation, see std startup.m, 293 std, 97–99 stem, 224 stop, 287 strcat, 277 string concatenation, 276 converting number to, 280 substrings, 278 struct, 43 sub2ind, 32, 33 subplot, 243 sum, 111–113, 166 surf, 226, 227 surfc, 228 surfl, 228 svd, 85 switch, 170 tan, 51 tanh, 53 tic, 286 title, 230 toc, 286 transpose, see matrix tril, 36 triu, 37 unidcdf, 209 unidpdf, 203 unidrnd, 194, 195 unifcdf, 208 unifpdf, 202 unifrnd, 196 unique, 143, 220

43

var, 100–102 variables assigning, 255 assigning in base environment from function, 256 evaluating from base environment within function, 257 names, 267 variance, see var vector boolean operations on, 139, 140 containing all indentical entries, 17 containing all zeros, 16 counts of binned values in, 144 counts of discrete values in, 143 cumulative sum of elements of, 115 differences between consecutive elements of, 121 minimum of values of, 123 norm, 89 position of first occurance of minimum value in, 130 reversing order of elements in, 25 size of, 134 sum of all elements, 111 truncating, 34 warning, 288 waterfall, 228 which, 291 while, 169 who, 258 whos, 259–261 xlabel, 231–233 ylabel, 231, 232 zeros, 16, 18

INDEX OF R COMMANDS AND CONCEPTS

44

Index of R commands and concepts
*, 79 /, 77 :, 12, 13 ;, 289 <-, 255 <<-, 256 =, 255 ?, 1, 2 [[, 41 #, 272 %%, 55, 298 &, 165, 166 ^, 46, 81 abs, 47, 65 acos, 52 acosh, 54 all, 167 any, 168 apply, 99, 101, 102, 125, 126 Arg, 66 array, 39 as.character, 280 as.numeric, 143 asin, 52 asinh, 54 atan, 52 atanh, 54 average, see mean barplot, 218 boolean tests scalar, 165 vector, 166–168 c, 6, 7 cbind, 23, 33 ceiling, 58 cellular automata animation, 248 chol, 87 class, 266 cloud, 225 coef, 153, 154, 156 colMeans, 95 colon, see : colormap building your own, 242 for image, 241 colSums, 112 column vector, 7 comments, 272 complex numbers, 64–69 Conj, 67 contour, 228 convolve, 145 cor, 106–110 cos, 51 cosh, 53 cov, 103–105 cubic splines, 158, 159, 161 natural, 157 periodic, 160 cummax, 120 cummin, 120 cumprod, 119 cumsum, 115–118 cumulative distribution functions binomial, 204 continuous uniform on interval (a, b), 208 discrete uniform from 1..n, 209 exponential, 206 normal, 207 Poisson, 205 curve, 238 data.frame, 43 dbinom, 198 dev.control, 245, 246, 248 dev.copy, 245, 246 dev.list, 212 dev.off, 213, 245–247 dev.set, 211 dexp, 200 diag, 20–22 diff, 121 differential equations, see lsoda dim, 35, 133, 136 dir, 251 dir.create, 249 dnorm, 201 dpois, 199 dunif, 202 echelon form, see matrix eig, 83 element-by-element matrix operations, see matrix else, 164 errbar, 222, 223 eval, 290 exp, 48

INDEX OF R COMMANDS AND CONCEPTS expand, 84 expand.grid, 228 expm, 114 file reading data from, 254 running commands in, 252 text reading data from, 253 saving data to, 254 filled.contour, 240 .First, 293 fix, 262 floor, 57 font size in plots, 233 for, 162 function multi-variable minimization, 150 minimization over first parameter only, 149 minimization over only some parameters, 151 single-variable minimization, 148 user-written, 171 returning multiple values, 172 graphics not being displayed from scripts/functions, 244 Greek letters in plot labels, 232 grid, 234 help, 1, 2 help.search, 5 help.start, 4 Hilbert, 38 hist, 144, 218–221 identity, see matrix if, 163–165 ifelse, 122 Im, 69 image, 239, 248 indexing matrix, 10 with a single index, 11 vector, 9 install.packages, 295 integrate, 152 inverse, see matrix jpeg, 247 kappa, 92 .Last, 294 .Last.value, 268 lattice package, 228, 240, 244, 295 layout, 243 legend, 235 length, 34, 134, 135 levelplot, 240, 244 library, 3, 295 lines, 236 lists, 40 extracting elements of, 41 lm, 153, 154, 156 log, 49 log10, 50 log2, 50 lower.tri, 37 ls, 258 ls.str, 259, 261 lsoda, 173–175

45

.Machine$double.eps, 282 match, 279 matplot, 237 matrix, 8 boolean operations on, 141, 142 changing shape of, 35 Cholesky factorization, 87 condition number, 91–93 containing all indentical entries, 19 containing all zeros, 18 converting row, column to single index, 32 converting single-index to row, column, 31 cumulative sums of all elements of, 118 cumulative sums of columns, 116 cumulative sums of rows, 117 diagonal, 21 echelon form, 74 eigenvalues and eigenvectors of, 83 equation solving, 73 exponential of, 114 extracting a column of, 26 extracting a rectangular piece of, 29 extracting a row of, 27 extracting specified rows and columns of, 30 “gluing” together, 23, 24 identity, 20 inverse, 75 lower-triangular portion of, 36 LU factorization, 84 minimum of values of, 124

INDEX OF R COMMANDS AND CONCEPTS minimum value of each column of, 125 minimum value of each row of, 126 modifying elements given lists of rows and columns, 33 multiplication, 70 element-by-element, 71 N -dimensional, 39 norm, 90 powers of, 80 QR factorization, 88 rank, 82 re-shaping its elements into a vector, 28 Schur decomposition, 86 singular value decomposition, 85 size of, 131–133, 135, 136 sum of all elements, 111 of columns of, 112 of rows of, 113 transpose, 72 upper-triangular portion of, 37 matrix, 8, 18, 19 max, see min mean, 94 min, 123–126, 129 Mod, 65 modulo arithmetic, 55, 298 multiple statements on one line, 289 names, 42, 143 ncol, 132 norm, 89, 90 nrow, 131 optim, 150, 151 optimization, 148–151 optimize, 148, 149 options digits=, 270 outer, 227 packages installing, 295 loading, 295 par, 233 par mfcol=, 243 mfrow=, 243 parse, 290 paste, 276, 277 pbinom, 204 pdf, 233, 245 perform some commands with probability p, 185 permutation of integers 1..n, 186

46 persp, 226, 227 pexp, 206 pie, 224 plot, 214–217 Greek letters in axis labels, 232 main=, 230 sub=, 230 xlab=, 231, 232 xlim=, 229 ylab=, 231, 232 ylim=, 229 pmin, 127, 128 pnorm, 60, 61, 207 points, 236 polynomial least-squares fitted, 154–156 multiplication, 145 roots of, 146 polyreg, 155 polyroot, 146 postscript, 246 ppois, 205 print, 244, 273, 274 probability density functions binomial, 198 continuous uniform on interval (a, b), 202 discrete uniform from 1..n, 203 exponential, 200 normal, 201 Poisson, 199 proc.time, 285, 286 punif, 208 q, 271 qnorm, 62, 63 qr, 82, 88 quartz, 210 quit, 271 rand, 183 random values Bernoulli, 182 binomial, 191 continuous uniform distribution on interval (a, b), 179, 196 continuous uniform distribution on interval (0,1), 176, 178 continuous uniform distribution on inteval (0,1), 177 discrete uniform distribution from a..b, 184 discrete uniform distribution from 1..k, 181, 194, 195 discrete uniform distribution, 180

INDEX OF R COMMANDS AND CONCEPTS exponential, 193 k unique values sampled from integers 1..n, 187 normal, 197 Poisson, 192 setting the seed, 190 rbind, 24 rbinom, 191 rcond, 91, 93 .RData, 293 Re, 68 read.table, 253, 254 rep, 16, 17 rev, 25 rexp, 193 rgb, 242 rm, 263–265 rnorm, 197 roots of general single-variable function, 147 polynomial, 146 round, 56 row vector, 6 rowMeans, 96 rpois, 192 .Rprofile, 293 runif, 176–182, 184, 196 sample, 186–189, 194, 195 sampling values from a vector, 188, 189 scan, 275, 284 Schur, 86 sd, 97–99 seq, 14, 15 set.seed, 190 setwd, 250 sign, 59 sin, 51 sinh, 53 solve, 73, 75, 76, 78 sort, 137, 138 source, 252 spline, 157, 158, 160 splines, see cubic splines split.screen, 243 sprintf, 281 sqrt, 45 standard deviation, see sd str, 260 string concatenation, 276 converting number to, 280 substrings, 278 substr, 278 sum, 111, 113, 166 svd, 85 switch, 170 symbols, 228 Sys.sleep, 283 t, 72 table, 143 tan, 51 tanh, 53 title, 230, 231 transpose, see matrix uniroot, 147 upper.tri, 36

47

var, 100–102, 104 variables assigning, 255 assigning in base environment from function, 256 evaluating from base environment within function, 257 names, 267 variance, see var vector boolean operations on, 139, 140 containing all indentical entries, 17 containing all zeros, 16 counts of binned values in, 144 counts of discrete values in, 143 cumulative sum of elements of, 115 differences between consecutive elements of, 121 minimum of values of, 123 norm, 89 position of first occurance of minimum value in, 130 reversing order of elements in, 25 size of, 134 sum of all elements, 111 truncating, 34 vector, 40 warning, 288 which, 140–142 which.max, see which.min which.min, 130 while, 169 windows, 210 wireframe, 228 write, 254 x11, 210

 
statystyka