Nothing
# Bandwidth selection through generalized cross-validation (GCV) for two-dimension cases.
# The codes in this Rscript refer to the source code of fdapace package with a few modifications.
GCVLwls2DV2 <- function(obsGrid, regGrid, ngrid=NULL, dataType='Sparse', error=rcov$error, kern,
rcov, h0=NULL, verbose=FALSE, CV=FALSE, useBW1SE = FALSE, Lt) {
# Returns: a list of length 2, containing the optimal bandwidth and the gcv score.
# obsGrid: observation points. sort(unique(c(unlist(Lt))))
## Get minimal bandwidth and range
r <- diff(range(obsGrid)) * sqrt(2) # sqrt(2) because the window is circular.
minBW <- GetMinb(Lt, obsGrid = obsGrid, dataType = dataType)
h0 <- minBW
if (is.null(h0)){
stop('the data is too sparse, no suitable bandwidth can be found! Try Gaussian Kernel instead!\n')
}
if (kern == 'gauss') {
h0 = h0 * 0.2;
}
## Get Candidate Bandwidths
h0 <- min(h0, r/4)
if (h0 < r/4) {
q <- (r / (4 * h0)) ^ (1/9)
} else if (h0 < r/2) {
q <- (r / (2 * h0)) ^ (1/9)
} else if (h0 < r) {
q <- (r / h0) ^ (1/9)
} else {
stop('Data is too sparse. The minimal bandwidth is the range of data')
}
bw <- (q ^ seq(0,9,length.out = 10)) * h0 # from h0 to r / 4
## Set GCV/CV options
opth <- h0
leave <- FALSE
iter <- 0
maxIter <- 1
minBWInvalid <- FALSE
while (!leave && iter < maxIter) {
if (minBWInvalid){
minBW <- bw[1]
}
Scores <- matrix(Inf, nrow = length(bw), ncol = 2); colnames(Scores) <- c('SUM','SD');
# try the bandwidths large to small in order to save time due to sparseness in the windows.
for (i in rev(seq_along(bw))) {
h <- bw[i]
Scores[i,'SUM'] <- getGCVscoresV2(h, kern, rcov$tPairs, rcov$cxxn, regGrid=regGrid, verbose=verbose)
if (is.infinite(Scores[i,'SUM'])) {
minBWInvalid <- TRUE
if (i < length(bw)) {
if (minBWInvalid) {
minBW <- bw[i + 1]
minBWInvalid <- FALSE
}
}
break; # This will help break out of the loop if the BW is too small to make sense
}
}
if(is.infinite(min(Scores[,'SUM']))){
opth <- max(bw)
optgcv <- Inf
} else {
ind <- which.min(Scores[,'SUM'])
opth <- bw[ind]
optgcv <- Scores[ind,'SUM']
}
## Check that what we found is coherent.
if (opth >= r - 1e-12) {
minBW <- r
leave <- TRUE
stop('Data is too sparse. The optimal bandwidth equals to the range of input time points. Try Gaussian kernel.')
}
if ( (abs(opth - max(bw)) > 1e-12) && !is.infinite(optgcv))
leave <- TRUE
else if (is.infinite(optgcv)) {
if (verbose)
warning('Data is too sparse, retry with larger bandwidths!')
h0 <- max(bw) * 1.01
} else if ( (abs(opth - max(bw)) > 1e-12) ) {
warning('Optimal bandwidth not found in the candidate bandwidths. Retry with larger bandwidths')
h0 <- max(bw)
}
if (!leave) {
iter <- iter + 1
}
} # The "while (!leave && iter < maxIter) ..." end here
ret <- list(h=opth, gcv=optgcv, minBW=minBW)
if (CV != FALSE)
names(ret)[2] <- 'cv'
return(ret)
}
GetMinb <- function(t, obsGrid, dataType='Sparse', npoly=1, minUniqPts=3, minPts=6) {
if (dataType == 'Sparse') {
dstar <- Minb(obsGrid, 2 + npoly) # rough 1D initial value
n_obs = length(obsGrid);
tmp1 = matrix( rep(0, n_obs^2), ncol = n_obs)
# Find the pair against which we have measurements in the same curve
for (i in 1:length(t)){
idx = match( t[[i]], obsGrid)
tmp1[idx, idx] = 1
}
res = tmp1 - diag(n_obs);
ids = matrix(res > 0);
b = matrix( rep(obsGrid, n_obs), nrow=n_obs);
# Use half of the largest difference between two consequative points in the same
# as curve as your candidate bandwith. We do no worry about the difference
# between to [t_j(end) - t_{1+j}(1)] because this will be negative. This bandwidth tends to be conservative (too large).
# dstar = max(dstar, max(diff(b[ids])/2)); # Original code
dstar = max(dstar, stats::quantile( diff(b[ids]), 0.95)/2 ); # Fix to avoid outliers
}else{
dstar = Minb(obsGrid, 2 + npoly) * 1.5;
}
return(dstar)
}
getGCVscoresV2 <- function(bw, kern, xin, yin, win=NULL, regGrid, RSS=NULL, verbose=FALSE) {
# RSS: for implementing GCV of binned rcov.
# browser()
if (is.null(win))
win <- rep(1, length(yin))
fit <- tryCatch(fdapace::Lwls2D(bw, kern, xin=xin, yin=yin, win=win, xout1=regGrid, xout2=regGrid), error=function(err) {
if (verbose) {
message('Invalid bandwidth. Try enlarging the window size.\n')
}
return(Inf)
})
# Catch
if (is.infinite(fit[1]))
return(Inf)
# workaround for degenerate case.
if (any(is.nan(fit)))
return(Inf)
obsFit <- pracma::interp2(regGrid, regGrid, fit, xin[, 1], xin[, 2])
# residual sum of squares
res <- sum((yin - obsFit) ^ 2 * win)
if (!is.null(RSS))
res <- res + sum(RSS)
# kernel at x=0
k0 <- KernelAt0(kern)
N <- sum(win)
r <- diff(range(xin[, 1]))
bottom <- max(1 - 3 * (1 / N) * (r * k0 / bw)^2, 0)
GCV <- res / bottom^2
return(GCV)
}
KernelAt0 <- function(kern) {
if (kern == 'quar')
k0 <- 0.9375
else if (kern == 'epan')
k0 <- 0.75
else if (kern == 'rect')
k0 <- 0.5
else if (kern == 'gausvar')
k0 <- 0.498677850501791
else if (kern == 'gauss')
k0 <- 0.398942280401433
else
stop('Unknown kernel')
return(k0)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.