Package 'spam'

Description

spam is a collection of functions for sparse matrix algebra.

Gereral overview

What is spam and what is it not:

While Matrix seems an overshoot of classes and SparseM focuses mainly on regression type problem, we provide a minimal set of sparse matrix functions fully functional for everyday spatial statistics life. There is however some emphasize on Markov chain Monte Carlo type calculations within the framework of (Gaussian) Markov random fields.

Emphasis is given on a comprehensive, simple, tutorial structure of the code. The code is S4 based but (in a tutorial spirit) the functions are in a S3 structure visible to the user (exported via NAMESPACE).

There exist many methods for sparse matrices that work identically as in the case of ordinary matrices. All the methods are discussed in the help and can be accessed directly via a *.spam concatenation to the function. For example, help(chol.spam) calls the help directly. We deliberately avoided aliases according to analogue helps from the base package.

Sparseness is used when handling large matrices. Hence, care has been used to provide efficient and fast routines. Essentially, the functions do not transform the sparse structure into full matrices to use standard (available) functionality, followed by a back transform. We agree, more operators, functions, etc. should eventually be implemented.

The packages fields and spam are closely linked.

Author(s)

Reinhard Furrer, with the help of Florian Gerber, Kaspar Moesinger and many others.
Some Fortran routines were written by Youcef Saad, Esmond G. Ng, Barry W. Peyton, Joseph W.H. Liu, Alan D. George.

References

Reinhard Furrer, Stephan R. Sain (2010). "spam: A Sparse Matrix R Package with Emphasis on MCMC Methods for Gaussian Markov Random Fields.", Journal of Statistical Software, 36(10), 1-25, doi:10.18637/jss.v036.i10.
Florian Gerber, Reinhard Furrer (2015). "Pitfalls in the Implementation of Bayesian Hierarchical Modeling of Areal Count Data: An Illustration Using BYM and Leroux Models.", Journal of Statistical Software, Code Snippets, 63(1), 1-32, doi:10.18637/jss.v063.c01.
F. Gerber, K. Moesinger, R. Furrer (2017), "Extending R packages to support 64-bit compiled code: An illustration with spam64 and GIMMS NDVI3g data.", Computer & Geoscience 104, 109-119, doi:10.1016/j.cageo.2016.11.015."

See Also

See spam.class for a detailed class description, spam and spam.ops for creation, coercion and algebraic operations. options.

Examples

## Citations:
citation('spam')
citation('spam', auto=TRUE)


## History of changes
## Not run: 
file.show(system.file("NEWS.md", package = "spam"))

## End(Not run)

Description

Constructing the adjacency graph of the administrative districts of Germany

Usage

adjacency.landkreis(loc)

Arguments

Details

The function is included as an example on how to construct adjacency matrices form a (common) adjacency structure. For the particular example, note that the nodes are not numbered consecutively and that they start from zero.

Value

a sparse matrix in spam format.

Author(s)

Reinhard Furrer

References

The adjacency data has been provided by Havard Rue and is also available in INLA.

See Also

germany.plot super-seeding map.landkreis for plotting.
Oral.

Examples

## Not run: 
loc <- system.file("demodata/germany.adjacency", package="spam")
display( adjacency.landkreis( loc))

## End(Not run)

Description

Utility to compare two spam objects testing 'near equality'. Depending on the type of difference, comparison is still made to some extent, and a report of the differences is returned.

Usage

## S3 method for class 'spam'
all.equal(target, current, tolerance = .Machine$double.eps^0.5,
    scale = NULL, check.attributes = FALSE,...)

Arguments

Details

Numerical comparisons for scale = NULL (the default) are typically on a relative difference scale unless the target values are close to zero or infinite. Specifically, the scale is computed as the average absolute value of target. If this scale is finite and exceeds tolerance, differences are expressed relative to it; otherwise, absolute differences are used.

If scale is numeric (and positive), absolute comparisons are made after scaling (dividing) by scale. Note that if all of scale is sufficiently close to 1 (specifically, within tolerance), the difference is still reported as being on an absolute scale.

Do not use all.equal.spam directly in if expressions: either use isTRUE( all.equal.spam(...)) or identical if appropriate.

Cholesky decomposition routines use this function to test for symmetry.

A method for matrix-spam objects is defined as well.

There is the additional catch of a zero matrix being represented by one zero element, see ‘Examples’ below.

Value

Either TRUE or a vector of 'mode' "character" describing the differences between target and current.

Author(s)

Reinhard Furrer

See Also

isSymmetric.spam and cleanup.

Examples

obj <- diag.spam(2)
obj[1,2] <- .Machine$double.eps

all.equal( diag.spam(2), obj)

all.equal( t(obj), obj)

all.equal( t(obj), obj*1.1)

# We can compare a spam to a matrix
all.equal(diag(2),diag.spam(2))

# the opposite does often not make sense,
# hence, it is not implemented.
all.equal(diag.spam(2),diag(2))


# A zero matrix contains one element:
str(spam(0))
# hence
all.equal.spam(spam(0,3,3), diag.spam(0,3) )
norm(spam(0,3,3) - diag.spam(0,3) )

Description

Returns a vector or array or list of values obtained by applying a function to margins of a sparse matrix.

Usage

apply.spam(X, MARGIN=NULL, FUN, ...)

Arguments

Details

This is a handy wrapper to apply a function to the (nonzero) elements of a sparse matrix. For example, it is possible to apply a covariance matrix to a distance matrix obtained by nearest.dist, see Examples.

A call to apply only coerces the sparse matrix to a regular one.

The basic principle is applying the function to @entries, or to the extracted columns or rows ([,i,drop=F] or [i,,drop=F]). It is important to note that an empty column contains at least one zero value and may lead to non intuitive results.

This function may evolve over the next few releases.

Value

Similar as a call to apply with a regular matrix. The most important cases are as follows. The result is a vector (MARGIN is length 1 and FUN is scalar) or a matrix (MARGIN is length 1 and FUN returns fixed length vectors, or MARGIN is length 2 and FUN is scalar) or a list (if FUN returns vectors of different lengths).

Author(s)

Reinhard Furrer

See Also

base:apply for more details on Value.

Examples

S <- as.spam(dist(1:5))
S <- apply.spam(S/2, NULL, exp)
# instead of 
# S@entries <- exp( S@entries/2) 

# Technical detail, a null matrix consists
# of one zero element.
apply.spam(S,c(1,2),pmax)
apply.spam(S,1,range)


# A similar example as for the base apply.
# However, no dimnames else we would get warnings. 
x <- as.spam(cbind(x1 = 3, x2 = c(0,0,0, 5:2)))
apply.spam(x, 2, mean, trim = .2)
col.sums <- apply.spam(x, 2, sum)
row.sums <- apply.spam(x, 1, sum)
rbind(cbind(x, row.sums), c(col.sums, sum(col.sums)))

apply.spam(x, 2, is.vector)

# Sort the columns of a matrix
# Notice that the result is a list due to the different
# lengths induced by the nonzero elements
apply.spam(x, 2, sort)

# Function with extra args:
cave <- function(x, c1, c2) c(mean(x[c1]), mean(x[c2]))
apply(x,1, cave,  c1=1, c2=c(1,2))

ma <- spam(c(1:4, 0, 0,0, 6), nrow = 2)
ma
apply.spam(ma, 1, table)  #--> a list of length 2
apply.spam(ma, 1, stats::quantile)# 5 x n matrix with rownames

Description

Returns the lower and upper bandwidth of a sparse matrix

Usage

bandwidth(A)

Arguments

Details

The matrix does not need to be diagonal. Values can be negative indicating the the matrix contains a band cinfined in the upper or lower triangular part.

Value

Integer vector containing the lower and upper bandwidth

Author(s)

Reinhard Furrer

See Also

diag.spam.

Examples

bandwidth(spam(c(0, 1), 3, 2))

bandwidth(spam(c(0, 0, 1, rep(0, 9)), 4, 3))

Description

Creates a sparse block-diagonal matrix.

Usage

bdiag.spam(...)

Arguments

Details

This is a small helper function to create block diagonal sparse matrices. In the two matrix case, bdiag.spam(A,B), this is equivalent to a complicated rbind(cbind(A, null), cbind(B, t(null))), where null is a null matrix of appropriate dimension.

It is recursively defined.

The arrays are coerced to sparse matrices first.

This function is similar to the function bdiag from the package Matrix. It is also similar to the function adiag from the package magic. However, here no padding is done and all the dimnames are stripped.

Value

Returns a spam matrix as described above.

Author(s)

Reinhard Furrer

See Also

diag.spam.

Examples

A <- diag.spam(2, 4)           # 2*I4
B <- matrix(1,3,3)
AB <- bdiag.spam(A,B)

# equivalent to:
ABalt <- rbind(cbind( A, matrix(0,nrow(A),ncol(B))),
               cbind( matrix(0,nrow(B),ncol(A)), B))
         
norm(AB-ABalt)


# Matrices do not need to be square:
bdiag.spam(1,2:5,6)

Description

Take a sequence of vector, matrix or spam object arguments and combine by columns or rows, respectively.

Usage

# cbind(\dots, force64 = getOption("spam.force64"), deparse.level = 0)
# rbind(\dots, deparse.level = 0)

Arguments

Details

rbind and cbind are not exactly symmetric in how the objects are processed. cbind calls a Fortran routine after the input has been coerced to spam objects. Whereas rbind calls a Fortran routine only in the case of to spam matrices. Note that row binding is essentially an concatenation of the slots due to the sparse storage format.

Only two objects at a time are processed. If more than two are present, a loop concatenates them successively.

A method is defined for a spam object as first argument.

Value

a spam object combining the ... arguments column-wise or row-wise. (Exception: if there are no inputs or all the inputs are NULL, the value is NULL.)

Author(s)

Reinhard Furrer

Examples

x <- cbind.spam(1:5,6)

y <- cbind(x, 7)

rbind( x, x)
# for some large matrices   t( cbind( t(x), t(x)))
# might be slightly faster:

Description

chol performs a Cholesky decomposition of a symmetric positive definite sparse matrix x of class spam.

Usage

# chol(x, \dots)

## S4 method for signature 'spam'
chol(x, pivot = "MMD", method = "NgPeyton",
       memory = list(), eps = getOption("spam.eps"), Rstruct=NULL,
       ..., verbose=FALSE)

# update.spam.chol.NgPeyton(object, x,...)
## S4 method for signature 'spam.chol.NgPeyton'
update(object, x,...)

Arguments

Details

chol performs a Cholesky decomposition of a symmetric positive definite sparse matrix x of class spam. Currently, there is only the block sparse Cholesky algorithm of Ng and Peyton (1993) implemented (method="NgPeyton").

To pivot/permute the matrix, you can choose between the multiple minimum degree (pivot="MMD") or reverse Cuthill-Mckee (pivot="RCM") from George and Lui (1981). It is also possible to furnish a specific permutation in which case pivot is a vector. For compatibility reasons, pivot can also take a logical in which for FALSE no permutation is done and for TRUE is equivalent to MMD.

Often the sparsity structure is fixed and does not change, but the entries do. In those cases, we can update the Cholesky factor with update.spam.chol.NgPeyton by suppling a Cholesky factor and the updated matrix. For U <- chol(A), update(U, Anew) and chol(Anew, Rstruct=U) are equivalent.

The option cholupdatesingular determines how singular matrices are handled by update. The function hands back an error ("error"), a warning ("warning") or the value NULL ("null").

The Cholesky decompositions requires parameters, linked to memory allocation. If the default values are too small the Fortran routine returns an error to R, which allocates more space and calls the Fortran routine again. The user can also pass better estimates of the allocation sizes to chol with the argument memory=list(nnzR=..., nnzcolindices=...). The minimal sizes for a fixed sparsity structure can be obtained from a summary call, see ‘Examples’.

The output of chol can be used with forwardsolve and backsolve to solve a system of linear equations.

Notice that the Cholesky factorization of the package SparseM is also based on the algorithm of Ng and Peyton (1993). Whereas the Cholesky routine of the package Matrix are based on CHOLMOD by Timothy A. Davis (C code).

Value

The function returns the Cholesky factor in an object of class spam.chol.method. Recall that the latter is the Cholesky factor of a reordered matrix x, see also ordering.

Note

Although the symmetric structure of x is needed, only the upper diagonal entries are used. By default, the code does check for symmetry (contrarily to base:::chol). However, depending on the matrix size, this is a time consuming test. A test is ignored if options("spam.cholsymmetrycheck") is set to FALSE.

If a permutation is supplied with pivot, options("spam.cholpivotcheck") determines if the permutation is tested for validity (defaults to TRUE).

Author(s)

Reinhard Furrer, based on Ng and Peyton (1993) Fortran routines

References

Ng, E. G. and Peyton, B. W. (1993) Block sparse Cholesky algorithms on advanced uniprocessor computers, SIAM J. Sci. Comput., 14, 1034–1056.

Gilbert, J. R., Ng, E. G. and Peyton, B. W. (1994) An efficient algorithm to compute row and column counts for sparse Cholesky factorization, SIAM J. Matrix Anal. Appl., 15, 1075–1091.

George, A. and Liu, J. (1981) Computer Solution of Large Sparse Positive Definite Systems, Prentice Hall.

See Also

det.spam, solve.spam, forwardsolve.spam, backsolve.spam and ordering.

Examples

# generate multivariate normals:
set.seed(13)
n <- 25    # dimension
N <- 1000  # sample size
Sigma <- .25^abs(outer(1:n,1:n,"-"))
Sigma <- as.spam( Sigma, eps=1e-4)

cholS <- chol( Sigma)
# cholS is the upper triangular part of the permutated matrix Sigma
iord <- ordering(cholS, inv=TRUE)

R <- as.spam(cholS)
mvsample <- ( array(rnorm(N*n),c(N,n)) %*% R)[,iord]
# It is often better to order the sample than the matrix
# R itself.

# 'mvsample' is of class 'spam'. We need to transform it to a
# regular matrix, as there is no method 'var' for 'spam' (should there?).
norm( var( as.matrix( mvsample)) - Sigma, type='m')
norm( t(R) %*% R - Sigma)


# To speed up factorizations, memory allocations can be optimized:
opt <- summary(cholS)
# here, some elements of Sigma may be changed...
cholS <- chol( Sigma, memory=list(nnzR=opt$nnzR,nnzcolindices=opt$nnzc))

Description

Creates a circulant matrix in spam format.

Usage

circulant.spam(x, n = NULL, eps = getOption("spam.eps"))

Arguments

Value

The circulant matrix in spam format.

Author(s)

Reinhard Furrer

See Also

circulant from package magic, toeplitz.spam

Examples

circulant.spam(c(1,.25,0,0,0))

Description

Eliminates an zeros in a sparse matrix.

Usage

cleanup(x, eps = getOption("spam.eps"))

Arguments

Details

A sparse matrix may still contain zeros. This function (aliased to as.spam) filters these values.
This often causes confusion when testing such matrices for symmetry or comparing apparently equal matrices with all.equal (see ‘Examples’ below.

Author(s)

Reinhard Furrer

See Also

isSymmetric.spam and all.equal.spam.

Examples

A <- diag.spam(2)
A[1,2] <- 0

all.equal(A, t(A))
isSymmetric.spam(A)
all.equal(cleanup(A), diag.spam(2))

Description

These functions manage the relations that allow coercing a spam object to a given class.

Methods

signature(from = "spam", to = "matrix")

this is essentially equivalent to as.matrix(object).

signature(from = "spam", to = "list")

this is essentially equivalent to triplet(object).

signature(from = "spam", to = "vector")

this is essentially equivalent to object@entries (structurebased=TRUE) or c(object).

signature(from = "spam", to = "logical")

the entries are forced to logicals (nonzeros only in case of structurebased=TRUE).

signature(from = "spam", to = "integer")

the entries are forced to integers (nonzeros only in case of structurebased=TRUE).

Examples

ifelse( diag.spam(2)*c(0,1), TRUE, FALSE)

Description

Coercion of spam matrices to proper vector objects

Usage

## S4 method for signature 'spam'
as.vector(x, mode = "any")

Arguments

Details

This coercion allows smooth transitions between different matrix formats, see example below.
The Cholesky factors are first transformed to a spam object.

Value

If structurebased=TRUE, the vector x@entries.
Conversely, if structurebased=FALSE, the result is identical to one with as.vector(as.matrix(x)).

Author(s)

Reinhard Furrer

See Also

spam.options

Examples

x <- diag(2)
ifelse( x, x, 1-x)
ifelse( x, as.vector(x), 1-as.vector(x))

x <- diag.spam(2)
options(spam.structurebased=FALSE)
ifelse( x, as.vector(x), 1-as.vector(x))
options(spam.structurebased=TRUE)
ifelse( x, as.vector(x), 1-as.vector(x))

Description

A few results of computational complexities for selected sparse algoritms in spam

Details

A Cholesky factorization of an n-matrix requires n^3/3 flops. In case of banded matrices (bandwidth p, p<<n) a factorization requires about 2np^2 flops. Forward- and backsolves for banded matrices require essentially 2np flops.

George and Liu (1981) proves that any reordering would require at least O(n^3/2) flops for the factorization and produce at least O(n log(n)) fill-ins for square lattices with a local neighbor hood.
They also show that algorithms based on nested dissection are optimal in the order of magnitude sense.

More to follow.

References

George, A. and Liu, J. (1981) Computer Solution of Large Sparse Positive Definite Systems, Prentice Hall.

See Also

det, solve, forwardsolve, backsolve and ordering.

Description

Modify slots of spam objects

Usage

rowpointers( x) <- value
colindices( x) <- value
entries( x) <- value

Arguments

Details

Various tests are performed. Thus much slower than direct assignment.
Slot dimension should be changed through pad or dim

Value

Modified spam object.

Author(s)

Reinhard Furrer

Examples

x <- diag.spam( 2)  
rowpointers( x) <- c(1,1,3)

# The last line is equivalent to 
x@rowpointers <- as.integer( c(1,1,3))

Description

Evaluate a covariance function.

Usage

covmat(h, theta, ... , type="sph")

cov.exp(h, theta, ... , eps= getOption("spam.eps"))
cov.sph(h, theta, ... , eps= getOption("spam.eps"))
cov.nug(h, theta, ... , eps= getOption("spam.eps"))
cov.wend1(h, theta, ... , eps= getOption("spam.eps"))
cov.wend2(h, theta, ... , eps= getOption("spam.eps"))
cov.wu1(h, theta, ... , eps= getOption("spam.eps"))
cov.wu2(h, theta, ... , eps= getOption("spam.eps"))
cov.wu3(h, theta, ... , eps= getOption("spam.eps"))
cov.mat(h, theta, ... , eps= getOption("spam.eps"))
cov.finnmat(h, theta, ... , eps= getOption("spam.eps"))
cov.mat12(h, theta, ... , eps= getOption("spam.eps"))
cov.mat32(h, theta, ... , eps= getOption("spam.eps"))
cov.mat52(h, theta, ... , eps= getOption("spam.eps"))

cor.sph(h, range, ... , eps= getOption("spam.eps"))

Arguments

Details

covmat is a wrapper that calls the other functions according to the argument type. The nomenclature is similar to precmat.
The parametrization is (range, [partial-sill = 1], [smoothness = 1], [nugget = 0]), where only the range needs to be specified. In case of negative parameter values, a warning is issued and the absolute value is retained. Although more cryptic, having all arguments as a single vector simplifies optimization with optim.
The parameters are and locations are up to precision epsilon. That means that all distances smaller than eps are considered zero; a nugget smaller than eps is ignored; a range smaller than eps represents a nugget model; etc.
cov.finnmat() is similar to cov.mat() but with the sqrt(8*smoothness)/range argument in the Bessel function (instead of 1/range). cov.mat12() is a wrapper to cov.exp() cov.mat32(), and cov.mat52() are fast version of cov.mat() with smoothness 3/2 and 5/2, respectively (factor 10).
cor.sph(,range) is a fast version of cov.sph(,c(range,1,0)).
Currently, the functions distinguish between a sparse spam object h and any other numeric type. In the future, this might change and appropriate methods will be implemented.

Value

Covariance function evaluated on h.

Author(s)

Reinhard Furrer

References

Any classical book about geostatistics.

See Also

precmat.

Examples

set.seed(123)
n <- 200
locs <- cbind(runif(n),runif(n))
h <- nearest.dist(locs, delta=sqrt(2), upper = NULL)
Sigma <- cov.sph(h, c(.3, 1, .1))

iidsample <- rnorm(n)
cholS <- chol.spam(as.spam(Sigma))
iorder <- iord <- ordering(cholS, inv = TRUE)
sample <- (iidsample %*% as.spam(cholS))[iorder]
plot(locs, col = fields::tim.colors(n = 256)[cut(sample, n)], pch = 20)

## Not run: 
h <- seq(0, to=1, length.out=100)
plot( h, cov.exp(h, c(1/3,1)), type='l', ylim=c(0,1))
type <- c("sph","wendland1","wendland2","wu1","wu2","wu3")
for (i in 1:6)
  lines( h, covmat(h, 1, type=type[i]), col=i+1)
legend('topright',legend=type, col=2:7, lty=1)


## End(Not run)

Description

Given matrices x and y as arguments, return a matrix cross-product. This is formally equivalent to (but usually slightly faster than) the call t(x) %*% y (crossprod.spam) or x %*% t(y) (tcrossprod.spam).

Usage

crossprod.spam(x, y = NULL, ...)
     
     tcrossprod.spam(x, y = NULL, ...)

Arguments

Value

A double matrix

Note

When x or y are not matrices, they are treated as column or row matrices.

Author(s)

Reinhard Furrer

Examples

crossprod.spam(diag.spam(2),1:2)

Description

det and determinant calculate the determinant of a symmetric, positive definite sparse matrix. determinant returns separately the modulus of the determinant, optionally on the logarithm scale, and the sign of the determinant.

Usage

det(x, ...)
determinant(x, logarithm = TRUE, ...)

Arguments

Details

If the matrix is not positive definite, the function issues a warning and returns NA.

The determinant is based on the product of the diagonal entries of a Cholesky factor, i.e. internally, a Cholesky decomposition is performed. By default, the NgPeyton algorithm with minimal degree ordering us used. To change the methods or supply additonal parameters to the Cholesky factorization function, it is possible to pass via chol.

The determinant of a Cholesky factor is also defined.

Value

For det, the determinant of x. For determinant, a list with components

Author(s)

Reinhard Furrer

References

Ng, E. G. and B. W. Peyton (1993) Block sparse Cholesky algorithms on advanced uniprocessor computers, SIAM J. Sci. Comput., 14, 1034–1056.

See Also

chol.spam

Examples

x <- spam( c(4,3,0,3,5,1,0,1,4), 3)
det( x)
determinant( x)

det( chol( x))

Description

Extract or replace the diagonal of a matrix, or construct a diagonal matrix.

Usage

## diag(x)
## diag(x=1, nrow, ncol, names = TRUE)
diag(x) <- value

diag.spam(x=1, nrow, ncol)
spam_diag(x=1, nrow, ncol)
diag.spam(x) <- value

Arguments

Details

Using diag(x) can have unexpected effects if x is a vector that could be of length one. Use diag(x, nrow = length(x)) for consistent behaviour.

Value

If x is a spam matrix then diag(x) returns the diagonal of x.

The assignment form sets the diagonal of the sparse matrix x to the given value(s).

diag.spam works as diag for spam matrices: If x is a vector (or 1D array) of length two or more, then diag.spam(x) returns a diagonal matrix whose diagonal is x. spam_diag is an alias for diag.spam and in the spirit of the result of diag is a spam object.

If x is a vector of length one then diag.spam(x) returns an identity matrix of order the nearest integer to x. The dimension of the returned matrix can be specified by nrow and ncol (the default is square).

The assignment form sets the diagonal of the matrix x to the given value(s).

Author(s)

Reinhard Furrer

See Also

upper.tri, lower.tri.

Examples

diag.spam(2, 4)           # 2*I4
smat <- diag.spam(1:5)
diag( smat)
diag( smat) <- 5:1

# The last line is equivalent to 
diag.spam( smat) <- 5:1

# Note that diag.spam( 1:5) <- 5:1 not work of course.

Description

Returns suitably lagged and iterated differences.

Usage

# diff.spam(x, lag = 1, differences = 1, ...)
## S4 method for signature 'spam'
diff(x, lag = 1, differences = 1, ...)

Arguments

Value

A spam matrix with elements similar to as.spam(diff(as.matrix(x), ...)).

Author(s)

Reinhard Furrer

See Also

diff in base, precmat.

Examples

# incidence matrix for a RW(3) model
D <- diff.spam(diag.spam(10), lag=1, differences=3)
t(D)%*%D

Description

Retrieve or set the dimension of an spam object.

Usage

# dim(x)
# dim(x) <- value

Arguments

Details

In older version of spam, the behavior of the replacement method was different and is now implemented in pad.spam.

Value

dim retrieves the dimension slot of the object. It is a vector of mode integer.

The replacemnt method changes the dimension of the object by rearranging.

Author(s)

Reinhard Furrer

See Also

pad.spam.

Examples

x <- diag(4)
dim(x)<-c(2,8)
x

s <- diag.spam(4)
dim(s) <- c(2,8)  # result is different than x

s <- diag.spam(4)
pad(s) <- c(7,3)  # any positive value can be used

Description

The function represents the nonzero entries in a simple bicolor plot.

Usage

display(x, ...)

Arguments

Details

spam.getOption("imagesize") determines if the sparse matrix is coerced into a matrix and the plotted with image.default or if the matrix is simply represented as a scatterplot with pch=".". The points are scaled according to cex*getOption("spam.cex")/(nrow + ncol). For some devices or for non-square matrices, cex needs probably some adjustment.

Author(s)

Reinhard Furrer

See Also

image, spam.options

Examples

set.seed(13)

smat <- spam_random(8)
par(mfcol=c(1,2), pty='s')
options(spam.imagesize = 1000)
display(smat)
options(spam.imagesize = 10)
display(smat, cex=.25)


# very large but very sparse matrix
smat <- spam_random(2^14, distribution=rnorm, density=1e-5, verbose=TRUE)
par(mfcol=c(1, 1), mai=c(.4,.4,.1,.1), pty='s')
display(smat)

Description

Functions to calculate eigenvalues and eigenvectors of sparse matrices. It uses the value of spam.options("inefficiencywarning") to dispatch between base::eigen() or the Implicitly Restarted Arnoldi Process, using 'ARPACK'.

eigen.spam is a wrapper function of eigen_approx and transforms its output to base::eigen like.

Usage

eigen.spam(x, nev = 10, symmetric, only.values = FALSE, control = list())
eigen_approx(x, nev, ncv, nitr, mode, only.values = FALSE, verbose = FALSE, f_routine)

Arguments

Details

mode = " ":

there are different modes available for this function, each mode returns a different range of eigenvalues. Also the available modes are dependent, whether the input matrix is symmetric or not:

"LM":

Eigenvalues with largest magnitude (sym, non sym), that is, largest eigenvalues in the Euclidean norm of complex numbers.

"SM":

Eigenvalues with smallest magnitude (sym, non sym), that is, smallest eigenvalues in the Euclidean norm of complex numbers.

"LR":

Eigenvalues with largest real part (non sym).

"SR":

Eigenvalues with smallest real part (non sym).

"LI":

Eigenvalues with largest imaginary part (non sym).

"SI":

Eigenvalues with smallest imaginary part (non sym).

"LA":

Eigenvalues with largest algebraic value (sym), that is, largest eigenvalues inclusive of any negative sign.

"SA":

Eigenvalues with smallest algebraic value (syn), that is, smallest eigenvalues inclusive of any negative sign.

ncv:

the largest number of basis vectors that will be used in the Implicitly Restarted Arnoldi Process. Work per major iteration is proportional to x@dimension[1]*ncv*ncv. The default is set if symmetric to min(x@dimension[1] + 1, max(2 * nev + 1, 200)) or else to min(x@dimension[1] - 1, max(2 * nev + 1, 100)). Note, this value should not be chosen arbitrary large, but slightly larger than nev. Otherwise it could lead to memory allocation problems.

nitr:

the maximum number of iterations. The default is set to ncv + 1000

spamflag = FALSE:

if TRUE, the Implicitly Restarted Arnoldi Process is used, independent of the dimension of the respective matrix (provided matrix is larger than 10x10).

verbose = FALSE:

print additional information.

cmplxeps:

threshold to determine whether a double value is zero, while transforming the ARPACK output to R class complex. The default is set to .Machine$double.eps.

Value

A vector of the length corresponding to the dimension of the input matrix. Containing the required nev eigenvalues. If requested also the corresponding eigenvectors. In the non symmetric case, the eigenvalues are returned in a matrix with a column containing the real parts and a column containing the imaginary parts of the eigenvalues. The eigenvectors are then returned in two matrices.

Note

The user is advised to choose the control options carefully, see ‘Details’ for more information.

Author(s)

Roman Flury, Reinhard Furrer

References

Lehoucq, R. B. and Sorensen, D. C. and Yang, C. (1997) ARPACK Users Guide: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods.

See Also

Option "inefficiencywarning" in spam.options and spam_random.

Examples

set.seed(81)
rspam <- spam_random(42^2, density = .0001, spd = TRUE)
SPD <- eigen.spam(rspam, nev = 18, control = list(mode = "SM"),
                  only.values = TRUE)

any(SPD$values <= 0, na.rm = TRUE)
isSymmetric(rspam)
# hence the matrix is symmetric positiv definit

rspam2 <- spam_random(50^2, density = .0001, spd = FALSE, sym = TRUE,
                      distribution = rpois, lambda = 2)
SNPD <- eigen.spam(rspam2, nev = 18, control = list(mode = "SM"),
                    only.values = TRUE)

any(SNPD$values <= 0, na.rm = TRUE)
isSymmetric(rspam2)
# hence the matrix is symmetric but not positiv definit

Description

These functions are simple wrappers to nearest.dist to be used in fields.

Usage

spam_rdist( x1, x2, delta = 1)
spam_rdist.earth( x1, x2, delta = 1, miles=TRUE, R=NULL)

Arguments

Details

These functions are wrappers to rdist and rdist.earth in fields. They are used to simplify the use of sparse matrices in functions like mKrig.

For great circle distance, the matrices x1 and x2 contain the degrees longitudes in the first and the degrees latitudes in the second column. delta is in degrees. Hence to restrict to distances smaller than delta.km, one has to specify delta=delta.km*360/(6378.388*2*pi).

Value

A spam object containing the distances spanned between zero and delta. The sparse matrix may contain many zeros (e.g., collocated data). However, to calculate covariances, these zeros are essential.

Author(s)

Reinhard Furrer

See Also

nearest.dist

Examples

## Not run: 
require(fields)
look <- mKrig(x,Y, Covariance="Wendland", dimension=2, k=1,
    cov.args=list( Distance='spam_rdist'))

## End(Not run)

Description

Transform between the spam sparse format to the matrix.csr format of SparseM and dgRMatrix format of Matrix

Usage

as.spam.matrix.csr(x)
as.dgRMatrix.spam(x)
as.dgCMatrix.spam(x)
as.spam.dgRMatrix(x)
as.spam.dgCMatrix(x)

Arguments

Details

We do not provide any S4 methods and because of the existing mechanism a standard S3 does not work.

The functions are based on require.

Notice that as.matrix.csr.spam should read as as."matrix.csr".spam.

Value

According to the call, a sparse matrix of class spam, matrix.csr, dgRMatrix or dgCMatrix.

Author(s)

Reinhard Furrer

See Also

triplet, Matrix or matrix.csr from package SparseM.

Examples

## Not run: 
S <- diag.spam(4)
R <- as.dgRMatrix.spam( S)
C <- as.dgCMatrix.spam( S)
as.spam.dgCMatrix(C)
slotNames(C)
slotNames(R)
# For column oriented sparse formats a transpose does not the job,
# as the slot names change.


# as.spam(R) does not work.

## End(Not run)

## Not run: 
# for transformations between SparseM and spam:

as.matrix.csr.spam <- function(x,...) {
  if (new("matrix.csr")) {
    newx <- new("matrix.csr")
    slot(newx,"ra",check=FALSE) <- x@entries
    slot(newx,"ja",check=FALSE) <- x@colindices
    slot(newx,"ia",check=FALSE) <- x@rowpointers
    slot(newx,"dimension",check=FALSE) <- x@dimension
    return(newx)
    }
  stop("function requires 'SparseM' package")
}
# then with `SparseM`:  as.matrix.csr.spam( spamobject )

## End(Not run)

## Not run: 
# a dataset contained in Matrix
data(KNex, package='Matrix')
summary( KN <- as.spam.dgCMatrix(KNex$mm) )


## End(Not run)

Description

Supplementary data used for the display of data from the administrative districts of Germany

Format

germany.info is a list with elements

n

544 (number of districts around 1990).

xrep,yrep

representative coordinates of the districts (vectors of length 544)

xlim,ylim

2-vectors defining the limits of the districts.

polyid

linking the polygons to the districts (599 vector).

id

linking the districts to Community Identification Number.

germany.poly defines the polygons. It is a 17965 by two matrix, each polygon separated by a row of NAs, each district by two rows.
germany defines the polygons in form of a list (backwards compatibility).

Details

The representative coordinates are calculated based on the mean value of the polygon coordinates. This creates sometimes strange values, e.g., district Leer.

Author(s)

Reinhard Furrer

References

The meta-data has been constructed based on (essentially) files from the package INLA, see demo(Bym).

See also https://de.wikipedia.org/wiki/Amtlicher_Gemeindeschl%C3%BCssel and https://en.wikipedia.org/wiki/Districts_of_Germany

See Also

germany.plot Oral.

Examples

# Plot the Bundeslaender:
germany.plot(germany.info$id%/%1000,col=rep(2:8,3), legend=FALSE)

Description

Displaying data over the administrative districts of Germany

Usage

germany.plot(vect,  col=NULL, zlim=range(vect), legend=TRUE, 
             main=NULL, cex.axis=1, cex.main=1.5, add=FALSE, ... )

Arguments

Details

The legend is only added, provided (a) function image.plot is available.

The perfect position of the legend is an art per se and depends on various par parameters. One possiblity for finer control is to not plot it and to manually call the function image.plot of fields.

Author(s)

Reinhard Furrer

References

See also https://de.wikipedia.org/wiki/Amtlicher_Gemeindeschl%C3%BCssel and https://de.wikipedia.org/wiki/Liste_der_Landkreise_in_Deutschland

See Also

Oral.

Examples

data( Oral)
germany.plot( Oral$Y/Oral$E)


# Plot the Bundeslaender:
germany.plot(germany.info$id%/%1000,col=rep(2:8,3), legend=FALSE)

Description

Multiplying specific submatrices of a spam matrix with different factors.

Usage

gmult(x, splits, fact)

Arguments

Value

Spam matrix, where each specified submatrix is multiplied with a factor.

Author(s)

Florian Gerber, Roman Flury

Examples

x <- spam(1, 15, 15)
print(x, minimal = FALSE)

splits <- c(1,2,8,ncol(x)+1) # divide matrix into 9 submatrices
fact <- array(1:9, c(3,3))   # multiply each submatrix with a different factor
oF <- gmult(x, splits, fact)
print(oF, minimal = FALSE)

Description

For two (MCMC) chains of the same length trace plots and a scatter plot are drawn.

Usage

grid_trace2(chain1, chain2 = NULL,
            xlim = NULL, ylim1 = NULL, ylim2=NULL,
            chain1_lab = "", chain2_lab = "", main = "",
            chain1_yaxis_at = NULL, chain2_yaxis_at = NULL,
            log = FALSE,
            cex_points = unit(0.2, "mm"),
            cex_main = unit(1.2, "mm"),
            lwd_lines = unit(0.1, "mm"),
            lwd_frame = unit(0.8, "mm"),
            draw = TRUE)

Arguments

Details

The figure is optimized for a plot of the format x11(width = 10, height = 3).

Value

A grob object.

Author(s)

Florian Gerber <[email protected]>

See Also

grid_zoom

Examples

grid_trace2(runif(500), runif(500),
            chain1_yaxis_at = seq(.2, 1, by = .2),
            chain1_lab = "chain1", chain2_lab = "chain2",
            main = "Uniform",
            lwd_lines = grid::unit(.5, "mm"))

Description

This function takes a grob object (e.g. created with package grid) and adds a zoom window.

Usage

grid_zoom(inputGrob = pointsGrob(runif(200),runif(200)),
          inputViewport = viewport(name='main'),
          x = 'topleft', y, just,
          ratio = c(.3,.4), zoom_xlim, zoom_ylim,
          rect = TRUE, rect_lwd = 1, rect_fill = 'gray92',
          draw =TRUE, zoom_fill = 'white',
          zoom_frame_gp = gpar(lwd = 1),
          zoom_gp = NULL, zoom_xaxis = xaxisGrob(main = FALSE),
          zoom_yaxis = NULL)

Arguments

Details

A zoom plot does only make sense if all objects of the inputGrob are specified in native units. Additional caution me be require for certain grobs: e.g. a zoom of a circleGrob() is problematic if the x and y axis are stretched by a different amount.

Value

A grob object.

Author(s)

Florian Gerber <[email protected]>

See Also

grid_trace2

Examples

require(grid)
## -- Example 1 --
set.seed(133)
grid_zoom(inputGrob = pointsGrob(runif(200), runif(200)),
          inputViewport = viewport(name = 'main'),
          zoom_xlim = c(.2, .3), zoom_ylim = c(.2, .3))


## -- Example 2 --
## initial plot
grid.newpage()
vp <- viewport(width=.8, height=.8, clip='on')
gt <- gTree(children=gList(polylineGrob(x=c((0:4)/10, rep(.5, 5), (10:6)/10, rep(.5, 5)),
              y=c(rep(.5, 5), (10:6/10), rep(.5, 5), (0:4)/10),
              id=rep(1:5, 4), default.units='native',
              gp=gpar(col=1:5, lwd=3)),
              pointsGrob(runif(1000), runif(1000),pch='.', gp=gpar(cex=3)),
              rectGrob(gp=gpar(lwd=3))))
pushViewport(vp)
grid.draw(gt)

## plot with zoom window
grid.newpage()
grid_zoom(inputGrob = gt,
          inputViewport = vp,
          x='topright', zoom_xlim=c(.6,.73), zoom_ylim=c(.3,.43),ratio=.4,
          zoom_xaxis = NULL, zoom_gp = gpar(cex=3))

Description

Returns the upper left or lower right part of a spam object.

Usage

## S4 method for signature 'spam'
head(x, n = 6L, m = n, ...)
## S4 method for signature 'spam'
tail(x, n = 6L, m = n, addrownums = TRUE, ...)

Arguments

Details

For matrices, 2-dim tables and data frames, head() (tail()) returns the first (last) n rows and m columns when n > 0 or all but the last (first) n rows when n < 0 (with similar behavior for m).

tail() will add row and column names of the form "[n,]" and "[,n]" to the result, so that it looks similar to the last lines and columns of x when printed. Setting addrownums = FALSE suppresses this behaviour.

A method for spam.chol.NgPeyton objects is exported as well.

Value

An regular matrix.

Author(s)

Reinhard Furrer

Examples

head( precmat.RW2( 10))
tail( precmat.season(n=10, season=3), n=4, m=10)

Description

The function creates a grid of colored rectangles with colors corresponding to the values of the spam matrix.

Usage

## S4 method for signature 'spam'
image(x, cex = NULL, ...)

Arguments

Details

getOption("spam.imagesize") determines if the sparse matrix is coerced into a matrix and the plotted similarly to image.default or if the matrix is simply represented as a scatterplot with pch=".". The points are scaled according to cex*getOption("spam.cex")/(nrow+ncol). For some devices or for non-square matrices, cex needs probably some adjustment.
The a zero matrix in spam format has as (1,1) entry the value zero and only missing entries are interpreted as NA in the scatter plot.

Author(s)

Reinhard Furrer

See Also

display and spam.options.
The code is based on image of graphics.

Examples

set.seed(13)

smat <- spam_random(8)

par(mfcol=c(1,2),pty='s')
options(spam.imagesize=1000)
image(smat) # or use better color schemes
options(spam.imagesize=10)
image(smat, cex=.25)

smat <- spam_random(2^14, distribution=rnorm, density=1e-5, verbose=TRUE)
par(mfcol=c(1,1), mai=c(.4,.4,.1,.1), pty='s')
image(smat)

Description

Read matrices stored in the Harwell-Boeing or MatrixMarket formats.

Usage

read.HB(file)
read.MM(file)

Arguments

Details

The names of files storing matrices in the Harwell-Boeing format usually end in ".rua" or ".rsa". Those storing matrices in the MatrixMarket format usually end in ".mtx".

Currently, only real assembled Harwell-Boeing can be read with read.HB. Reading MatrixMarket formats is more flexible. However, as entries of spam matrices are of mode double, integers matrices are coerced to doubles, patterns lead to matrices containing ones and complex are coerced to the real part thereof. In these aforementioned cases, a warning is issued.

MatrixMarket also defines an array format, in which case a (possibly) dense spam object is return (retaining only elements which are larger than options('spam.eps'). A warning is issued.

Value

A sparse matrix of class spam.

Note

The functions are based on readHB and readMM from the library Matrix to build the connection and read the raw data. At present, read.MM is more flexible than readMM.

Author(s)

Reinhard Furrer based on Matrix functions by Douglas Bates [email protected] and Martin Maechler [email protected]

References

https://math.nist.gov/MatrixMarket/

https://sparse.tamu.edu/

Examples

## Not run: 
image(read.MM(gzcon(url(
  "ftp://math.nist.gov/pub/MatrixMarket2/Harwell-Boeing/bcspwr/bcspwr01.mtx.gz"))))

## End(Not run)

## Not run: 
## Datasets supplied within Matrix
str(read.MM(system.file("external/pores_1.mtx",package = "Matrix")))
str(read.HB(system.file("external/utm300.rua", package = "Matrix")))
str(read.MM(system.file("external/lund_a.mtx", package = "Matrix")))
str(read.HB(system.file("external/lund_a.rsa", package = "Matrix")))

## End(Not run)

Description

Efficient function to test if 'object' is symmetric or not.

Usage

# isSymmetric.spam(object, ...)
## S3 method for class 'spam'
isSymmetric(object, tol = 100 * .Machine$double.eps, ...)

Arguments

Details

symmetry is assessed by comparing the sparsity structure of object and t(object) via the function all.equal.spam. If a difference is detected, the matrix is cleaned with cleanup and compared again.

Value

logical indicating if object is symmetric or not.

Author(s)

Reinhard Furrer

See Also

all.equal.spam, cleanup.

Examples

obj <- diag.spam(2)
isSymmetric(obj)

obj[1,2] <- .Machine$double.eps
isSymmetric(obj)
all.equal(obj, t(obj))

Description

Computes the generalised kronecker product of two arrays, X and Y.

Usage

kronecker.spam(X, Y, FUN = "*", make.dimnames = FALSE, ...)

Arguments

Details

The sparsity structure is determined by the ordinary %x%. Hence, the result of kronecker(X, Y, FUN = "+") is different depending on the input.

Value

An array A with dimensions dim(X) * dim(Y).

Author(s)

Reinhard Furrer

Examples

# Starting with non-spam objects, we get a spam matrix
kronecker.spam( diag(2), array(1:4, c(2, 2)))

kronecker( diag.spam(2), array(1:4, c(2, 2)))

# Notice the preservation of sparsity structure:
kronecker( diag.spam(2), array(1:4, c(2, 2)), FUN="+")

Description

The R package spam can handle sparse matrices with up to 2^31-1 non-zero elements. For matrices with more non-zero elements it is necessary to load the spam64 package in addition.

Details

With the help of the R package dotCall64 spam interfaces either the compiled code with 32-bit integers provided in spam or the compiled code with 64-bit integers provided in spam64.
To mimick 64-bit behavior, set options(spam.force64 = TRUE). The subsequent matrix indices are then stored in double format.

Author(s)

Reinhard Furrer, Florian Gerber, Kaspar Moesinger, Daniel Gerber

References

F. Gerber, K. Moesinger, R. Furrer (2017), Extending R packages to support 64-bit compiled code: An illustration with spam64 and GIMMS NDVI3g data, Computer & Geoscience 104, 109-119, https://doi.org/10.1016/j.cageo.2016.11.015.

See Also

spam64-package, dotCall64.

Examples

## Not run: 
## the following matrices are very large, and hence,
## require much memory and cpu time.
library("spam64")
s1 <- spam(1, ncol=2^30)        # 32-bit matrix
s1

s2 <- cbind(s1, s1)             # 64-bit matrix
s2

s3 <- spam(1, ncol=2^31)        # 64-bit matrix
s3

## End(Not run)

Description

Returns the lower or upper triangular structure or entries of a sparse matrix.

Usage

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

Arguments

Details

Often not only the structure of the matrix is required but the entries as well. For compatibility, the default is only a structure consisting of ones (representing TRUEs). Setting the flag getOption( "spam.trivalues") to TRUE, the function returns the actual values.

See Also

spam.options and diag

Examples

smat <- spam( c( 1,2,0,3,0,0,0,4,5),3)
upper.tri( smat)
upper.tri( smat, diag=TRUE)

options(spam.trivalues=TRUE)
upper.tri( smat)

Description

Creates precision matrices for gridded GMRF.

Usage

precmat.GMRFreglat(n,m, par, model = "m1p1",  eps = getOption("spam.eps"))

Arguments

Details

The function should be illustrative on how to create precision matrices for gridded GMRF. Hence, no testing (positive definiteness is done).

The model specification "m" determines the complexity and "p" the number of parameters.

Please see the examples on the meaning of the different models.

Value

A spam matrix of dimension prod(dims)xprod(dims).

Author(s)

Reinhard Furrer

See Also

precmat, toeplitz.spam, kronecker.spam

Examples

as.matrix(precmat.GMRFreglat(4, 3, c(.4),         'm1p1'))
as.matrix(precmat.GMRFreglat(4, 3, c(.4,.3),      'm1p2'))
as.matrix(precmat.GMRFreglat(4, 3, c(.4,.3,.2),   'm2p3'))
as.matrix(precmat.GMRFreglat(4, 3, c(.4,.3,.2,.1),'m2p4'))

# up to the diagonal, the following are equivalent:
cleanup( precmat.IGMRFreglat(3,4) -
             precmat.GMRFreglat(3,4,1, 'm1p1'))

Description

Displaying data over the administrative districts of Germany

Usage

map.landkreis(data, col=NULL, zlim=range(data), add=FALSE,
              legendpos=c( 0.88,0.9,0.05,0.4))

Arguments

Details

The function germany.plot super-seeds map.landkreis (it is several factors faster).

The perfect position of the legend is an art per se and depends on various par parameters. See also the source code of the function image.plot of fields.

Author(s)

Reinhard Furrer

References

The code of map.landkreis is very similar to germany.map from the package INLA.

See Also

germany.plot super-seeding map.landkreis.

Examples

## Not run: 
data( Oral)
par( mfcol=c(1,2))
germany.plot( log( Oral$Y), legend=TRUE)
map.landkreis( log( Oral$Y))

## End(Not run)

Description

Applies the Math group functions to spam objects

Usage

# ceiling(x)
# floor(x)
# exp(x, base = exp(1))
# log(x, base = exp(1))
# sqrt(x)

# abs(x)
# cumprod(x)
# cumsum(x)

# cos(x)
# sin(x)
# tan(x)

# acosh(x)

Arguments

Details

It is important to note that the zero entries do not enter the evaluation when structurebased=FALSE. The operations are performed on the stored non-zero elements. This may lead to differences if compared with the same operation on a full matrix.

Value

If the option spam.structurebased=TRUE, all functions operate on the vector x@entries and return the result thereof.
Conversely, if structurebased=FALSE, the result is identical to one with as.matrix(x) input and an as.spam purger.

Author(s)

Reinhard Furrer

See Also

Summary.spam, Ops.spam and Math2.spam

Examples

getGroupMembers("Math")

mat <- matrix(c( 1,2,0,3,0,0,0,4,5),3)
smat <- as.spam( mat)
cos( mat)
cos( smat)

options(spam.structurebased=FALSE)
cos( smat)

sqrt( smat)

Description

Applies the Math2 group functions to 'spam' objects

Usage

## S4 method for signature 'spam'
round(x, digits = 0)
## S4 method for signature 'spam'
signif(x, digits = 6)

Arguments

Value

All functions operate on the vector x@entries and return the result thereof.

Author(s)

Reinhard Furrer

See Also

Ops.spam and Math.spam

Examples

getGroupMembers("Math2")

set.seed(12)
smat <- diag.spam( rnorm(15))
round(smat, 3)

Description

Maximum likelihood estimates of a simple spatial model

Usage

neg2loglikelihood.spam(y, X, distmat, Covariance,
                 beta, theta, Rstruct = NULL, cov.args = NULL, ...)

neg2loglikelihood(y, X, distmat, Covariance,
                 beta, theta, cov.args = NULL, ...)

neg2loglikelihood.nomean(y, distmat, Covariance,
                 theta, cov.args = NULL, ...)

mle.spam(y, X, distmat, Covariance,
     beta0, theta0, thetalower, thetaupper, optim.control=NULL,
     Rstruct = NULL, hessian = FALSE, cov.args = NULL, ...)

mle(y, X, distmat, Covariance,
     beta0, theta0, thetalower, thetaupper, optim.control=NULL,
     hessian = FALSE, cov.args = NULL, ...)

mle.nomean.spam(y, distmat, Covariance,
     theta0, thetalower, thetaupper, optim.control=NULL,
     Rstruct = NULL, hessian = FALSE, cov.args = NULL, ...) 

mle.nomean(y, distmat, Covariance,
     theta0, thetalower, thetaupper, optim.control=NULL,
     hessian = FALSE, cov.args = NULL, ...)

Arguments

Details

We provide functions to calculate the negative-2-log-likelihood and maximum likelihood estimates for the model

y ~ N_n( X beta, Sigma(h;theta) )

in the case of a sparse or ordinary covariance matrices.

In the case of the *.spam versions, the covariance function has to return a spam object. In the other case, the methods are correctly overloaded and work either way, slightly slower than the *.spam counterparts though.

When working on the sphere, the distance matrix has to be transformed by

h -> R / 2 sin(h/2)

where R is the radius of the sphere.

The covariance function requires that the first argument is the distance matrix and the second the parameters. One can image cases in which the covariance function does not take the entire distance matrix but only some partial information thereof. (An example is the use of a kronecker type covariance structure.) In case of a sparse covariance construction where the argument Rstruct is not given, the first parameter element needs to be the range parameter. (This results from the fact, that a sparse structure is constructed that is independent of the parameter values to exploit the fast Choleski decomposition.)

In the zero-mean case, the neg2loglikelihood is calculated by setting the parameters X or beta to zero.

Value

The negative-2-loglikelihood or the output from the function optim.

Author(s)

Reinhard Furrer

See Also

covmat, rmvnorm.spam

Examples

# True parameter values:
truebeta <- c(1,2,.2)    # beta = (intercept, linear in x, linear in y)
truetheta <- c(.5,2,.02) # theta = (range, sill, nugget)



# We now define a grid, distance matrix, and a sample:
x <- seq(0,1,l=5)
locs <- expand.grid( x, x)
X <- as.matrix( cbind(1,locs))  # design matrix

distmat <- nearest.dist( locs, upper=NULL) # distance matrix
Sigma <- cov.sph( distmat, truetheta)    # true covariance matrix


set.seed(15)
y <- c(rmvnorm.spam(1,X %*% truebeta,Sigma)) # construct sample

# Here is the negative 2 log likelihood:
neg2loglikelihood.spam( y, X, distmat, cov.sph,
                       truebeta, truetheta)

# We pass now to the mle:
res <- mle.spam(y, X, distmat, cov.sph,
         truebeta, truetheta,thetalower=c(0,0,0),thetaupper=c(1,Inf,Inf))

# Similar parameter estimates here, of course:
mle.nomean.spam(y-X%*%res$par[1:3], distmat, cov.sph,
         truetheta, thetalower=c(0,0,0), thetaupper=c(1,Inf,Inf))

Description

This function computes and returns specific elements of distance matrix computed by using the specified distance measure.

Usage

nearest.dist( x, y=NULL, method = "euclidean",
             delta = 1, upper = if (is.null(y)) FALSE else NULL,
             p = 2, miles = TRUE, R = NULL, fortran = FALSE)

Arguments

Details

For great circle distance, the matrices x and y contain the degrees longitudes in the first and the degrees latitudes in the second column. delta is in degrees. Hence to restrict to distances smaller than delta.km, one has to specify delta=delta.km*360/(6378.388*2*pi).

The distances are calculated based on spherical law of cosines. Care is needed for ‘zero’ distances due to the final acosin: acos(1-1e-16), especially with an actual radius.

Default value of Earth's radius is 3963.34miles (6378.388km).

There are many other packages providing distance functions. Especially for great circle distances there are considerable differences between the implementations. For high precision results, sp::spDists is a good candidate and distances of large amount of locations can be processed in parallel with the parallelDist package.

The formerly depreciated arguments eps and diag are now eliminated.

x and y can be any object with an existing as.matrix method.

The Fortran code is based on a idea of Doug Nychka.

Value

A spam object containing the distances spanned between zero and delta. The sparse matrix may contain many zeros (e.g., collocated data). However, to calculate covariances, these zeros are essential.

Author(s)

Annina Cincera (C++ code), Reinhard Furrer

See Also

spam_rdist

Examples

# Note that upper=T and using t(X)+X is quicker than upper=NULL;
#     upper=T marginally slower than upper=F.

# To compare nearest.dist with dist, use as.dist(...)
nx <- 4
x <- expand.grid(as.double(1:nx),as.double(1:nx))
sum( ( as.dist(nearest.dist( x, delta=nx*2))-
          dist(x)                            )^2)

# Create nearest neighbor structures:
par(mfcol=c(1,2))
x <- expand.grid(1:nx,1:(2*nx))
display( nearest.dist( x, delta=1))
x <- expand.grid(1:(2*nx),1:nx)
display( nearest.dist( x, delta=1))

Description

Allow the user to set and examine a variety of options which affect the way in which R computes and displays sparse matrix results.

Details

Invoking options() with no arguments returns a list with the current values of the options. To access the value of a single option, one should use getOption("spam.eps"), e.g., rather than options("spam.eps") which is a list of length one.

Of course, printing is still subordinate to getOption("max.print") or similar options.

Value

For getOption, the current value set for option x, or NULL if the option is unset.

For options(), a list of all set options sorted by category. For options(name), a list of length one containing the set value, or NULL if it is unset. For uses setting one or more options, a list with the previous values of the options changed (returned invisibly).

Options used for the package spam

A short description with the default values follows.

spam.eps=.Machine$double.eps:

values smaller than this are considered as zero. This is only used when creating spam objects.

spam.drop=FALSE:

default parameter for drop when subsetting

spam.printsize=100:

the max number of elements of a matrix which we display as regular matrix.

spam.imagesize=10000:

the max number of elements of a matrix we display as regular matrix with image or display. Larger matrices are represented as dots only.

spam.cex=1200:

default dot size for image or display.

spam.structurebased=FALSE:

should operations be carried out on the nonzero entries (the structure) or including the zeros.

spam.inefficiencywarning=1e6:

issue a warning when inefficient operations are performed and the matrix exceeds the specified size. Valid value is a postive integer or a logical. TRUE corresponds to 1 (always), FALSE to Inf.

spam.trivalues=FALSE:

a flag whether to return the structure (FALSE) or the values themselves (TRUE) when returning the upper and lower triangular part of a matrix.

spam.listmethod="PE":

algorithm for spam.list. Default is suggestion by Paul Eilers (thanks). Any other specification uses a bubble sort algorithm which is only slightly faster for very sparse matrices.

spam.dopivoting=TRUE:

default parameter for "solve" routines. FALSE would solve the system without using the permutation.

spam.NAOK=FALSE:

logical determines if NA, NaN and Inf are allowed to Fortan. Setting to TRUE allows to work with these but full functionality has not been tested.

spam.safemodevalidity=TRUE:

logical determines if sanity check is peformed when constructing sparse matrices. Default is safer but somewhat slower.

spam.cholsymmetrycheck=TRUE:

for the Cholesky factorization, verify if the matrix is symmetric.

spam.cholpivotcheck=TRUE:

for the Cholesky factorization, when passing a permutation, should a minimum set of checks be performed?

spam.cholupdatesingular="warning":

for a Cholesky update, what happens if the matrix is singular: "warning" only and returning the not updated factor, "error" or return simply "NULL".

spam.cholincreasefactor=c(1.25,1.25):

If not enought memory could be allocated, these are the steps to increase it.

spam.nnznearestdistnnz=c(400^2,400):

Memory allocation parameters for nearest.dist.

spam.nearestdistincreasefactor=1.25:

If not enought memory could be allocated, this is the step to increase it.

See Also

Functions influenced by these options include: print.spam, display.spam, image.spam, upper.tri.spam, chol.spam, nearest.dist, etc.
powerboost

Examples

smat <- diag.spam( 1:8)
smat
options(spam.printsize=49)
smat

# List all spam options:
options()[grep("spam",names(options()))]

# Reset to default values:
options(spam.eps=.Machine$double.eps,
        spam.drop=FALSE,
        spam.printsize=100,
        spam.imagesize=10000,
        spam.cex=1200,
        spam.structurebased=FALSE,
        spam.inefficiencywarning=1e6,
        spam.trivalues=FALSE,
        spam.listmethod="PE",
        spam.NAOK=FALSE,
        spam.safemodevalidity=TRUE,
        spam.dopivoting=TRUE,
        spam.cholsymmetrycheck=TRUE,
        spam.cholpivotcheck=TRUE,
        spam.cholupdatesingular="warning",
        spam.cholincreasefactor=c(1.25,1.25),
        spam.nearestdistincreasefactor=1.25,
        spam.nearestdistnnz=c(400^2,400))

Description

Oral cavity cancer counts in 544 districts in Germany over 1986-1990.

Format

Oral is a dataframe with 3 columns.

Y

observed counts

E

expected counts

SMR

standardized mortality ratios

germany is a list of 544 elements, each describing an individual polygon of the district.

Details

The expected counts depend on the number of people in the region and their age distribution.
The regions are ordered according the supplied polygon description and adjacency graph.

There is a similar dataset data(Germany) with larynx cancer cases from the package INLA as well, with an additional smoking covariate.

Source

The data is available from the package INLA distributed from https://www.r-inla.org.

References

Knorr-Held, L. and Rasser, G. (2000) Bayesian Detection of Clusters and Discontinuities in Disease Maps, Biometrics, 56, 13–21.

See Also

germany.plot.

Description

Extract the (inverse) permutation used by the Cholesky decomposition

Usage

ordering( x, inv=FALSE)

Arguments

Details

Recall that calculating a Cholesky factor from a sparse matrix consists of finding a permutation first, then calculating the factors of the permuted matrix. The ordering is important when working with the factors themselves.

The ordering from a full/regular matrix is 1:n.

Note that there exists many different algorithms to find orderings.

See the examples, they speak more than 10 lines.

Author(s)

Reinhard Furrer

See Also

chol.spam, solve.spam.

Examples

# Construct a pd matrix S to work with (size n)
n <- 100    # dimension
S <- .25^abs(outer(1:n,1:n,"-"))
S <- as.spam( S, eps=1e-4)
I <- diag(n)  # Identity matrix

cholS <- chol( S)
ord <- ordering(cholS)
iord <- ordering(cholS, inv=TRUE)

R <- as.spam( cholS ) # R'R = P S P', with P=I[ord,],
  # a permutation matrix (rows permuted).
RtR <- t(R) %*% R

# the following are equivalent:
as.spam( RtR -            S[ord,ord],    eps=1e-15)
as.spam( RtR[iord,iord] - S,             eps=1e-15)
as.spam( t(R[,iord]) %*% R[,iord] - S, eps=1e-15)
# we use 'eps' to avoid issues close to machine precision

# trivially:
as.spam( t(I[iord,]) - I[ord,])  # (P^-1)' = P  
as.spam( t(I[ord,]) - I[,ord])  # 
as.spam( I[iord,] - I[,ord])
as.spam( I[ord,]%*%S%*%I[,ord] - S[ord,ord] )
   # pre and post multiplication with P and P' is ordering

Description

Resets the dimension of a spam matrix by truncation or zero padding.

Usage

pad(x) <- value

Arguments

Details

It is important to notice the different behavior of the replacement method for ordinary arrays and spam objects (see ‘Examples’). Here, the elements are not simply rearranged but an entirely new matrix is constructed. If the new column dimension is smaller than the original, the matrix is also cleaned (with option("spam.eps") as filter).

Value

A (spam) matrix of dimension value where trunction or padding has been used.

Author(s)

Reinhard Furrer

See Also

dim.spam.

Examples

x <- diag(4)
dim(x)<-c(2,8)
x

s <- diag.spam(4)
pad(s) <- c(7,3)  # any positive value can be used

s <- diag.spam(4)
pad(s) <- c(2,8)  # result is different than x

Description

Row and/or column permutes a (sparse) matrix.

Usage

permutation.spam(A, P=NULL, Q=NULL, ind=FALSE, check=TRUE)

Arguments

Details

If P and Q are permutation matrices, the result is PAQ. However, it is also possible to specify the indices and to perform in a very efficient way A[rowind, colind], see examples.

A row permutation is much faster than a colum permutation. For very large matrices, a double transpose might be faster.

The spam option spam.checkpivot determines if the permutation is verified.

Value

A permuted matrix.

Author(s)

Reinhard Furrer

See Also

ordering, spam.options.

Examples

A <- spam(1:12,3)
P <- c(3,1,2)
Q <- c(2,3,1,4)

permutation(A,P,Q)-A[order(P),order(Q)]

permutation(A,P,Q,ind=TRUE)-A[P,Q]

Description

Sets several options for speed-up.

Usage

powerboost(flag)

Arguments

Details

The options turn checking off ("safemode", "cholsymmetrycheck" and "cholpivotcheck") and switch to single precision for "eps".

Value

NULL in any case.

Author(s)

Reinhard Furrer, after receiving too much C.mc.st adds.

See Also

spam.options.

Description

Fast ways to create sparse precision matrices for various IGMRF.

Usage

precmat(n, season=12, m=n, A=NULL, order=1, ... , type="RW1")
precmat.RW1(n)
precmat.RW2(n)
precmat.RWn(n, order=3)
precmat.season(n, season=12)
precmat.IGMRFreglat(n, m, order=1, anisotropy=1)
precmat.IGMRFirreglat(A, eps=getOption("spam.eps"))

Arguments

Details

precmat is a wrapper that calls the other functions according to the argument type.
Implements many of the precision matrices discussed in Chapter 3 of Rue and Held (2005). For example, precmat.RW1, precmat.RW2 and precmat.season are given in equations (3.22), (3.40) and (3.59); precmat.IGMRFreglat on page 107. Note that for the latter we reverse the order of the dimension here!
If adjacency matrix is a regular matrix, it is coerced to a spam object. Only the structure is used. Make sure, that the diagonal is empty.

Value

A sparse precision matrix.

Author(s)

Reinhard Furrer

References

Rue and Held (2005).

See Also

precmat.GMRFreglat, rmvnorm.prec, adjacency.landkreis.

Examples

n <- 10
Q <- precmat.RW2( n)
# rmvnorm.prec(1, Q=Q) # does not work, because the matrix is singular.

Q%*%cbind(1,1:n)

Description

Printing (non-zero elements) of sparse matrices and summarizing the sparsity structure thereof.

Usage

## S4 method for signature 'spam'
print(x, ...)
## S4 method for signature 'spam'
summary(object, ...)

Arguments

Details

getOption('spam.printsize') determines if the sparse matrix is coerced into a matrix and the printed as an array or if only the non-zero elements of the matrix are given.

Value

NULL for print, because the information is printed with cat there is no real need to pass any object back.
A list containing the non-zero elements and the density for summary for class spam.
A list containing the non-zero elements of the factor, the density and the fill-in for summary for class spam.chol.NgPeyton.

Author(s)

Reinhard Furrer

See Also

display or image for a graphical visualization; spam.options

Examples

set.seed(13)
smat <- spam_random(8)
par(mfcol=c(1,2),pty='s')
options(spam.printsize=1000)
print(smat)
options(spam.printsize=10)
print(smat)
summary(smat)
summary(smat)$nnz

smat@entries[1:5] <- 0
print(smat, minimal = FALSE)
print(smat, minimal = FALSE, rowpointer = TRUE)
smat@rowpointers
print_nnzpos(smat)

Description

Creates random spam matrix given the dimension and other parameters.

Usage

spam_random(nrow = 1L, ncol = nrow, density = 0.5, distribution = NULL, digits = NULL,
            sym = FALSE, spd = FALSE, verbose = FALSE, ...)

Arguments

Details

To create a random spam64 matrix, set options(spam.force64 = TRUE).

Value

A random matrix in spam format.

Author(s)

Florian Gerber, Roman Flury, Reinhard Furrer

See Also

spam-class and display.spam

Examples

set.seed(42)
rspam <- spam_random(500, digits = 2, distribution = rnorm, sd = 2, mean = 10, density = .01)
display.spam(rspam, cex = 2)

Description

Fast and intuitive ways to draw from a Gaussian random field.

Usage

rgrf( n,
   locs, nx, ny=nx, xlim=c(0,1), ylim=c(0,1), tau=0,
   Covariance, theta, beta=0, X, 
   method=c('chol'),  method.args=list(sparse=FALSE), 
   eps = getOption("spam.eps"), drop=TRUE, attributes=TRUE, ...)

Arguments

Details

If no locations are given, the function constructs these according a regular or a regular perturbed grid. The perturbation is determined by tau, which has to be greater than zero (no perturbation) and strictly smaller than 1/2 (max perturbation).

The regular grid has spacing (here for x) dx=diff(xlim)/nx and runs from xlim[1]+dx/2 to xlim[2]-dx/2. The locations are at least (1/nx-2*tau*dx) separated.

Currently, the only method implemented is a Cholesky factorization routine, (much as in rmvnorm).

The rdist() from the fields package is awefully fast. Unless one has very sparse covariance matrices, a sparse approach is not bringing a lot of improvements.

The methods may use different covariance construction approaches and thus the nesting of cov.args in method.args.

Author(s)

Reinhard Furrer

See Also

rgrf, chol and ordering.

Examples

require(fields)
# Regular grid with constant mean:
nx <- 10
field <- rgrf(1, nx=nx,  Covariance="cov.wend2", theta=c(.5, 1), beta=5)
quilt.plot(cbind(attr(field,"locs"),z=field), nx=nx, ny=nx)
points(attr(field,"locs"))

# Irregluar grid:
field <- rgrf(1, nx=10, tau=0.3, Covariance="cov.mat", theta=c(.2, 1, 1.5))
fields::quilt.plot(attr(field,"locs"), field)

Description

Fast ways to draw multivariate normals when the variance or precision matrix is sparse.

Usage

rmvnorm(n, mu=rep.int(0, dim(Sigma)[1]), Sigma, ..., mean, sigma)
rmvnorm.spam(n, mu=rep.int(0, dim(Sigma)[1]), Sigma, Rstruct=NULL, ..., mean, sigma)
rmvnorm.prec(n, mu=rep.int(0, dim(Q)[1]), Q, Rstruct=NULL, ...)
rmvnorm.canonical(n, b, Q, Rstruct=NULL, ...)

Arguments

Details

All functions rely on a Cholesky factorization of the covariance or precision matrix.
The functions rmvnorm.prec and rmvnorm.canonical do not require sparse precision matrices Depending on the the covariance matrix Sigma, rmvnorm or rmvnorm.spam is used. If wrongly specified, dispatching to the other function is done.
Default mean is zero. Side note: mean is added via sweep() and no gain is accieved by distinguishing this case.
Often (e.g., in a Gibbs sampler setting), the sparsity structure of the covariance/precision does not change. In such setting, the Cholesky factor can be passed via Rstruct in which only updates are performed (i.e., update.spam.chol.NgPeyton instead of a full chol).

Author(s)

Reinhard Furrer

References

See references in chol.

See Also

rgrf, chol and ordering.

Examples

# Generate multivariate from a covariance inverse:
# (usefull for GRMF)
set.seed(13)
n <- 25    # dimension
N <- 1000  # sample size
Sigmainv <- .25^abs(outer(1:n,1:n,"-"))
Sigmainv <- as.spam( Sigmainv, eps=1e-4)


Sigma <- solve( Sigmainv)  # for verification
iidsample <- array(rnorm(N*n),c(n,N))

mvsample <- backsolve( chol(Sigmainv), iidsample)
norm( var(t(mvsample)) - Sigma, type="m")

# compare with:
mvsample <- backsolve( chol(as.matrix( Sigmainv)), iidsample, n)
   #### ,n as patch
norm( var(t(mvsample)) - Sigma, type="m")


# 'solve' step by step:
b <- rnorm( n)
R <- chol(Sigmainv)
norm( backsolve( R, forwardsolve( R, b))-
      solve( Sigmainv, b) )
norm( backsolve( R, forwardsolve( R, diag(n)))- Sigma )

Description

Fast way to draw conditional multivariate normals when the covariance matrix is sparse.

Usage

rmvnorm.conditional(n, y, mu = rep.int(0, dim(SigmaXX)[1]+dim(SigmaYY)[1]),
                    SigmaXX, SigmaYY, SigmaXY, noise, RstructYY = NULL, ...)

Arguments

Details

Quite often, we want to draw condional observations $X|y$ from the model $Y=X+e$, where $X$ has covariance matrix SigmaXX and $e$ has white noise.

Covariance of $Y$ can be specified by SigmaYY or SigmaXX+diag(noise,). If $Y$ and $X$ do not have the same dimensions, SigmaXY needs to be specified.

The function also implmements a general multivariate model, where the we only observe part of the vector. The components are first $X$ then $Y$.

The function rmvnorm.cond() is a wrapper to rmvnorm.conditional() and included to increase similarities with other packages.

Author(s)

Reinhard Furrer

See Also

rmvnorm.spam.

Examples

set.seed(12)
N <- 300
y <- c(5, -5, -5, 5)
SigmaXX <- as.spam(.95^abs(outer(1:N, 1:N, "-")), eps=1e-4)
sel <- c(10, 100, 120, 300)        # where we observe y
SigmaXY <- SigmaXX[, sel]
SigmaYY <- SigmaXX[sel,sel] + diag.spam(.01, length(y)) # some noise
x <- rmvnorm.conditional(3, y, SigmaXX=SigmaXX, SigmaXY=SigmaXY,
                         SigmaYY=SigmaYY)
# unconditional sample:
ux <- rmvnorm(1, Sigma=SigmaXX)
matplot(t(rbind(x, ux)), type='l', lty=1)
points(sel, y, pch=19)

Description

Fast ways to draw multivariate normals with linear constrains when the variance or precision matrix is sparse.

Usage

rmvnorm.const(n, mu = rep.int(0, dim(Sigma)[1]), Sigma, Rstruct = NULL,
              A = array(1, c(1,dim(Sigma)[1])), a=0, U=NULL,  ...)
rmvnorm.prec.const(n, mu = rep.int(0, dim(Q)[1]), Q, Rstruct = NULL,
              A = array(1, c(1,dim(Q)[1])), a=0, U=NULL,  ...)
rmvnorm.canonical.const(n, b, Q, Rstruct = NULL,
              A = array(1, c(1,dim(Q)[1])), a=0, U=NULL, ...)