Software tutorial/Vectors and matrices

From Statistics for Engineering
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
← Programming in R: loops and flow control (previous step) Tutorial index Next step: Building a least squares model in R →


<rst> <rst-options: 'toc' = False/> <rst-options: 'reset-figures' = False/> Creating matrices and vectors

An vector of zeros is created using the ``numeric`` function. This is used, for example, to initialize a vector of zeros before calling a :ref:`for loop <r-programming-loops-for-loop>` to fill the entries in that vector. This is called *preallocation*, and is good programming practice in any language: rather preallocate the matrix for storage than incrementally add entries to an existing matrix.

.. code-block:: s

numeric(6) [1] 0 0 0 0 0 0

A matrix or vector can be created using the ``matrix`` function. The advantage of the ``matrix`` function is that it can set the values to a non-zero, or missing values (``NA``):

.. code-block:: s

# Matrix of zeros: matrix(data=0, nrow=2, ncol=3) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0

# Column vector of missing values: matrix(nrow=4) [,1] [1,] NA [2,] NA [3,] NA [4,] NA

# Row vector of seven 7's matrix(data=7, ncol=7) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 7 7 7 7 7 7 7

I prefer to allocate matrices with missing entries, ``NA``, because it is easy to check after the loop if you have missed any entries. For example, for matrix ``X``, we can check: ``any(is.na(X))``

There is yet another way to create a matrix or vector of only *zero entries*:

.. code-block:: s

# 3x1 vector of zeros mat.or.vec(3, 1)

# 2x3 matrix of zeros mat.or.vec(2, 3)


Creating a matrix from several vectors

If you have a some equal-length vectors you might want to stack them together to create a matrix.

Stacking vectors side-by-side: ``cbind`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The ``cbind`` function will put each vector in a column and return the matrix. Let's say we wanted to use only a subset of the columns from a given data frame:

.. code-block:: s

bio <- read.csv('http://openmv.net/file/bioreactor-yields.csv') bio temperature duration speed baffles yield 1 82 260 4300 No 51 2 90 260 3700 Yes 30 .. ... ... .... ... .. 13 101 260 4400 No 42 14 92 260 4900 Yes 38

data <- cbind(bio$temperature, bio$ speed) data [,1] [,2] [1,] 82 4300 [2,] 90 3700 ..... ... .... [13,] 101 4400 [14,] 92 4900


Stacking top-to-bottom: ``rbind`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Use the ``rbind`` function to join existing matrices or vectors, one on top of another:

Create a matrix with 2 rows and 4 columns filled with random numbers, and another matrix containing 3 rows and 4 columns, filled with the number 2.6. Then to stack them:

.. code-block:: s

X <- matrix(rnorm(2*4), 2, 4) Y <- matrix(data=2.6, nrow=3, ncol=4)

rbind(X,Y)

[,1] [,2] [,3] [,4] [1,] 0.6096289 -0.3730975 -1.0376647 -0.01033531 [2,] 0.6305808 -0.6768284 0.4918309 -0.44430265 [3,] 2.6000000 2.6000000 2.6000000 2.60000000 [4,] 2.6000000 2.6000000 2.6000000 2.60000000 [5,] 2.6000000 2.6000000 2.6000000 2.60000000

Reshaping a vector into a matrix

.. TODO:KGD : link to PID book here later on

There are occasions when we need to create subgroups of data from a vector. This happens, for example, when constructing control charts:

  • we need to create the subgroups from *non-overlapping* segments in the vector
  • we need to do calculations on these subgroups

One way to construct the subgroups is to use 2 nested :ref:`for-loops <r-programming-loops-for-loop>` to extract the necessary data. This works, but can get messy. There is an alternative way: by rearranging your vector into a matrix and then doing your calculations on each column of the matrix.

Let the raw data be given by this random vector:

.. code-block:: s

N.raw = 5000 raw <- rnorm(N.raw, mean=220, sd=50) plot(raw)

Let's say our subgroup size is 9. Now 5000 is not neatly divisible by 9, but this doesn't matter. We create a matrix with 9 rows and as many columns as we can (roughly 5000/9).

.. code-block:: s

N.sub = 9 subgroups <- matrix(raw, N.sub, N.raw/N.sub) # Warning message: # In matrix(raw, N.sub, N.raw/N.sub) : # data length [5000] is not a sub-multiple or multiple of the number of rows [9]

Notice that R gives a warning, not an error message. You can see what size matrix it created: 9 rows and 555 columns, so there are 4995 elements in the array. That's usually good enough: we throw away the last 5 entries that don't form a complete subgroup.

.. code-block:: s

dim(subgroups) [1] 9 555

nrow(subgroups) [1] 9

The next step is to calculate the mean and standard deviation of each subgroup, i.e. for each column in the new matrix. For that we will use the ``apply(...)`` command, which applies any R function to a matrix along a particular direction.

  • ``apply(X, 1, sd)`` will apply the ``sd`` function to matrix ``X`` across the row direction (that's what the ``1`` is for: the first dimension)
  • ``apply(X, 2, sd)`` will apply the ``sd`` function to matrix ``X`` across the column direction (that's what the ``2`` is for: the second array dimension)

In our example the subgroups appear in each column. So calculate the standard deviation and mean for each column:

.. code-block:: s

subgroups.S <- apply(subgroups, 2, sd) subgroups.xbar <- apply(subgroups, 2, mean)

length(subgroups.S) [1] 555

# "subgroups.S" is not a matrix, so the "dim" function isn't valid dim(subgroups.S) NULL

    • Note**: there are built-in R functions, called ``colMeans(X)`` and ``rowMeans(X)`` which will perform the specific task of applying the ``mean`` function along the column or row direction. The ``apply`` function is more general though.


Matrix operations

Once you've got your matrices created it is time to work with them: addition and subtraction, multiplication, transposes, inverses and determinants can all be calculated in R.

Matrix addition and subtraction ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code-block:: s

X <- matrix(data=56, nrow=4, ncol=5) Y <- matrix(data=32, nrow=4, ncol=5) X+Y [,1] [,2] [,3] [,4] [,5] [1,] 88 88 88 88 88 [2,] 88 88 88 88 88 [3,] 88 88 88 88 88 [4,] 88 88 88 88 88

X-Y [,1] [,2] [,3] [,4] [,5] [1,] 24 24 24 24 24 [2,] 24 24 24 24 24 [3,] 24 24 24 24 24 [4,] 24 24 24 24 24

Matrix multiplication ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    • Element-by-element multiplication**

.. code-block:: s

a <- matrix(data=5, nrow=4, ncol=5) b <- matrix(data=6, nrow=4, ncol=5)

a*b [,1] [,2] [,3] [,4] [,5] [1,] 30 30 30 30 30 [2,] 30 30 30 30 30 [3,] 30 30 30 30 30 [4,] 30 30 30 30 30


    • The usual algebraic matrix multiplication** - which is probably what you are expecting when multiplying matrices:

.. code-block:: s

a <- matrix(data=5, nrow=4, ncol=5) b <- matrix(data=6, nrow=5, ncol=3) dim(a) [1] 4 5 dim(b) [1] 5 3

a %*% b [,1] [,2] [,3] [1,] 150 150 150 [2,] 150 150 150 [3,] 150 150 150 [4,] 150 150 150

Matrix transpose: ``t(X)`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code-block:: s

X <- matrix(data=c(56,12), nrow=4, ncol=5) X [,1] [,2] [,3] [,4] [,5] [1,] 56 56 56 56 56 [2,] 12 12 12 12 12 [3,] 56 56 56 56 56 [4,] 12 12 12 12 12

t(X) [,1] [,2] [,3] [,4] [1,] 56 12 56 12 [2,] 56 12 56 12 [3,] 56 12 56 12 [4,] 56 12 56 12 [5,] 56 12 56 12

Solving a system of equations :math:`Ax = b`: ``solve(A, b)`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Use the ``solve(a, b, ...)`` function in R. The ``solve()`` function solves the linear system of equations :math:`Ax = b` and returns :math:`x`.

.. Also see the discussion :ref:`in the QR decomposition section <r-vectors-matrices-qr-decomposition>` later on.

.. code-block:: s

A <- matrix(data = c(2,4,6, 5, 2, 7, 6, 2, 2), 3, 3) b <- matrix(data = c(6, 2, 2), nrow=3) x <- solve(A, b)

x [,1] [1,] 0 [2,] 0 [3,] 1

# Should be a vector of zeros: check <- A %*% x - b


Matrix inverse: ``solve(A)`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Recall that only square matrices, :math:`n \times n`, can be inverted. We use the ``solve(...)`` function in R, which solves the system :math:`Ax = b`, but if we replace ``b`` with successive columns from the identity matrix, we can successively build up a solution to :math:`AA^{-1} = I`.

Fortunately in R, all this work is done for you automatically. The ``solve()`` function creates an identity matrix if you don't supply the ``b`` argument. So to invert a matrix :math:`X`:

.. code-block:: s

n = 4 A <- matrix(rnorm(n*n), n, n) A.inv <- solve(A) [,1] [,2] [,3] [,4] [1,] 0.6564071 0.09397735 0.769666 -0.836499 [2,] 6.5235936 3.83526471 15.547712 -24.572593 [3,] 3.4061013 2.48309826 8.647189 -13.329604 [4,] 0.7537863 0.97065805 1.759342 -3.509522

# Should be the identity matrix check <- A %*% A.inv

# which it is, within machine precision: [,1] [,2] [,3] [,4] [1,] 1.000000e+00 0.000000e+00 -1.110223e-16 -8.881784e-16 [2,] -9.714451e-17 1.000000e+00 1.387779e-17 5.551115e-17 [3,] -5.551115e-17 0.000000e+00 1.000000e+00 0.000000e+00 [4,] -2.775558e-17 -4.440892e-16 5.551115e-17 1.000000e+00

The `generalized matrix inverse <http://en.wikipedia.org/wiki/Generalized_inverse>`_ can also be calculated, after loading the built-in ``MASS`` package, use the ``ginv(...)`` function.

Matrix determinant: ``det(X)`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Recall that the determinant is only computable for a square, :math:`n \times n` matrix.

.. code-block:: s

N = 4 X <- diag(x=3, N, N) det(X) [1] 81


Singular value decomposition: ``svd(X)`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The singular value decomposition (SVD) is useful for solving certain principal component analysis models, amongst other uses. The SVD is given by:

.. math::

\bf X = U D V'

where matrices :math:`\bf{U}` and :math:`\bf{V}` are orthogonal, and matrix :math:`\bf{D}` is a diagonal matrix, containing the singular values. For efficiency, the entries in :math:`\bf{D}` are returned as a vector.

Use the example below as a guide to using the ``svd`` function in R:

.. code-block:: s

X <- matrix(data=rnorm(10*5), 10, 5) decomp <- svd(X) d.vector <- decomp$d U <- decomp$u V <- decomp$v

Other matrix operations ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Other important matrix operations that can be performed efficiently in R are:

.. rubric:: The ``Matrix`` package

Enhanced matrix capability is provided by the ``Matrix`` package, which is not loaded by default. To start using it, load the library as usual: ``library(Matrix)``. Type ``help(package="Matrix")`` for more details. </rst>

← Programming in R: loops and flow control (previous step) Tutorial index Next step: Building a least squares model in R →
Retrieved from "https://learnche.org/wiki_4C3/index.php?title=Software_tutorial/Vectors_and_matrices&oldid=1597"