Extract specific Matern parameters from a parameter sequence

This function is used to obtain specific Matern parameters (e.g., range or smoothness) from the covparams slot of a PrestoGPModel object.

Usage

create_param_sequence(P, ns = 1)

Arguments

P: Number of outcome variables
ns: Number of scale parameters

Value

A matrix with five rows and two columns as described below:

Row 1:: Starting and ending indices for the sigma parameter(s)
Row 2:: Starting and ending indices for the scale parameter(s)
Row 3:: Starting and ending indices for the smoothness parameter(s)
Row 4:: Starting and ending indices for the nugget(s)
Row 5:: Starting and ending indices for the correlation parameter(s)

Details

This function is intended for advanced users who want to specify the input Matern parameters for functions such as vecchia_Mlikelihood or createUMultivariate. To extract the Matern parameters from a fitted PrestoGP model, it is strongly recommended to use link{get_theta} instead.

References

Apanasovich, T.V., Genton, M.G. and Sun, Y. "A valid Matérn class of cross-covariance functions for multivariate random fields with any number of components", Journal of the American Statistical Association (2012) 107(497):180-193.
Genton, M.G. "Classes of kernels for machine learning: a statistics perspective", The Journal of Machine Learning Research (2001) 2:299-312.

Examples

# Space/elevation model
data(soil250, package="geoR")
y2 <- soil250[,7]               # predict pH level
X2 <- as.matrix(soil250[,c(4:6,8:22)])
# columns 1+2 are location coordinates; column 3 is elevation
locs2 <- as.matrix(soil250[,1:3])

soil.vm2 <- new("VecchiaModel", n_neighbors = 10)
# fit separate scale parameters for location and elevation
soil.vm2 <- prestogp_fit(soil.vm2, y2, X2, locs2, scaling = c(1, 1, 2))
#> 
#> Estimating initial beta... 
#> Estimation of initial beta complete 
#> 
#> Beginning iteration 1 
#> Estimating theta... 
#> Estimation of theta complete 
#> Estimating beta... 
#> Estimation of beta complete 
#> Iteration 1 complete 
#> Current penalized negative log likelihood: -254.4925 
#> Current MSE: 0.009739877 
#> Beginning iteration 2 
#> Estimating theta... 
#> Estimation of theta complete 
#> Estimating beta... 
#> Estimation of beta complete 
#> Iteration 2 complete 
#> Current penalized negative log likelihood: -254.4925 
#> Current MSE: 0.01017229 

pseq <- create_param_sequence(1, 2)
soil2.params <- soil.vm2@covparams
# sigma
soil2.params[pseq[1,1]:pseq[1,2]]
#> [1] 0.005850106
# scale parameters
soil2.params[pseq[2,1]:pseq[2,2]]
#> [1] 7.04081881 0.07812963
# smoothness parameter
soil2.params[pseq[3,1]:pseq[3,2]]
#> [1] 1.167065
# nugget
soil2.params[pseq[4,1]:pseq[4,2]]
#> [1] 0.004038009

# Multivariate model
ym <- list()
ym[[1]] <- soil250[,4] # predict sand/silt portion of the sample
ym[[2]] <- soil250[,5]
ym[[3]] <- soil250[,6]
Xm <- list()
Xm[[1]] <- Xm[[2]] <- Xm[[3]] <- as.matrix(soil250[,7:22])
locsm <- list()
locsm[[1]] <- locsm[[2]] <- locsm[[3]] <- as.matrix(soil250[,1:3])

soil.mvm <-  new("MultivariateVecchiaModel", n_neighbors = 10)
soil.mvm <- prestogp_fit(soil.mvm, ym, Xm, locsm)
#> 
#> Estimating initial beta... 
#> Estimation of initial beta complete 
#> 
#> Beginning iteration 1 
#> Estimating theta... 
#> Estimation of theta complete 
#> Estimating beta... 
#> Estimation of beta complete 
#> Iteration 1 complete 
#> Current penalized negative log likelihood: 1077.458 
#> Current MSE: 2.38111 
#> Beginning iteration 2 
#> Estimating theta... 
#> Estimation of theta complete 
#> Estimating beta... 
#> Estimation of beta complete 
#> Iteration 2 complete 
#> Current penalized negative log likelihood: 1026.008 
#> Current MSE: 3.738481 
#> Beginning iteration 3 
#> Estimating theta... 
#> Estimation of theta complete 
#> Estimating beta... 
#> Estimation of beta complete 
#> Iteration 3 complete 
#> Current penalized negative log likelihood: 1020.135 
#> Current MSE: 3.961958 
#> Beginning iteration 4 
#> Estimating theta... 
#> Estimation of theta complete 
#> Estimating beta... 
#> Estimation of beta complete 
#> Iteration 4 complete 
#> Current penalized negative log likelihood: 1018.501 
#> Current MSE: 4.09541 
#> Beginning iteration 5 
#> Estimating theta... 
#> Estimation of theta complete 
#> Estimating beta... 
#> Estimation of beta complete 
#> Iteration 5 complete 
#> Current penalized negative log likelihood: 1018.501 
#> Current MSE: 4.220799 

pseq <- create_param_sequence(3, 2)
soil.params <- soil.mvm@covparams
# sigmas
soil.params[pseq[1,1]:pseq[1,2]]
#> [1] 0.7178631 5.8137600 9.5754928
# scale parameters
scale.seq <- pseq[2,1]:pseq[2,2]
# scale parameter for location, outcome 1
soil.params[scale.seq[1]]
#> [1] 13.09338
# scale parameter for elevation, outcome 1
soil.params[scale.seq[2]]
#> [1] 19.73042
# scale parameter for location, outcome 2
soil.params[scale.seq[3]]
#> [1] 18.48061
# scale parameter for elevation, outcome 2
soil.params[scale.seq[4]]
#> [1] 0.4919777
# scale parameter for location, outcome 3
soil.params[scale.seq[5]]
#> [1] 0.5017857
# scale parameter for elevation, outcome 3
soil.params[scale.seq[6]]
#> [1] 0.5202431
# smoothness parameters
soil.params[pseq[3,1]:pseq[3,2]]
#> [1] 0.1079127 0.3079058 0.2777398
# nuggets
soil.params[pseq[4,1]:pseq[4,2]]
#> [1] -0.1754727 -0.3717984 -0.8519659
# correlation
soil.corr <- diag(2) / 2
soil.corr[upper.tri(soil.corr)] <- soil.params[pseq[5,1]:pseq[5,2]]
#> Warning: number of items to replace is not a multiple of replacement length
soil.corr <- soil.corr + t(soil.corr)