MSS.snow: Metaheuristic Stochastic Search

Description

Locate WIDALS hyperparameters

Usage

MSS.snow(FUN.source, current.best, p.ndx.ls, f.d, sds.mx, k.glob, k.loc.coef, X = NULL)

Value

Nothing. After completion, the best hyperparameters, GP, are assigned to the Global Environment.

Arguments

FUN.source: Search function definitions (see Details). A path to source code, or function, e.g., fun.load.widals.a.
current.best: An initial cost. A scalar. Setting to NA will cause MSS.snow to make an initial pass over the data to create an initial cost to beat.
p.ndx.ls: Hyperparameter indices (of GP) to search. A list of vectors. For example, list( c(1,2), c(3,4,5) ) will instruct MSS.snow, for each local search, to search over the first two hyperparameters as a pair, then to search the last three as a group.
f.d: Local search functions. A list of functions (one for each element of GP). Typically, for WIDALS, all five will be dlog.norm.
sds.mx: The standard deviations for f.d. An k.glob x q matrix, where q is the number of hyperparameters, i.e., the length of GP.
k.glob: The number of global searches. A scalar integer.
k.loc.coef: The coeficient for the number of local searches to make. A scalar integer.
X: A placeholder for values to be passed between functions inside MSS.snow (see Details).

Details

This function requires the presence of a number of values and functions out-of-scope. It is assumed that these are available in the Global Environment. They are: run.parallel (boolean), FUN.MH (a function that creates, for a given GP, a cost), FUN.GP (a function that applies constraints to GP), FUN.I (a function that does something when local searches have reduced the cost), FUN.EXIT (a function that does something when MSS.snow is done).

Examine the code for fun.load.widals.a for an example of the four functions described above. Note that these four functions may themselves require objects out-of-scope.

In general, for a given R session, special care should be taken concerning the naming and assigning of the following objects: Z (the space-time data), Z.na (a boolean matrix indicating missing values in Z), locs (site locations), Hs (spacial covariates), Ht (temporal covariates), Hst.ls (space-time covariates), lags (temporal lag vector), b.lag (the ALS lag), cv (cross-validation switch), xgeodesic (boolean), ltco (weight cut-off), GP (hyperparameter vector), run.parralel (boolean), stnd.d (boolean), train.rng (time index vector), test.rng (time index vector).

Examples

Run this code


##### simulate a state-space system (using pkg SSsimple)

if (FALSE) {

### using dontrun because of excessive run time for CRAN submission

set.seed(9999)

library(SSsimple)


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

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

udom <- (0:300)/100
plot( udom,    R.scale * exp(-d.alpha*udom) ,  type="l", col="red" ) #### see the covariogram

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

set.seed(9999)
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

	
	
########### now make assignments required by MSS.snow
	


##### randomly remove five sites to serve as interpolation points
ndx.interp <- sample(1:n.all, size=5) 
ndx.support <- I(1:n.all)[ -ndx.interp ] ##### support sites



########### what follows are important assignments, 
########### since MSS.snow and the four helper functions
########### will look for these in the Global Environment 
########### to commence fitting the model (as noted in Details above)
train.rng <- 30:(tau) ; test.rng <- train.rng

Z <- Z.all[ , ndx.support ] 
Hs <- Hs.all[ ndx.support, , drop=FALSE] 
locs <- locs.all[ndx.support, , drop=FALSE] 

Ht <- NULL
Hst.ls <- NULL

lags <- c(0) 
b.lag <- c(-1) 
cv <- -2
xgeodesic <- FALSE
stnd.d <- FALSE
ltco <- -10
GP <- c(1/10, 1, 20, 20, 1) ### -- initial hyperparameter values
run.parallel <- TRUE 

if( cv==2 ) { rm.ndx <- create.rm.ndx.ls( nrow(Hs), 14 ) } else { rm.ndx <- 1:nrow(Hs) }
rgr.lower.limit <- 10^(-7) ; d.alpha.lower.limit <- 10^(-3) ; rho.upper.limit <- 10^(4)


############## tell snowfall to use two threads for local searches
sfInit(TRUE, cpus=2)
fun.load.widals.a()


######## now, finally, search for best fit over support
######## Note that p.ndx.ls and f.d are produced inside fun.load.widals.a()
MSS.snow(fun.load.widals.a, NA, p.ndx.ls, f.d, matrix(1/10, 10, length(GP)), 10, 7)
sfStop()

######## we can use these hyperparameters to interpolate to the 
######## deliberately removed sites, and measure MSE, RMSE
Z0.hat <- widals.predict(Z, Hs, Ht, Hst.ls, locs, lags, b.lag, 
Hs0=Hs.all[ ndx.interp, , drop=FALSE ], 
Hst0.ls=NULL, locs0=locs.all[ ndx.interp, , drop=FALSE],
geodesic = xgeodesic, wrap.around = NULL, GP, stnd.d = stnd.d, ltco = ltco)

resids.wid <- ( Z.all[ , ndx.interp ] - Z0.hat )
mse.wid <- mean( resids.wid[ test.rng, ]^2 )
mse.wid
sqrt(mse.wid)






########################################### Simulated Imputation with WIDALS
Z.all <- xsssim$Z
Z.missing <- Z.all

Z.na.all <- matrix( sample(c(TRUE, FALSE), size=n.all*tau, prob=c(0.01, 0.99), replace=TRUE), 
tau, n.all)
Z.missing[ Z.na.all ] <- NA


Z <- Z.missing
Z[ is.na(Z) ] <- mean(Z, na.rm=TRUE)
X <- list("Z.fill"=Z)

Z.na <- Z.na.all
Hs <- Hs.all
locs <- locs.all
Ht <- NULL
Hst.ls <- NULL
lags <- c(0)
b.lag <- c(-1)
cv <- -2
xgeodesic <- FALSE
ltco <- -10
if( cv==2 ) { rm.ndx <- create.rm.ndx.ls( nrow(Hs), 14 ) } else { rm.ndx <- 1:nrow(Hs) }

GP <- c(1/10, 1, 20, 20, 1)

rgr.lower.limit <- 10^(-7) ; d.alpha.lower.limit <- 10^(-3) ; rho.upper.limit <- 10^(4)

run.parallel <- TRUE

sfInit(TRUE, cpus=2)
fun.load.widals.fill()

MSS.snow(fun.load.widals.fill, NA, p.ndx.ls, f.d, 
seq(2, 0.01, length=10)*matrix(1/10, 10, length(GP)), 10, 7, X=X)
sfStop()

sqrt(mean(( (Z.all[train.rng, ] - Z.fill[train.rng, ])^2 )[ Z.na[ train.rng, ] ]))



    
    
    
    

############################################ Now Try with ALS alone

Z.all <- xsssim$Z

GP <- c(1/10, 1) ### -- initial hyperparameter values

############## tell snowfall to use two threads for local searches
sfInit(TRUE, cpus=2)
fun.load.hals.a()

######## now, finally, search for best fit over support
######## Note that p.ndx.ls and f.d are produced inside fun.load.widals.a()
MSS.snow(fun.load.hals.a, NA, p.ndx.ls, f.d, matrix(1/10, 10, length(GP)), 10, 7)
sfStop()

######## we can use these hyperparameters to interpolate to the deliberately removed sites, 
######## and measure MSE, RMSE
hals.obj <- H.als.b(Z, Hs, Ht, Hst.ls, rho=GP[1], reg=GP[2], b.lag = b.lag, 
Hs0 = Hs.all[ ndx.interp, , drop=FALSE ], Ht0 = NULL, Hst0.ls = NULL)
Z0.hat <- hals.obj$Z0.hat

resids.als <- ( Z.all[ , ndx.interp ] - Z0.hat )
mse.als <- mean( resids.als[ test.rng, ]^2 )
mse.als
sqrt(mse.als)



########################################### Simulated Imputation with ALS
Z.all <- xsssim$Z
Z.missing <- Z.all

set.seed(99)
Z.na.all <- matrix( sample(c(TRUE, FALSE), size=n.all*tau, prob=c(0.03, 0.97), replace=TRUE), 
tau, n.all)
Z.missing[ Z.na.all ] <- NA


Z <- Z.missing
Z[ is.na(Z) ] <- 0 #mean(Z, na.rm=TRUE)
X <- list("Z.fill"=Z)
    
Z.na <- Z.na.all

Hs <- Hs.all

GP <- c(1/10, 1) ### -- initial hyperparameter values

sfInit(TRUE, cpus=2)
fun.load.hals.fill()

MSS.snow(fun.load.hals.fill, NA, p.ndx.ls, f.d, 
seq(3, 0.01, length=10)*matrix(1, 10, length(GP)), 10, 7, X=X)
sfStop()

sqrt(mean(( (Z.all[train.rng, ] - Z.fill[train.rng, ])^2 )[ Z.na[ train.rng, ] ]))

}