The first question concerns the computation in S-Plus of the square root of a symmetric nonnegative-definite square matrix.
1. Write an S-function, call it msqrt, with argument A which:
- checks that A is a square matrix and exits if not;
- checks that A is symmetric and exits if not;
- diagonalizes the matrix by computing the eigenvalues and the matrix of eigenvectors
- checks that all the eigenvalues are nonnegative and exits, if not;
- returns a symmetric matrix of the same size as A, with the same eigenvectors, the eigenvalue corresponding to a given eigenvector being the square root of the corresponding eigenvalue of A
The matrix returned by such a function msqrt is called the square root of the matrix A and it will be denoted by A1/2.
The second question concerns the generation of normal random vectors in S-Plus. In other words, we write an S-Plus function to play the role of the function rnorm in the case of multidimensional random vectors. Such a function does exist in the S-Plus distribution. It is called mvrnorm. The goal of this second question is to understand how such a generation method works
2. Write an S-function, call it vnorm, with arguments Mu, Sigma and N which:
- checks that Mu is a vector, exits if not, and otherwise reads its dimension, say L;
- checks that Sigma is an LxL symmetric matrix with nonnegative eigenvalues and exits, if not;
- creates a numeric array with N rows and L columns and fills it with independent random numbers with the standard normal distribution N(0, 1);
- treats each row of this array as a vector, and multiplies it by the square root of the matrix Sigma (as computed in question I.1 above) and adds the vector Mu to the result;
- returns the random array modified in this way. The array produced by the function vnorm is a sample of size N of L-dimensional random vectors (arranged as rows of the matrix outputted by vnorm) with the normal distribution with mean Mu and variance/covariance matrix Sigma. Indeed, this function implements the following simple fact reviewed during the lectures:
If X is an L-dimensional normal vector with mean 0 and covariance matrix given by the L×L identity matrix (i.e. if all the L entries of X are independent N(0, 1) random variables), then:
is an L-dimensional normal vector with mean µ and variance/covariance matrix Σ.