Note: you only need to hand-in your R/Splus printouts for Tasks 1 to 5.
Basic matrix operations in R/Splus:
A matrix is a rectangular array of numbers. Its dimension is described as “(number
of horizontal rows) by (number of vertical columns)”. Elements of a matrix are specified in
row/column coordinates. For example, in
a 5x4 matrix, the element x[5,4] is in
the 5th row, 4th column, which is the bottom right-hand
corner; x[5,1] is in the 5th row, 1st column,
which is the bottom left-hand corner; and so on. A matrix is square if the number of
rows and number of columns are equal. A
1x1 matrix is a single number, usually called a scalar. A 1xn or nx1 matrix is called a
vector.
A zero matrix is a matrix in which each element is
zero. Two matrices are equal if they
have the same dimensions and the corresponding elements in each matrix are
equal.
R/Splus has a number of facilities
for handling matrices. Some knowledge of
manipulating matrices will be necessary as we begin the study of Markov chains.
Task 1. Create the following matrices.
a). A 5x5 zero matrix: Z<-matrix(0,5,5)
b). Some 2x2 matrices, using
various input techniques:
A<-cbind(c(0,0),c(1,2))
B<-rbind(c(1,1),c(3,4))
D<-matrix(c(2,3,5,4),2,2)
E<-matrix(0,2,2)
E[1,1]<-3
E[1,2]<-7
d). Miscellaneous matrices:
one<-rbind(c(1,1))
two<-rbind(1:4,2:5)
identity<-diag(c(1,1))
e). Check the dimensions of a matrix: dim(A)
Task 2. Matrix arithmetic.
a). Scalar multiplication. What is 3*B?
b). Matrix addition.
What is Z+A? What is B+A?
What is A+B? What is A-B?
c). Matrix multiplication (denoted %*% in R/Splus). Compute the following:
Z%*%A
A%*%B
one%*%B
identity%*%B
identity%*%D
identity%*%E
B%*%two
two%*%B
Task 3. Recall that for scalars
a and e, ae=0
implies that a=0 or e=0.
Also, the cancellation law says that if a is
not equal to zero, then ab=ad implies
that b=d. The commutative law
for multiplication says that ab=ba. Consider the matrices you created in Task 1.
a). What is AE? How does this compare to the scalar result?
b). Does AB=AD? How does this compare to the scalar result?
c). Does AB=BA? How
does this compare to the scalar result?
Task 4. Solving
systems of linear equations. A system of linear equations can be rewritten
in matrix form to facilitate its solution.
For example, the system
x+y+2z=9
2x+4y-3z=1
3x+6y-5z=0
can be rewritten as
G (x,y,z)’=b
where the apostrophe ’ is
read transpose. The transpose
converts the 1x3 vector into a 3x1 vector.
The matrix G in the above expression is
G<-rbind(c(1,1,2),c(2,4,-3),c(3,6,-5))
and the vector b is
b<-cbind(c(9,1,0))
We can solve for (x,y,z)’
in the above equation by computing the inverse of the matrix G and
multiplying this inverse by b. (A
matrix must be square in order to have this kind of inverse, but not all square
matrices have inverses.) This is
analogous to dividing in the scalar case.
The inverse of the matrix G is computed in R/Splus
as
solve(G).
What is solve(G)%*%G? What is G%*%solve(G)? Find the solution
to the above system of equations. What
is solve(A) where A is defined in Task 1?
Task 5. The Weather Example discussed in class (see also Example 11.1, page 406 of text).
a). Set up the transition matrix P in R/Splus.
b). Given that today (Friday) is “Rain”, what is the probability that it will be “Snow” next Wednesday?
c). Thomas has been locked in a room with no windows for 10
days and he has no idea of today’s weather.
He guesses that, for today’s weather, the probabilities for “Rain”,
“Nice” and “Snow” are 0.2, 0.5 and 0.3 respectively. Based on his guess, what is the chance that
next Wednesday is “Nice”?
Task 6. A stochastic
process is a collection of random variables, usually indexed by time. It can be one-dimensional {X(1),X(2),X(3),…},
two-dimensional {(X(1),Y(1)), (X(2),Y(2)), (X(3),Y(3)), …}, or
multi-dimensional. Usually there is
dependence in the random variables from one time step to the next. Try out the following functions. Each uses randomly selected matrix operations
to simulate a two-dimensional stochastic process (X(t),Y(t))
and plot it in the x-y plane. The following two files can be found in
the subfolder st321/Spring2004 in drive G:, filename
“labsession7”.
First function:
“triangle”<-function(n = 20, p = c(1/3, 1/3)) {
if(sum(p) < 1) {
xy
<- matrix(1, 2, n)
a1 <- cbind(c(0.5, 0), c(0, 0.5))
b1 <- cbind(c(0, 0))
b2 <- cbind(c(0.5, 0))
b3 <- cbind(c(0.25,0.5))
for(i in 2:n) {
u <- runif(1)
if(u < p[1])
{
xy[, i] <- a1 %*% cbind(xy[, i
- 1]) + b1
}
if(u > p[1]
& u < sum(p)) {
xy[, i] <- a1 %*% cbind(xy[, i
- 1]) + b2
}
if(u >
sum(p)) {
xy[, i] <- a1 %*% cbind(xy[, i
- 1]) + b3
}
}
plot(xy[1, ], xy[2, ], xlab = "", ylab =
"", axes = F, pch = 17, col
= 8)
}
else{
cat(“sum(p)is
greater than one”)
}
return()
}
Second function:
"fern"<-function(n
= 20) {
xy <- matrix(1, 2, n)
a1 <- cbind(c(0, 0), c(0, 0.25))
b1 <- cbind(c(0, 0))
a2 <- cbind(c(0.85, -0.04), c(0.04,
0.85))
b2 <- cbind(c(0, 1.6))
a3 <- cbind(c(0.2, 0.26), c(-0.26,
0.22))
b3 <- cbind(c(0, 0.8))
a4 <- cbind(c(-0.15, 0.26), c(0.28,
0.24))
b4 <- cbind(c(0, 1))
for(i in 2:n) {
u <- runif(1)
if(u < 0.01) {
xy[, i] <- a1 %*% cbind(xy[, i - 1]) + b1
}
if(u > 0.01 & u < 0.86) {
xy[, i] <- a2 %*% cbind(xy[, i - 1]) + b2
}
if(u > 0.86 & u < 0.93) {
xy[, i] <- a3 %*% cbind(xy[, i - 1]) + b3
}
if(u > 0.93) {
xy[, i] <- a4 %*% cbind(xy[, i - 1]) + b4
}
}#end for loop
plot(xy[1, ], xy[2, ], xlab =
"", ylab = "", axes = F, pch = 18, col = 4)
return()
}
Task 7. It will be convenient if we have a function in R/Splus that will raise the power of a matrix. The following function, also in the file “labsession7”, executes this task. Copy this function and try the following commands.
mtxpow<-function(X,n) {
# raise the
matrix X to power n
# assume n>=2
temp<-X
for (i in 1:(n-1)) {
temp<-temp%*%X
}
return(temp)
}
mtxpow(A,2)
mtxpow(A,5)