crispify: Observation-Space Stochastic Correction
In widals: Weighting by Inverse Distance with Adaptive Least Squares

Description Usage Arguments Details Value See Also Examples

Improve observation-space predictions using 'left over' spacial correlation between model residuals

1 2	crispify(locs1, locs2, Z.delta, z.lags.vec, geodesic, alpha, flatten, self.refs, lags, stnd.d = FALSE, log10cutoff = -16)

`locs1`	Locations of supporting sites. An n x 3 matrix, first column is spacial x, second column is spacial y, third column contains relative temporal 'distance'. If the `geodesic` is `TRUE`, make sure latitude is in the first column.
`locs2`	Locations of interpolation sites. An n* x 3 matrix, where n* is the number of interpolation sites. See `locs1` above.
`Z.delta`	Observed residuals. A τ x x matrix.
`z.lags.vec`	Temporal lags. An integer vector or scalar.
`geodesic`	Use geodesic distance? Boolean. If true, distance (used internally) is in units kilometers.
`alpha`	The WIDALS distance rate hyperparameter. A scalar non-negative number.
`flatten`	The WIDALS 'flattening' hyperparameter. A scalar non-negative number. Typically between 0 and some number slightly greater than 1. When 0, no crispification.
`self.refs`	Which sites are self-referencing? An integer vector of (zero-based) lag indices, OR a scalar set to `-1`. This argument only has meaning when `locs1` is identical to `locs2`. If the `lags` argument is, say, 0, then it would be pointless to smooth predictions with existing values. In this case, we can set `self.refs = 0`. If `locs1` is NOT the same as `locs2`, then set this argument to `-1`.
`lags`	Temporal lags. An integer vector or scalar. E.g., if the data's time increment is daily, then `lags = c(-1,0,1)` would have `crispify` smooth today's predictions using yesterdays, today's, and tomorrow's observed residuals.
`stnd.d`	Spacial compression. Boolean.
`log10cutoff`	Weight threshold. A scalar number. A value of, e.g., -10, will instruct `crispify` to ignore weights less than 10^(-10) when smoothing.

This function is called inside widals.predict and widals.snow. It may be useful for the user in building their own WIDALS model extensions.

A τ x x matrix.

widals.predict, widals.snow.

######### here's an itty-bitty example

######### simulate itty-bitty data

tau <- 21 #### number of time points

d.alpha <- 2
R.scale <- 1
sigma2 <- 0.01
F <- 1
Q <- 0

n.all <- 14 ##### number of spacial locations

set.seed(9999)


library(SSsimple)

locs.all <- cbind(runif(n.all, -1, 1), runif(n.all, -1, 1)) #### random location of sensors
D.mx <- distance(locs.all, locs.all, FALSE) #### distance matrix

#### create measurement variance using distance and covariogram
R.all <- exp(-d.alpha*D.mx) + diag(sigma2, n.all)

Hs.all <- matrix(1, n.all, 1) #### constant mean function

##### use SSsimple to simulate system
xsssim <- SS.sim(F=F, H=Hs.all, Q=Q, R=R.all, length.out=tau, beta0=0)
Z.all <- xsssim$Z ###### system observation matrix


######## suppose use the global mean as a prediction

z.mean <- mean(Z.all)

Z.delta <- Z.all - z.mean


z.lags.vec <- rep(0, n.all)

geodesic <- FALSE
alpha <- 5
flatten <- 1

## emmulate cross-validation, i.e., 
## don't use observed site values to predict themselves (zero-based)
self.refs <- 0 
lags <- 0

locs1 <- cbind(locs.all, rep(0, n.all))
locs2 <- cbind(locs.all, rep(0, n.all))

Z.adj <- crispify(locs1, locs2, Z.delta, z.lags.vec, geodesic, alpha, 
    flatten, self.refs, lags, stnd.d = FALSE, log10cutoff = -16) 

Z.adj

Z.hat <- z.mean + Z.adj

sqrt( mean( (Z.all - Z.hat)^2 ) )


######### set flatten to zero -- this means no crispification

Z.adj <- crispify(locs1, locs2, Z.delta, z.lags.vec, geodesic, alpha, 
    flatten=0, self.refs, lags, stnd.d = FALSE, log10cutoff = -16) 

Z.adj

Z.hat <- z.mean + Z.adj

sqrt( mean( (Z.all - Z.hat)^2 ) )