Description Usage Arguments Details Value See Also Examples
Locate WIDALS hyperparameters
1 |
FUN.source |
Search function definitions (see Details). A path to source code, or function, e.g., |
current.best |
An initial cost. A scalar. Setting to NA will cause |
p.ndx.ls |
Hyperparameter indices (of |
f.d |
Local search functions. A list of functions (one for each element of GP). Typically, for WIDALS, all five will be |
sds.mx |
The standard deviations for |
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 |
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).
Nothing. After completion, the best hyperparameters, GP
, are assigned to the Global Environment.
Hals.fastcv.snow
, Hals.snow
, widals.snow
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 | ##### simulate a state-space system (using pkg SSsimple)
## Not run:
### 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, ] ]))
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.