Package 'fastmatrix'

Description

Multiplication of 3-dimensional arrays was first introduced by Bates and Watts (1980). More extensions and technical details can be found in Wei (1998).

Usage

array.mult(a, b, x)

Arguments

Details

Let $\bold{X} = (x_{tij})$ be a 3-dimensional $n\times p\times q$ where indices $t, i$ and $j$ indicate face, row and column, respectively. The product $\bold{Y} = \bold{AXB}$ is an $n\times r\times s$ array, with $\bold{A}$ and $\bold{B}$ are $r\times p$ and $q\times s$ matrices respectively. The elements of $\bold{Y}$ are defined as:

$y_{tkl} = \sum\limits_{i=1}^p\sum\limits_{j=1}^q a_{ki}x_{tij}b_{jl}$

Value

array.mult returns a 3-dimensional array of dimension $n\times r\times s$.

References

Bates, D.M., Watts, D.G. (1980). Relative curvature measures of nonlinearity. Journal of the Royal Statistical Society, Series B 42, 1-25.

Wei, B.C. (1998). Exponential Family Nonlinear Models. Springer, New York.

See Also

array, matrix, bracket.prod.

Examples

x <- array(0, dim = c(2,3,3)) # 2 x 3 x 3 array
x[,,1] <- c(1,2,2,4,3,6)
x[,,2] <- c(2,4,4,8,6,12)
x[,,3] <- c(3,6,6,12,9,18)

a <- matrix(1, nrow = 2, ncol = 3)
b <- matrix(1, nrow = 3, ncol = 2)

y <- array.mult(a, b, x) # a 2 x 2 x 2 array
y

Description

Force a square matrix x to be symmetric

Usage

asSymmetric(x, lower = TRUE)

Arguments

Value

a square symmetric matrix.

Examples

a <- matrix(1:16, ncol = 4)
isSymmetric(a) # FALSE
a <- asSymmetric(a) # copy lower triangle into upper triangle

Description

Computes the Bezier curve based on $n+1$ control points using the De Casteljau's method.

Usage

bezier(x, y, ngrid = 200)

Arguments

Details

Given $\bold{p}_0,\bold{p}_1,\dots,\bold{p}_n$ control points the Bezier curve is given by $B(t)$ defined as

$B(t) = \left({\begin{array}{c} x(t) \\ y(t) \end{array}}\right) = \sum\limits_{k=0}^n {n \choose k} t^k (1 - t)^k\bold{p}_k$

where $t\in[0,1]$. To evaluate the Bezier curve the De Casteljau's method is used.

Value

A list containing xgrid and ygrid elements used to plot the Bezier curve.

Examples

# a tiny example
x <- c(1.0, 0.25, 1.25, 2.5, 4.00, 5.0)
y <- c(0.5, 2.00, 3.75, 4.0, 3.25, 1.0)
plot(x, y, type = "o")
z <- bezier(x, y, ngrid = 50)
lines(z$xgrid, z$ygrid, lwd = 2, lty = 2, col = "red")

# other simple example
x <- c(4,6,4,5,6,7)
y <- 1:6
plot(x, y, type = "o")
z <- bezier(x, y, ngrid = 50)
lines(z$xgrid, z$ygrid, lwd = 2, lty = 2, col = "red")

Description

Bracket product of a matrix and a 3-dimensional array.

Usage

bracket.prod(a, x)

Arguments

Details

Let $\bold{X} = (x_{tij})$ be a 3-dimensional $n\times p\times q$ array and $\bold{A}$ an $m\times n$ matrix, then $\bold{Y} = [\bold{A}][\bold{X}]$ is called the bracket product of $\bold{A}$ and $\bold{X}$, that is an $m \times p\times q$ with elements

$y_{tij} = \sum\limits_{k=1}^n a_{tk}x_{kij}$

Value

bracket.prod returns a 3-dimensional array of dimension $m\times p\times q$.

References

Wei, B.C. (1998). Exponential Family Nonlinear Models. Springer, New York.

See Also

array, matrix, array.mult.

Examples

x <- array(0, dim = c(2,3,3)) # 2 x 3 x 3 array
x[,,1] <- c(1,2,2,4,3,6)
x[,,2] <- c(2,4,4,8,6,12)
x[,,3] <- c(3,6,6,12,9,18)

a <- matrix(1, nrow = 3, ncol = 2)

y <- bracket.prod(a, x) # a 3 x 3 x 3 array
y

Description

Conjugate gradients (CG) method is an iterative algorithm for solving linear systems with positive definite coefficient matrices.

Usage

cg(a, b, maxiter = 200, tol = 1e-7)

Arguments

Value

a vector with the approximate solution, the iterations performed are returned as the attribute 'iterations'.

Warning

The underlying C code does not check for symmetry nor positive definitiveness.

References

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

Hestenes, M.R., Stiefel, E. (1952). Methods of conjugate gradients for solving linear equations. Journal of Research of the National Bureau of Standards 49, 409-436.

See Also

jacobi, seidel, solve

Examples

a <- matrix(c(4,3,0,3,4,-1,0,-1,4), ncol = 3)
b <- c(24,30,-24)
z <- cg(a, b)
z # converged in 3 iterations

Description

function cholupdate, where R = chol(A) is the original Cholesky factorization of $\bold{A}$, returns the upper triangular Cholesky factor of $\bold{A} + \bold{xx}^T$, with $\bold{x}$ a column vector of appropriate dimension.

Usage

cholupdate(r, x)

Arguments

References

Golub, G.H., Van Loan, C.F. (2013). Matrix Computations, 4th Edition. John Hopkins University Press.

See Also

chol

Examples

a <- matrix(c(1,1,1,1,2,3,1,3,6), ncol = 3)
r <- chol(a)
x <- c(0,0,1)
b <- a + outer(x,x)
r1 <- cholupdate(r, x)
r1
all(r1 == chol(b)) # TRUE

Description

Forms a symmetric circulant matrix using a backwards shift of its first column

Usage

circulant(x)

Arguments

Value

A symmetric circulant matrix.

Examples

x <- c(2,3,5,7,11,13)
circulant(x)

Description

This function provides the minimum information required to create the commutation matrix.

The commutation matrix is a square matrix of order $mn$ that, for an $m\times n$ matrix $\bold{A}$, transform vec$(\bold{A}$) to vec$(\bold{A}^T)$.

Usage

comm.info(m = 1, n = m, condensed = TRUE)

Arguments

Details

This function returns a list containing two vectors that represent an element of the commutation matrix and is accesed by the indexes in vectors row and col. This information is used by function comm.prod to do some operations involving the commutation matrix without forming it. This information also can be obtained using function commutation.

Value

A list containing the following elements:

References

Magnus, J.R., Neudecker, H. (1979). The commutation matrix: some properties and applications. The Annals of Statistics 7, 381-394.

See Also

commutation, comm.prod

Examples

z <- comm.info(m = 3, n = 2, condensed = FALSE)
z # where are the ones in commutation matrix of order '3,2'?

K32 <- commutation(m = 3, n = 2, matrix = TRUE)
K32 # only recommended if m and n are very small

Description

Given the row and column dimensions of a commutation matrix $\bold{K}$ of order $mn$ and a conformable matrix $\bold{x}$, performs one of the matrix-matrix operations:

• $\bold{Y} = \bold{KX}$, if side = "left" and transposed = FALSE, or

• $\bold{Y} = \bold{K}^T\bold{X}$, if side = "left" and transposed = TRUE, or

• $\bold{Y} = \bold{XK}$, if side = "right" and transposed = FALSE, or

• $\bold{Y} = \bold{XK}^T$, if side = "right" and transposed = TRUE.

The main aim of comm.prod is to do this matrix multiplication without forming the commutation matrix.

Usage

comm.prod(m = 1, n = m, x = NULL, transposed = FALSE, side = "left")

Arguments

Details

Underlying Fortran code only uses information provided by comm.info to performs the matrix multiplication. The commutation matrix is never created.

See Also

commutation

Examples

K42 <- commutation(m = 4, n = 2, matrix = TRUE)
x <- matrix(1:24, ncol = 3)
y <- K42 %*% x

z <- comm.prod(m = 4, n = 2, x) # K42 is not stored
all(z == y) # matrices y and z are equal!

Description

This function returns the commutation matrix of order $mn$ which transforms, for an $m\times n$ matrix $\bold{A}$, vec$(\bold{A})$ to vec$(\bold{A}^T)$.

Usage

commutation(m = 1, n = m, matrix = FALSE, condensed = FALSE)

Arguments

Details

This function is a wrapper function for the function comm.info. This function provides the minimum information required to create the commutation matrix. If option matrix = FALSE the commutation matrix is stored in two vectors containing the coordinate list of indexes for rows and columns. Option condensed = TRUE only returns vector of indexes for the rows of commutation matrix.

Warning: matrix = TRUE is not recommended, unless the order m and n be small. This matrix can require a huge amount of storage.

Value

Returns an $mn$ by $mn$ matrix (if requested).

References

Magnus, J.R., Neudecker, H. (1979). The commutation matrix: some properties and applications. The Annals of Statistics 7, 381-394.

Magnus, J.R., Neudecker, H. (2007). Matrix Differential Calculus with Applications in Statistics and Econometrics, 3rd Edition. Wiley, New York.

See Also

comm.info

Examples

z <- commutation(m = 100, condensed = TRUE)
object.size(z) # 40.6 Kb of storage

z <- commutation(m = 100, condensed = FALSE)
object.size(z) # 80.7 Kb of storage

K100 <- commutation(m = 100, matrix = TRUE) # time: < 2 secs
object.size(K100) # 400 Mb of storage, do not request this matrix!

# a small example
K32 <- commutation(m = 3, n = 2, matrix = TRUE)
a <- matrix(1:6, ncol = 2)
v <- K32 %*% vec(a)
all(vec(t(a)) == as.vector(v)) # vectors are equal!

Description

This function is a constructor for the corAR1 correlation matrix representing an autocorrelation structure of order 1.

Usage

corAR1(rho, p = 2)

Arguments

Value

Returns a $p$ by $p$ matrix.

Examples

R <- corAR1(rho = 0.8, p = 5)

Description

This function is a constructor for the corCS correlation matrix representing a compound symmetry structure corresponding to uniform correlation.

Usage

corCS(rho, p = 2)

Arguments

Value

Returns a $p$ by $p$ matrix.

Examples

R <- corCS(rho = 0.8, p = 5)

Description

Returns a list containing the mean and covariance matrix of the data.

Usage

cov.MSSD(x)

Arguments

Details

This procedure uses the Holmes-Mergen method using the difference between each successive pairs of observations also known as Mean Square Successive Method (MSSD) to estimate the covariance matrix, which is given by

$\bold{S}_{HD} = \frac{1}{2(n-1)} \sum\limits_{i=2}^n (\bold{x}_i - \bold{x}_{i-1})(\bold{x}_i - \bold{x}_{i-1})^T.$

Value

A list containing the following named components:

References

Holmes, D.S., Mergen, A.E. (1993). Improving the performance of the $T^2$ control chart. Quality Engineering 5, 619-625.

See Also

cov and var.

Examples

x <- cbind(1:10, c(1:3, 8:5, 8:10))
z0 <- cov(x)
z0
z1 <- cov.MSSD(x)
z1

Description

Returns a list containing estimates of the weighted mean and covariance matrix of the data.

Usage

cov.weighted(x, weights = rep(1, nrow(x)))

Arguments

Details

The covariance matrix is divided by the number of observations, which arise for instance, when we use the class of elliptical contoured distributions. Thus,

$W_n = \sum\limits_{i=1}^n w_i, \qquad \overline{\bold{x}}_n = \frac{1}{W_n} \sum\limits_{i=1}^n w_i\bold{x}_i \qquad \bold{S}_n = \frac{1}{n} \sum\limits_{i=1}^n w_i (\bold{x}_i - \overline{\bold{x}}_n)(\bold{x}_i - \overline{\bold{x}}_n)^T.$

This differs from the behaviour of function cov.wt.

Value

A list containing the following named components:

References

Clarke, M.R.B. (1971). Algorithm AS 41: Updating the sample mean and dispersion matrix. Applied Statistics 20, 206-209.

See Also

cov.wt, cov and var.

Examples

x <- cbind(1:10, c(1:3, 8:5, 8:10))
z0 <- cov.weighted(x) # all weights are 1
D2 <- Mahalanobis(x, center = z0$mean, cov = z0$cov)
p <- ncol(x)
wts <- (p + 1) / (1 + D2) # nice weights!
z1 <- cov.weighted(x, weights = wts)
z1

Description

Given the order of two duplication matrices and a conformable matrix $\bold{X}$, this function performs the operation: $\bold{Y} = \bold{D}_n^T\bold{X}\bold{D}_k$, where $\bold{D}_n$ and $\bold{D}_k$ are duplication matrices of order $n$ and $k$, respectively.

Usage

dupl.cross(n = 1, k = n, x = NULL)

Arguments

Details

This function calls dupl.prod to performs the matrix multiplications required but without forming any duplication matrices.

See Also

dupl.prod

Examples

D2 <- duplication(n = 2, matrix = TRUE)
D3 <- duplication(n = 3, matrix = TRUE)
x <- matrix(1, nrow = 9, ncol = 4)
y <- t(D3) %*% x %*% D2

z <- dupl.cross(n = 3, k = 2, x) # D2 and D3 are not stored
all(z == y) # matrices y and z are equal!

x <- matrix(1, nrow = 9, ncol = 9)
z <- dupl.cross(n = 3, x = x) # same matrix is used to pre- and post-multiplying x
z # print result

Description

This function provides the minimum information required to create the duplication matrix.

Usage

dupl.info(n = 1, condensed = TRUE)

Arguments

Details

This function returns a list containing two vectors that represent an element of the duplication matrix and is accesed by the indexes in vectors row and col. This information is used by function dupl.prod to do some operations involving the duplication matrix without forming it. This information also can be obtained using function duplication

Value

A list containing the following elements:

See Also

duplication, dupl.prod

Examples

z <- dupl.info(n = 3, condensed = FALSE)
z # where are the ones in duplication of order 3?

D3 <- duplication(n = 3, matrix = TRUE)
D3 # only recommended if n is very small

Description

Given the order of a duplication and a conformable matrix $\bold{X}$, performs one of the matrix-matrix operations:

• $\bold{Y} = \bold{DX}$, if side = "left" and transposed = FALSE, or

• $\bold{Y} = \bold{D}^T\bold{X}$, if side = "left" and transposed = TRUE, or

• $\bold{Y} = \bold{XD}$, if side = "right" and transposed = FALSE, or

• $\bold{Y} = \bold{XD}^T$, if side = "right" and transposed = TRUE,

where $\bold{D}$ is the duplication matrix of order $n$. The main aim of dupl.prod is to do this matrix multiplication without forming the duplication matrix.

Usage

dupl.prod(n = 1, x, transposed = FALSE, side = "left")

Arguments

Details

Underlying C code only uses information provided by dupl.info to performs the matrix multiplication. The duplication matrix is never created.

See Also

duplication

Examples

D4 <- duplication(n = 4, matrix = TRUE)
x <- matrix(1, nrow = 16, ncol = 2)
y <- crossprod(D4, x)

z <- dupl.prod(n = 4, x, transposed = TRUE) # D4 is not stored
all(z == y) # matrices y and z are equal!

Description

This function returns the duplication matrix of order $n$ which transforms, for a symmetric matrix $\bold{A}$, vech$(\bold{A})$ into vec$(\bold{A})$.

Usage

duplication(n = 1, matrix = FALSE, condensed = FALSE)

Arguments

Details

This function is a wrapper function for the function dupl.info. This function provides the minimum information required to create the duplication matrix. If option matrix = FALSE the duplication matrix is stored in two vectors containing the coordinate list of indexes for rows and columns. Option condensed = TRUE only returns vector of indexes for the columns of duplication matrix.

Warning: matrix = TRUE is not recommended, unless the order n be small. This matrix can require a huge amount of storage.

Value

Returns an $n^2$ by $n(n + 1)/2$ matrix (if requested).

References

Magnus, J.R., Neudecker, H. (1980). The elimination matrix, some lemmas and applications. SIAM Journal on Algebraic Discrete Methods 1, 422-449.

Magnus, J.R., Neudecker, H. (2007). Matrix Differential Calculus with Applications in Statistics and Econometrics, 3rd Edition. Wiley, New York.

See Also

dupl.info

Examples

z <- duplication(n = 100, condensed = TRUE)
object.size(z) # 40.5 Kb of storage

z <- duplication(n = 100, condensed = FALSE)
object.size(z) # 80.6 Kb of storage

D100 <- duplication(n = 100, matrix = TRUE)
object.size(D100) # 202 Mb of storage, do not request this matrix!

# a small example
D3 <- duplication(n = 3, matrix = TRUE)
a <- matrix(c( 1, 2, 3,
               2, 3, 4,
               3, 4, 5), nrow = 3)
upper <- vech(a)
v <- D3 %*% upper
all(vec(a) == as.vector(v)) # vectors are equal!

Description

Equilibrate a rectangular or symmetric matrix using 2-norm.

Usage

equilibrate(x, scale = TRUE)

Arguments

Value

For scale = TRUE, the equilibrated matrix. The scalings and an approximation of the condition number, are returned as attributes "scales" and "condition". If x is a rectangular matrix, only the columns are equilibrated.

Examples

x <- matrix(c(1, 1, 1,
              1, 2, 1,
              1, 3, 1,
              1, 1,-1,
              1, 2,-1,
              1, 3,-1), ncol = 3, byrow = TRUE)
z <- equilibrate(x)
apply(z, 2, function(x) sum(x^2)) # all 1

xx <- crossprod(x)
equilibrate(xx)

Description

This function returns the Frank matrix of order $n$.

Usage

frank(n = 1)

Arguments

Details

A Frank matrix of order $n$ is a square matrix $\bold{F}_n = (f_{ij})$ defined as

$f_{ij} = \left\{ {\begin{array}{ll} n - j + 1, & i \le j, \\ n - j, & i = j + 1, \\ 0, & i \ge j + 2 \end{array}} \right.$

Value

Returns an $n$ by $n$ matrix.

References

Frank, W.L. (1958). Computing eigenvalues of complex matrices by determinant evaluation and by methods of Danilewski and Wielandt. Journal of the Society for Industrial and Applied Mathematics 6, 378-392.

Examples

F5 <- frank(n = 5)
det(F5) # equals 1

Description

It calculates the geometric mean using a Fused-Multiply-and-Add (FMA) compensated scheme for accurate computation of floating-point product.

Usage

geomean(x)

Arguments

Details

If x contains any non-positive values, geomean returns NA and a warning message is displayed.

The geometric mean is a measure of central tendency, which is defined as

$G = \sqrt[n]{x_1 x_2 \ldots x_n} = \Big(\prod_{i=1}^n x_i\Big)^{1/n}.$

This procedure calculates the product required in the geometric mean safely using a compensated scheme as proposed by Graillat (2009).

Value

The geometric mean of the sample, a non-negative number.

References

Graillat, S. (2009). Accurate floating-point product and exponentiation. IEEE Transactions on Computers 58, 994-1000.

Oguita, T., Rump, S.M., Oishi, S. (2005). Accurate sum and dot product. SIAM Journal on Scientific Computing 26, 1955-1988.

See Also

mean, median.

Examples

set.seed(149)
x <- rlnorm(1000)
mean(x)    # 1.68169
median(x)  # 0.99663
geomean(x) # 1.01688

Description

This function returns the Hadamard or element-wise product of two matrices x and y, that have the same dimensions.

Usage

hadamard(x, y = x)

Arguments

Value

A matrix with the same dimension of x (and y) which corresponds to the element-by-element product of the two matrices.

References

Styan, G.P.H. (1973). Hadamard products and multivariate statistical analysis, Linear Algebra and Its Applications 6, 217-240.

Examples

x <- matrix(rep(1:10, times = 5), ncol = 5)
y <- matrix(rep(1:5, each = 10), ncol = 5)
z <- hadamard(x, y)
z

Description

Performs large-sample methods for testing equality of $p \ge 2$ correlated variables.

Usage

harris.test(x, test = "Wald")

Arguments

Value

A list of class 'harris.test' with the following elements:

References

Harris, P. (1985). Testing the variance homogeneity of correlated variables. Biometrika 72, 103-107.

Examples

x <- iris[,1:4]
z <- harris.test(x, test = "robust")
z

Description

This function returns the Helmert matrix of order $n$.

Usage

helmert(n = 1)

Arguments

Details

A Helmert matrix of order $n$ is a square matrix defined as

$\bold{H}_n = \left[ {\begin{array}{ccccc} 1/\sqrt{n} & 1/\sqrt{n} & 1/\sqrt{n} & \dots & 1/\sqrt{n} \\ 1/\sqrt{2} & -1/\sqrt{2} & 0 & \dots & 0 \\ 1/\sqrt{6} & 1/\sqrt{6} & -2/\sqrt{6} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{1}{\sqrt{n(n-1)}} & \frac{1}{\sqrt{n(n-1)}} & \frac{1}{\sqrt{n(n-1)}} & \dots & -\frac{(n-1)}{\sqrt{n(n-1)}} \end{array}} \right].$

Helmert matrix is orthogonal and is frequently used in the analysis of variance (ANOVA).

Value

Returns an $n$ by $n$ matrix.

References

Lancaster, H.O. (1965). The Helmert matrices. The American Mathematical Monthly 72, 4-12.

Gentle, J.E. (2007). Matrix Algebra: Theory, Computations, and Applications in Statistics. Springer, New York.

Examples

n <- 1000
set.seed(149)
x <- rnorm(n)

H <- helmert(n)
object.size(H) # 7.63 Mb of storage
K <- H[2:n,]
z <- c(K %*% x)
sum(z^2) # 933.1736

# same that
(n - 1) * var(x)

Description

Returns TRUE if the given matrix is lower or upper triangular matrix.

Usage

is.lower.tri(x, diag = FALSE)
is.upper.tri(x, diag = FALSE)

Arguments

Value

Check if a matrix is lower or upper triangular. You can also include diagonal to the check.

See Also

lower.tri, upper.tri

Examples

x <- matrix(rnorm(10 * 3), ncol = 3)
  R <- chol(crossprod(x))

  is.lower.tri(R)
  is.upper.tri(R)

Description

Jacobi method is an iterative algorithm for solving a system of linear equations.

Usage

jacobi(a, b, start, maxiter = 200, tol = 1e-7)

Arguments

Details

Let $\bold{D}$, $\bold{L}$, and $\bold{U}$ denote the diagonal, lower triangular and upper triangular parts of a matrix $\bold{A}$. Jacobi's method solve the equation $\bold{Ax} = \bold{b}$, iteratively by rewriting $\bold{Dx} + (\bold{L} + \bold{U})\bold{x} = \bold{b}$. Assuming that $\bold{D}$ is nonsingular leads to the iteration formula

$\bold{x}^{(k+1)} = -\bold{D}^{-1}(\bold{L} + \bold{U})\bold{x}^{(k)} + \bold{D}^{-1}\bold{b}$

Value

a vector with the approximate solution, the iterations performed are returned as the attribute 'iterations'.

References

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

seidel

Examples

a <- matrix(c(5,-3,2,-2,9,-1,3,1,-7), ncol = 3)
b <- c(-1,2,3)
start <- c(1,1,1)
z <- jacobi(a, b, start)
z # converged in 15 iterations

Description

Performs an omnibus test for univariate normality.

Usage

JarqueBera.test(x, test = "DH")

Arguments

Value

A list of class 'JarqueBera.test' with the following elements:

References

Doornik, J.A., Hansen, H. (2008). An omnibus test for univariate and multivariate normality. Oxford Bulletin of Economics and Statistics 70, 927-939.

Gel, Y.R., Gastwirth, J.L. (2008). A robust modification of the Jarque-Bera test of normality. Economics Letters 99, 30-32.

Jarque, C.M., Bera, A.K. (1980). Efficient tests for normality, homoscedasticity and serial independence of regression residuals. Economics Letters 6, 255-259.

Urzua, C.M. (1996). On the correct use of omnibus tests for normality. Economics Letters 53, 247-251.

Examples

set.seed(149)
x <- rnorm(100)
z <- JarqueBera.test(x, test = "DH")
z

set.seed(173)
x <- runif(100)
z <- JarqueBera.test(x, test = "DH")
z

Description

Computes the kronecker product of two matrices, x and y.

Usage

kronecker.prod(x, y = x)

Arguments

Details

Let $\bold{X}$ be an $m\times n$ and $\bold{Y}$ a $p\times q$ matrix. The $mp\times nq$ matrix defined by

$\left[{\begin{array}{ccc} x_{11}\bold{Y} & \dots & x_{1n}\bold{Y} \\ \vdots & & \vdots \\ x_{m1}\bold{Y} & \dots & x_{mn}\bold{Y} \end{array}}\right],$

is called the Kronecker product of $\bold{X}$ and $\bold{Y}$.

Value

An array with dimensions dim(x) * dim(y).

References

Magnus, J.R., Neudecker, H. (2007). Matrix Differential Calculus with Applications in Statistics and Econometrics, 3rd Edition. Wiley, New York.

See Also

kronecker function from base package is based on outer. Our C version is slightly faster.

Examples

# block diagonal matrix:
a <- diag(1:3)
b <- matrix(1:4, ncol = 2)
kronecker.prod(a, b)

# examples with vectors
ones <- rep(1, 4)
y <- 1:3
kronecker.prod(ones, y) # 12-dimensional vector
kronecker.prod(ones, t(y)) # 3 x 3 matrix

Description

Given $\bold{A}$ an $n$ by $n$ real matrix and an $n$-vector $\bold{b}$, this function constructs the Krylov matrix $\bold{K}$, where

$\bold{K} = [\bold{b},\bold{Ab},\dots,\bold{A}^{m-1}\bold{b}].$

Usage

krylov(a, b, m = ncol(a))

Arguments

Value

Returns an $n$ by $m$ matrix.

Examples

a <- matrix(c(1, 3, 2, -5, 1, 7, 1, 5, -4), ncol = 3, byrow = TRUE)
b <- c(1, 1, 1)
k <- krylov(a, b, m = 4)
k

Description

Functions to compute measures of multivariate skewness $(b_{1p})$ and kurtosis $(b_{2p})$ proposed by Mardia (1970),

$b_{1p} = \frac{1}{n^2}\sum\limits_{i=1}^n\sum\limits_{j=1}^n ((\bold{x}_i - \overline{\bold{x}})^T\bold{S}^{-1}(\bold{x}_j - \overline{\bold{x}}))^3,$

and

$b_{2p} = \frac{1}{n}\sum\limits_{i=1}^n ((\bold{x}_i - \overline{\bold{x}})^T \bold{S}^{-1}(\bold{x}_j - \overline{\bold{x}}))^2.$

Usage

kurtosis(x)

skewness(x)

Arguments

References

Mardia, K.V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika 57, 519-530.

Mardia, K.V., Zemroch, P.J. (1975). Algorithm AS 84: Measures of multivariate skewness and kurtosis. Applied Statistics 24, 262-265.

Examples

setosa <- iris[1:50,1:4]
kurtosis(setosa)
skewness(setosa)

Description

Compute the LDL decomposition of a real symmetric matrix.

Usage

ldl(x)

Arguments

Value

The factorization has the form $\bold{X} = \bold{LDL}^T$, where $\bold{D}$ is a diagonal matrix and $\bold{L}$ is unitary lower triangular.

The LDL decomposition of $\bold{x}$ is returned as a list with components:

References

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

chol

Examples

a <- matrix(c(2,-1,0,-1,2,-1,0,-1,1), ncol = 3)
z <- ldl(a)
z # information of LDL factorization

# computing det(a)
prod(z$d) # product of diagonal elements of D

# a non-positive-definite matrix
m <- matrix(c(5,-5,-5,3), ncol = 2)
try(chol(m)) # fails
ldl(m)

Description

lu computes the LU factorization of a matrix.

Usage

lu(x)
## Default S3 method:
lu(x)

## S3 method for class 'lu'
solve(a, b, ...)

is.lu(x)

Arguments

Details

The LU factorization plays an important role in many numerical procedures. In particular it is the basic method to solve the equation $\bold{Ax} = \bold{b}$ for given matrix $\bold{A}$, and vector $\bold{b}$.

solve.lu is the method for solve for lu objects.

is.lu returns TRUE if x is a list and inherits from "lu".

Unsuccessful results from the underlying LAPACK code will result in an error giving a positive error code: these can only be interpreted by detailed study of the Fortran code.

Value

The LU factorization of the matrix as computed by LAPACK. The components in the returned value correspond directly to the values returned by DGETRF.

Note

To compute the determinant of a matrix (do you really need it?), the LU factorization is much more efficient than using eigenvalues (eigen). See det.

LAPACK uses column pivoting and does not attempt to detect rank-deficient matrices.

References

Anderson. E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. Sorensen, D. (1999). LAPACK Users' Guide, 3rd Edition. SIAM.

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

extractL, extractU, constructX for reconstruction of the matrices, lu2inv

Examples

a <- matrix(c(3,2,6,17,4,18,10,-2,-12), ncol = 3)
z <- lu(a)
z # information of LU factorization

# computing det(a)
prod(diag(z$lu)) # product of diagonal elements of U

# solve linear equations
b <- matrix(1:6, ncol = 2)
solve(z, b)

Description

Returns the original matrix from which the object was constructed or the components of the factorization.

Usage

constructX(x)
extractL(x)
extractU(x)

Arguments

Value

constructX returns $\bold{X}$, the original matrix from which the lu object was constructed (because of the pivoting the $\bold{X}$ matrix is not exactly the product between $\bold{L}$ and $\bold{U}$).

extractL returns $\bold{L}$. This may be pivoted.

extractU returns $\bold{U}$.

See Also

lu.

Examples

a <- matrix(c(10,-3,5,-7,2,-1,0,6,5), ncol = 3)
z <- lu(a)
L <- extractL(z)
L
U <- extractU(z)
U
X <- constructX(z)
all(a == X)

Description

Invert a square matrix from its LU factorization.

Usage

lu2inv(x)

Arguments

Value

The inverse of the matrix whose LU factorization was given.

Unsuccessful results from the underlying LAPACK code will result in an error giving a positive error code: these can only be interpreted by detailed study of the Fortran code.

Source

This is an interface to the LAPACK routine DGETRI. LAPACK is from https://netlib.org/lapack/ and its guide is listed in the references.

References

Anderson. E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. Sorensen, D. (1999). LAPACK Users' Guide, 3rd Edition. SIAM.

Golub, G.H., Van Loan, C.F. (1996). Matrix Computations, 3rd Edition. John Hopkins University Press.

See Also

lu, solve.

Examples

a <- matrix(c(3,2,6,17,4,18,10,-2,-12), ncol = 3)
z <- lu(a)
a %*% lu2inv(z)

Description

Returns the squared Mahalanobis distance of all rows in $\bold{x}$ and the vector $\bold{\mu}$ = center with respect to $\bold{\Sigma}$ = cov. This is (for vector $\bold{x}$) defined as

$D^2 = (\bold{x} - \bold{\mu})^T \bold{\Sigma}^{-1} (\bold{x} - \bold{\mu})$

Usage

Mahalanobis(x, center, cov, inverted = FALSE)

Arguments

Details

Unlike function mahalanobis, the covariance matrix is factorized using the Cholesky decomposition, which allows to assess if cov is positive definite. Unsuccessful results from the underlying LAPACK code will result in an error message.

See Also

cov, mahalanobis

Examples

x <- cbind(1:6, 1:3)
xbar <- colMeans(x)
S <- matrix(c(1,4,4,1), ncol = 2) # is negative definite
D2 <- mahalanobis(x, center = xbar, S)
all(D2 >= 0) # several distances are negative

## next command produces the following error:
## Covariance matrix is possibly not positive-definite
## Not run: D2 <- Mahalanobis(x, center = xbar, S)

Description

Computes the inner product between two rectangular matrices calling BLAS.

Usage

matrix.inner(x, y = x)

Arguments

Value

a real value, indicating the inner product between two matrices.

Examples

x <- matrix(c(1, 1, 1,
              1, 2, 1,
              1, 3, 1,
              1, 1,-1,
              1, 2,-1,
              1, 3,-1), ncol = 3, byrow = TRUE)
y <- matrix(1, nrow = 6, ncol = 3)
matrix.inner(x, y)

# must be equal
matrix.norm(x, type = "Frobenius")^2
matrix.inner(x)

Description

Computes a matrix norm of x using LAPACK. The norm can be the one ("1") norm, the infinity ("inf") norm, the Frobenius norm, the maximum modulus ("maximum") among elements of a matrix, as determined by the value of type.

Usage

matrix.norm(x, type = "Frobenius")

Arguments

Details

As function norm in package base, method of matrix.norm calls the LAPACK function DLANGE.

Note that the 1-, Inf- and maximum norm is faster to calculate than the Frobenius one.

Value

The matrix norm, a non-negative number.

Examples

# a tiny example
x <- matrix(c(1, 1, 1,
              1, 2, 1,
              1, 3, 1,
              1, 1,-1,
              1, 2,-1,
              1, 3,-1), ncol = 3, byrow = TRUE)
matrix.norm(x, type = "Frobenius")
matrix.norm(x, type = "1")
matrix.norm(x, type = "Inf")

# an example not that small
n <- 1000
x <- .5 * diag(n) + 0.5 * matrix(1, nrow = n, ncol = n)
matrix.norm(x, type = "Frobenius")
matrix.norm(x, type = "1")
matrix.norm(x, type = "Inf")
matrix.norm(x, type = "maximum") # equal to 1

Description

Compute the Cholesky factorization of a real symmetric but not necessarily positive definite matrix.

Usage

mchol(x)

Arguments

Value

The lower triangular factor of modified Cholesky decomposition, i.e., the matrix $\bold{L}$ such that $\bold{X} + \bold{E} = \bold{LL}^T$, where $\bold{E}$ is a nonnegative diagonal matrix that is zero if $\bold{X}$ es sufficiently positive definite.

References

Gill, P.E., Murray, W., Wright, M.H. (1981). Practical Optimization. Academic Press, London.

Nocedal, J., Wright, S.J. (1999). Numerical Optimization. Springer, New York.

See Also

chol, ldl

Examples

# a non-positive-definite matrix
a <- matrix(c(4,2,1,2,6,3,1,3,-.004), ncol = 3)
try(chol(a)) # fails
z <- mchol(a)
z # triangular factor

# modified 'a' matrix
tcrossprod(z)

Description

It calculates the mediancenter (or geometric median) of multivariate data.

Usage

mediancenter(x)

Arguments

Details

The mediancenter for a sample of multivariate observations is computed using a steepest descend method combined with bisection. The mediancenter invariant to rotations of axes and is useful as a multivariate generalization of the median of univariate sample.

Value

A list containing the following named components:

References

Gower, J.C. (1974). Algorithm AS 78: The mediancentre. Applied Statistics 23, 466-470.

See Also

cov.wt, median.

Examples

x <- cbind(1:10, c(1:3, 8:5, 8:10))
z <- mediancenter(x)$median # degenerate solution
xbar <- colMeans(x)
plot(x, xlab = "", ylab = "")
points(x = xbar[1], y = xbar[2], pch = 16, col = "red")
points(x = z[1], y = z[2], pch = 3, col = "blue", lwd = 2)

Description

Computes a $p$-norm of vector $\bold{x}$. The norm can be the one ($p = 1$) norm, Euclidean ($p = 2$) norm, the infinity ($p$ = Inf) norm. The underlying C or Fortran code is inspired on ideas of BLAS Level 1.

Usage

minkowski(x, p = 2)

Arguments

Details

Method of minkowski for $p$ = Inf calls idamax BLAS function. For other values, C or Fortran subroutines using unrolled cycles are called.

Value

The vector $p$-norm, a non-negative number.

Examples

# a tiny example
x <- rnorm(1000)
minkowski(x, p = 1)
minkowski(x, p = 1.5)
minkowski(x, p = 2)
minkowski(x, p = Inf)

x <- x / minkowski(x)
minkowski(x, p = 2) # equal to 1

Description

It calculates up to fourth central moments (or moments about the mean), and the skewness and kurtosis coefficients using an online algorithm.

Usage

moments(x)

Arguments

Details

The $k$-th central moment is defined as

$m_k = \frac{1}{n}\sum_{i=1}^n (x_i - \overline{x})^k.$

In particular, the second central moment is the variance of the sample. The sample skewness and kurtosis are defined, respectively, as

$b_1 = \frac{m_3}{m_2^{3/2}}, \qquad b_2 = \frac{m_4}{m_2^2}.$

Value

A list containing second, third and fourth central moments, and skewness and kurtosis coefficients.

References

Spicer, C.C. (1972). Algorithm AS 52: Calculation of power sums of deviations about the mean. Applied Statistics 21, 226-227.

See Also

var.

Examples

set.seed(149)
x <- rnorm(1000)
z <- moments(x)
z

Description

Returns an object of class "ols" that represents a linear model fit.

Usage

ols(formula, data, subset, na.action, method = "qr", tol = 1e-7, maxiter = 100,
  x = FALSE, y = FALSE, contrasts = NULL, ...)

Arguments

Value

ols returns an object of class "ols".

The function summary is used to obtain and print a summary of the results. The generic accessor functions coefficients, fitted.values and residuals extract various useful features of the value returned by ols.

An object of class "ols" is a list containing at least the following components:

See Also

ols.fit, lm, lsfit

Examples

# tiny example of regression
y <- c(1, 3, 3, 2, 2, 1)
x <- matrix(c(1, 1,
              2, 1,
              3, 1,
              1,-1,
              2,-1,
              3,-1), ncol = 2, byrow = TRUE)
f0 <- ols(y ~ x) # intercept is included by default
f0 # printing results (QR method was used)

f1 <- ols(y ~ x, method = "svd") # using SVD method instead
f1

Description

This function is a switcher among various numerical fitting functions (ols.fit.cg, ols.fit.chol, ols.fit.qr, ols.fit.svd and ols.fit.sweep). The argument method does the switching: "qr" for ols.fit.qr, etc. This should usually not be used directly unless by experienced users.

Usage

ols.fit(x, y, method = "qr", tol = 1e-7, maxiter = 100)

Arguments

Value

a list with components:

See Also

ols.fit.cg, ols.fit.chol, ols.fit.qr, ols.fit.svd, ols.fit.sweep.

Examples

set.seed(151)
n <- 100
p <- 2
x <- matrix(rnorm(n * p), n, p) # no intercept!
y <- rnorm(n)
fm <- ols.fit(x = x, y = y, method = "chol")
fm

Description

Fits a linear model, returning the bare minimum computations.

Usage

ols.fit.cg(x, y, tol = 1e-7, maxiter = 100)
ols.fit.chol(x, y)
ols.fit.qr(x, y)
ols.fit.svd(x, y)
ols.fit.sweep(x, y)

Arguments

Value

The bare bones of an ols object: the coefficients, residuals, fitted values, and some information used by summary.ols.

See Also

ols, ols.fit, lm

Examples

set.seed(151)
n <- 100
p <- 2
x <- matrix(rnorm(n * p), n, p) # no intercept!
y <- rnorm(n)
z <- ols.fit.chol(x, y)
z

Description

The power method seeks to determine the eigenvalue of maximum modulus, and a corresponding eigenvector.

Usage

power.method(x, only.value = FALSE, maxiter = 100, tol = 1e-8)

Arguments

Value

When only.value is not true, as by default, the result is a list with components "value" and "vector". Otherwise only the dominan eigenvalue is returned. The performed number of iterations to reach convergence is returned as attribute "iterations".

See Also

eigen for eigenvalues and eigenvectors computation.

Examples

n <- 1000
x <- .5 * diag(n) + 0.5 * matrix(1, nrow = n, ncol = n)

# dominant eigenvalue must be (n + 1) / 2
z <- power.method(x, only.value = TRUE)

Description

Random vector generation uniformly in the unitary ball.

Usage

rball(n = 1, p = 2)

Arguments

Details

The function rball is an interface to C routines, which make calls to subroutines from BLAS.

Value

If n = 1 a p-dimensional vector, otherwise a matrix of n rows of random vectors.

References

Hormann, W., Leydold, J., Derflinger, G. (2004). Automatic Nonuniform Random Variate Generation. Springer, New York.

See Also

runif

Examples

# generate the sample
z <- rball(n = 500)

# scatterplot of a random sample of 500 points uniformly distributed
# in the unitary ball
par(pty = "s")
plot(z, xlab = "x", ylab = "y")
title("500 points in the ball", font.main = 1)

Description

Fit a linear model by ridge regression, returning an object of class "ridge".

Usage

ridge(formula, data, subset, lambda = 1.0, method = "GCV", ngrid = 200, tol = 1e-07,
  maxiter = 50, na.action, x = FALSE, y = FALSE, contrasts = NULL, ...)

Arguments

Details

ridge function fits in linear ridge regression without scaling or centering the regressors and the response. In addition, If an intercept is present in the model, its coefficient is penalized.)

Value

A list with the following components:

References

Golub, G.H., Heath, M., Wahba, G. (1979). Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 21, 215-223.

Hoerl, A.E., Kennard, R.W., Baldwin, K.F. (1975). Ridge regression: Some simulations. Communication in Statistics 4, 105-123.

Hoerl, A.E., Kennard, R.W. (1970). Ridge regression: Biased estimation of nonorthogonal problems. Technometrics 12, 55-67.

Lawless, J.F., Wang, P. (1976). A simulation study of ridge and other regression estimators. Communications in Statistics 5, 307-323.

Lee, T.S (1987). Algorithm AS 223: Optimum ridge parameter selection. Applied Statistics 36, 112-118.

See Also

lm, ols

Examples

z <- ridge(GNP.deflator ~ ., data = longley, lambda = 4, method = "grid")
z # ridge regression on a grid over seq(0, 4, length = 200)

z <- ridge(GNP.deflator ~ ., data = longley)
z # ridge parameter selected using GCV (default)

Description

Random number generation from the multivariate normal (Gaussian) distribution.

Usage

rmnorm(n = 1, mean = rep(0, nrow(Sigma)), Sigma = diag(length(mean)))

Arguments

Details

The function rmnorm is an interface to C routines, which make calls to subroutines from LAPACK. The matrix decomposition is internally done using the Cholesky decomposition. If Sigma is not non-negative definite then there will be a warning message.

Value

If $n = 1$ a vector of the same length as mean, otherwise a matrix of $n$ rows of random vectors.

References

Devroye, L. (1986). Non-Uniform Random Variate Generation. Springer-Verlag, New York.

See Also

rnorm

Examples

# covariance parameters
Sigma <- matrix(c(10,3,3,2), ncol = 2)
Sigma

# generate the sample
y <- rmnorm(n = 1000, Sigma = Sigma)
var(y)

# scatterplot of a random bivariate normal sample with mean
# vector zero and covariance matrix 'Sigma'
par(pty = "s")
plot(y, xlab = "", ylab = "")
title("bivariate normal sample", font.main = 1)

# QQ-plot of transformed distances
z <- WH.normal(y)
par(pty = "s")
qqnorm(z, main = "Transformed distances QQ-plot")
abline(c(0,1), col = "red", lwd = 2, lty = 2)

Description

Random vector generation uniformly on the sphere.

Usage

rsphere(n = 1, p = 2)

Arguments

Details

The function rsphere is an interface to C routines, which make calls to subroutines from BLAS.

Value

If n = 1 a p-dimensional vector, otherwise a matrix of n rows of random vectors.

References

Devroye, L. (1986). Non-Uniform Random Variate Generation. Springer-Verlag, New York.

See Also

runif

Examples

# generate the sample
z <- rsphere(n = 200)

# scatterplot of a random sample of 200 points uniformly distributed
# on the unit circle
par(pty = "s")
plot(z, xlab = "x", ylab = "y")
title("200 points on the circle", font.main = 1)

Description

Compute the scaled condition number of a rectangular matrix.

Usage

scaled.condition(x, scales = FALSE)

Arguments

Value

The columns of a rectangular matrix x are equilibrated (but not centered), then the scaled condition number is computed following the guidelines of Belsley (1990). If requested, the column scalings are returned as the attribute 'scales'.

References

Belsley, D.A. (1990). Conditioning Diagnostics: Collinearity and Weak Data in Regression. Wiley, New York.

Examples

x <- matrix(c(1, 1, 1,
              1, 2, 1,
              1, 3, 1,
              1, 1,-1,
              1, 2,-1,
              1, 3,-1), ncol = 3, byrow = TRUE)
scaled.condition(x)

Description

Gauss-Seidel method is an iterative algorithm for solving a system of linear equations.

Usage

seidel(a, b, start, maxiter = 200, tol = 1e-7)

Arguments

Details

Let $\bold{D}$, $\bold{L}$, and $\bold{U}$ denote the diagonal, lower triangular and upper triangular parts of a matrix $\bold{A}$. Gauss-Seidel method solve the equation $\bold{Ax} = \bold{b}$, iteratively by rewriting $(\bold{L} + \bold{D})\bold{x} + \bold{Ux} = \bold{b}$. Assuming that $\bold{L} + \bold{D}$ is nonsingular leads to the iteration formula

$\bold{x}^{(k+1)} = -(\bold{L} + \bold{D})^{-1}\bold{U}\bold{x}^{(k)} + (\bold{L} + \bold{D})^{-1}\bold{b}$