widals.snow: 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)

Value

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

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.

Examples

Run this code

	
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, ]) )