widals.snow {widals}

Fit WIDALS

Description

Locate the WIDALS hyperparameters

Usage

widals.snow(j, rm.ndx, Z, Hs, Ht, Hst.ls, locs, lags, b.lag, cv = 0, 
geodesic = FALSE, wrap.around = NULL, GP.mx, stnd.d = FALSE, ltco = -16)

Arguments

j

Index used by snowfall. A scalar integer. Which row of GP.mx to use for the ALS hyperparameters, GP.

rm.ndx

A list of vectors of indices to remove for k-fold cross-validation.

Z

Data. A \tau x n numeric matrix.

Hs

Spacial covariates. An n x p_s numeric matrix.

Ht

Temporal covariates. A \tau x p_t numeric matrix.

Hst.ls

Space-time covariates. A list of length \tau, each element containing a n x p_st numeric matrix.

locs

Locations of supporting sites. An n x 2 numeric matrix, first column is spacial x, second column is spacial y. If the geodesic is TRUE, make sure latitude is in the first column.

lags

Temporal lags for stochastic smoothing. An integer vector or scalar. E.g., if the data's time increment is daily, then lags = c(-1,0,1) would tell the enclosed function crispify smooth today's predictions using yesterdays, today's, and tomorrow's observed residuals.

b.lag

ALS lag. A scalar integer, typically -1 (a-prior), or 0 (a-posteriori).

cv

Cross-validation switch. Currently takes on a value of -2 or 2. See Details below.

geodesic

Use geodesic distance? Boolean. If true, distance (used internally) is in units kilometers.

wrap.around

**Unused.

GP.mx

Hyperparameters. A k.glob x 2 non-negative matrix. See MSS.snow.

stnd.d

Spacial compression. Boolean.

ltco

Weight threshold. A scalar number. A value of, e.g., -10, will instruct crispify to ignore weights less than 10^(-10) when smoothing.

Details

When the cv is set to 2, then this function uses spacial k-fold validation, according to the site indices present in rm.ndx. When cv is set to -2, self-referencing sites are given zero-weight, i.e., a site's value is not allowed to contribute to its predicted value.

Value

A \tau x n matrix. The WIDALS predictions at locs.

See Also

crispify, H.als.b, widals.predict.

Examples

	
set.seed(99999)

library(SSsimple)

tau <- 100
n.all <- 35

Hs.all <- matrix(rnorm(n.all), nrow=n.all)
Ht <- matrix(rnorm(tau*2), nrow=tau)
Hst.ls.all <- list()
for(i in 1:tau) { Hst.ls.all[[i]] <- matrix(rnorm(n.all*2), nrow=n.all) }

Hst.combined <- list()
for(i in 1:tau) { 
    Hst.combined[[i]] <- cbind( Hs.all, matrix(Ht[i, ], nrow=n.all, ncol=ncol(Ht), 
    byrow=TRUE), Hst.ls.all[[i]] ) 
}

locs.all <- cbind(runif(n.all, -1, 1), runif(n.all, -1, 1))
D.mx.all <- distance(locs.all, locs.all, FALSE)
R.all <- exp(-2*D.mx.all) + diag(0.01, n.all)

######## use SSsimple to simulate
sssim.obj <- SS.sim.tv( 0.999, Hst.combined, 0.01, R.all, tau )


n <- n.all
locs <- locs.all

Z.all <- sssim.obj$Z
Z <- Z.all


Hst.ls <- Hst.ls.all
Hs <- Hs.all

test.rng <- 20:tau

################  WIDALS, true cross-validation

rm.ndx <- create.rm.ndx.ls(n, 10)

cv <- 2
lags <- c(0)
b.lag <- 0

GP <- c(1/8, 1/12, 5, 0, 1)
GP.mx <- matrix(GP, ncol=length(GP))
Zwid <- widals.snow(j=1, rm.ndx, Z, Hs, Ht, Hst.ls, locs, lags, b.lag, cv = cv, 
geodesic = FALSE, wrap.around = NULL, GP.mx, stnd.d = FALSE, ltco = -16) 

errs.sq <- (Z - Zwid)^2
sqrt( mean(errs.sq[test.rng, ]) )


################  WIDALS, pseudo cross-validation

rm.ndx <- I(1:n)

cv <- -2
lags <- c(0)
b.lag <- -1

GP <- c(1/8, 1/12, 5, 0, 1)
GP.mx <- matrix(GP, ncol=length(GP))
Zwid <- widals.snow(j=1, rm.ndx, Z, Hs, Ht, Hst.ls, locs, lags, b.lag, cv = cv, 
geodesic = FALSE, wrap.around = NULL, GP.mx, stnd.d = FALSE, ltco = -16) 

errs.sq <- (Z - Zwid)^2
sqrt( mean(errs.sq[test.rng, ]) )