Nothing
#····································································
# np.svar.R (npsp package)
#····································································
# np.svar S3 class and methods
# np.svar() S3 generic
# np.svar.default(x, y, h, maxlag, nlags, minlag, degree,
# drv, hat.bin, ncv, ...)
# np.svar.svar.bin(x, h, degree, drv, hat.bin, ncv, ...)
# np.svariso(x, y, h, maxlag, nlags, minlag, degree,
# drv, hat.bin, ncv, ...)
# np.svariso.hcv(x, y, maxlag, nlags, minlag, degree,
# drv, hat.bin, loss, ncv, warn, ...)
# np.svariso.corr(lp, x, h, maxlag, nlags, minlag, degree, drv, hat.bin,
# tol, max.iter, plot, ylim)
#
# (c) R. Fernandez-Casal
#
# NOTE: Press Ctrl + Shift + O to show document outline in RStudio
#····································································
# PENDENTE:
# - svarisohcv o final da documentacion
# - engadir aniso, 2iso, niso
#····································································
#····································································
#' Local polynomial estimation of the semivariogram
#'
#' Estimates a multidimensional semivariogram (and its first derivatives)
#' using local polynomial kernel smoothing of linearly binned semivariances.
#' @aliases np.svar-class
#' @param x object used to select a method. Usually a matrix with the
#' coordinates of the data locations (columns correspond with dimensions and
#' rows with data).
#' @param ... further arguments passed to or from other methods.
#' @details Currently, only isotropic semivariogram estimation is supported.
#'
#' If parameter \code{nlags} is not specified is set to \code{101}.
# If parameter \code{nlags} is not specified is set to \code{max(50, rule.svar(x))}.
#'
#' The computation of the hat matrix of the binned semivariances (\code{hat.bin = TRUE})
#' allows for the computation of approximated estimation variances (e.g. in \code{\link{fitsvar.sb.iso}}).
#'
#' A multiplicative triweight kernel is used to compute the weights.
#' @return Returns an S3 object of class \code{np.svar} (locpol svar + binned svar + grid par.),
#' extends \code{\link{svar.bin}}, with the additional (some optional) 3 components:
#' \item{est}{vector or array with the
#' local polynomial semivariogram estimates. }
#' \item{locpol}{a list of 6 components:
#' \itemize{
#' \item{\code{degree} degree of the local polinomial used.}
#' \item{\code{h} smoothing matrix.}
#' \item{\code{rm} mean of residual semivariances.}
#' \item{\code{rss} sum of squared residual semivariances.}
#' \item{\code{ncv} number of cells ignored in each direction.}
#' \item{\code{hat} (if requested) hat matrix of the binned semivariances.}
#' \item{\code{nrl0} (if appropriate) number of cells with \code{binw > 0}
#' and \code{est == NA}.}
#' }}
#' \item{deriv}{(if requested) matrix of estimated first semivariogram derivatives.}
#' @seealso \code{\link{svar.bin}}, \code{\link{data.grid}}, \code{\link{locpol}}.
#' @references
#' Fernandez Casal R., Gonzalez Manteiga W. and Febrero Bande M. (2003)
#' Space-time dependency modeling using general classes of flexible stationary
#' variogram models, \emph{J. Geophys. Res.}, \bold{108}, 8779,
#' doi:10.1029/2002JD002909.
#'
#' Garcia-Soidan P.H., Gonzalez-Manteiga W. and Febrero-Bande M. (2003)
#' Local linear regression estimation of the variogram,
#' \emph{Stat. Prob. Lett.}, \bold{64}, 169-179.
#' @export
np.svar <- function(x, ...) {
#····································································
UseMethod("np.svar")
} # S3 generic function
# Non parametric pilot estimation of an isotropic semivariogram
# Returns an S3 object of class "np.svar" (extends "svar.bin")
#····································································
#' @rdname np.svar
#' @aliases iso.np.svar
#' @method np.svar default
#' @param y vector of data (response variable).
#' @inheritParams locpol.default
#' @param maxlag maximum lag. Defaults to 55\% of largest lag.
#' @param nlags number of lags. Defaults to 101.
#' @param minlag minimun lag.
#' @param hat.bin logical; if \code{TRUE}, the hat matrix of the binned semivariances is returned.
# @param cov.bin covariance matrix of the binned semivariances.
# Defaults to identity.
#' @export
np.svar.default <- function(x, y, h = NULL, maxlag = NULL, nlags = NULL,
minlag = maxlag/nlags, degree = 1,
drv = FALSE, hat.bin = TRUE, ncv = 0, ...) {
# binning cells without data are set to missing.
# Devuelve estimador np del semivariograma y rejilla binning
# Interfaz para la rutina de fortran "svar_iso_np"
#····································································
y <- as.numeric(y)
ny <- length(y) # number of data
x <- as.matrix(x)
if ( !identical(ny, nrow(x)) )
stop("arguments 'y' and 'x' have incompatible dimensions")
drv <- as.logical(drv)
degree <- as.integer(degree)
if(!(degree %in% 0:2))
stop("argument 'degree' must be 0, 1 or 2")
if (drv && degree==0)
stop("'degree' must be greater than or equal to 1 if 'drv == TRUE'")
# Remove missing values
ok <- complete.cases(x, y) # observations having no missing values across x and y
# ok <- !rowSums(!is.finite(x)) & is.finite(y)
# if (!all(ok)) {
# warning("not finite values removed")
if (any(!ok)) {
warning("missing values removed")
x <- x[ok,]
y <- y[ok]
ny <- length(y)
}
nd <- ncol(x) # number of dimensions
if (is.null(maxlag))
maxlag <- 0.55*sqrt(sum(diff(apply(x, 2, range))^2)) # 55% of largest lag
if (is.null(nlags)) nlags <- 101 # dimension of the binning grid
if(is.null(h)) {
stop("argument 'h' (bandwith) must be provided")
# h <- 1 # bandwith matrix PENDENTE
} else if (!is.numeric(h) || length(h)!= 1L)
stop("bandwith 'h' is not a numeric value")
hat.bin <- as.logical(hat.bin)
hat <- if (hat.bin) double(nlags*nlags) else NA_real_
deriv <- if (drv) rep(NA_real_, nlags) else NA_real_
# Let's go FORTRAN!
# subroutine svar_iso_np( nd, x, ny, y, nlags, minlag, maxlag,
# bin_lag, bin_med, bin_y, bin_w,
# h, lpe, degree, ideriv, deriv, ihat, hatlp,
# ndelcv, rm, rss, nrl0)
ret <-.Fortran("svar_iso_np", nd = as.integer(nd), x = as.double(t(x)),
ny = as.integer(ny), y = as.double(y), nlags = as.integer(nlags),
minlag = as.double(minlag), maxlag = as.double(maxlag),
lag = double(1), med = double(1), biny = double(nlags),
binw = double(nlags), h = as.double(h),
elp = as.double(rep(NA_real_, nlags)), degree = as.integer(degree),
ideriv = as.integer(drv), deriv = deriv, ihat = as.integer(hat.bin),
hat = hat, ncv = as.integer(ncv), rm = double(1), rss = double(1),
nrl0 = integer(1), NAOK = TRUE, PACKAGE = "npsp")
if (ret$nrl0 > 0)
warning("Not enough data in some neighborhoods ('NRL < NINDRL'): ", ret$nrl0)
# Construir o resultado
is.na(ret$biny) <- ret$binw == 0 # biny[binw == 0] <- NA
names(ret$min) <- "h"
result <- with( ret,
data.grid(est = elp, biny = biny, binw = binw,
grid = grid.par(n = nlags, min = minlag, lag = lag, dimnames = "h")) )
result$data <- list(x = x, y = y, med = ret$med)
result$svar <- list(type = "isotropic", estimator = "classical")
result$locpol <- with( ret,
list( degree = degree, h = h, rm = rm, rss = rss, ncv = ncv ))
if (hat.bin) result$locpol$hat <- matrix(ret$hat, nrow = nlags)
if (ret$nrl0 > 0) {
warning("Not enough data in some neighborhoods ('NRL < NINDRL'): ", ret$nrl0)
result$locpol$nrl0 <- ret$nrl0
}
if (drv) result$deriv <- ret$deriv
oldClass(result) <- c("np.svar", "svar.bin", "bin.data", "bin.den", "data.grid")
return(result)
#····································································
} # svarisonp, iso.np.svar, np.svar.default
#····································································
# np.svar.svar.bin(x, h, degree, drv, hat.bin, ncv, ...) ----
#····································································
#' @rdname np.svar
#' @method np.svar svar.bin
#' @export
np.svar.svar.bin <- locpol.svar.bin
#····································································
# np.svariso(x, y, h, maxlag, nlags, minlag, degree, ----
# drv, hat.bin, ncv, ...)
#····································································
#' @rdname np.svar
#' @export
np.svariso <- np.svar.default
#····································································
#' @rdname np.svar
#' @inheritParams h.cv.svar.bin
#' @details \code{np.svariso.hcv} calls \code{\link{h.cv}} to obtain an "optimal"
#' bandwith (additional arguments \code{...} are passed to this function).
#' Argument \code{ncv} is only used here at the bandwith selection stage
#' (estimation is done with all the data).
#' @export
np.svariso.hcv <- function(x, y, maxlag = NULL, nlags = NULL, minlag = maxlag/nlags,
degree = 1, drv = FALSE, hat.bin = TRUE,
loss = c("MRSE", "MRAE", "MSE", "MAE"), ncv = 1, warn = FALSE, ...) {
#····································································
loss <- match.arg(loss)
if (is.null(maxlag))
maxlag <- 0.55*sqrt(sum(diff(apply(x, 2, range))^2)) # 55% of largest lag
if (is.null(nlags)) nlags <- 101 # dimension of the binning grid
bin <- svariso(x, y, maxlag = maxlag, nlags = nlags, minlag = minlag,
estimator = "classical")
hopt <- h.cv.svar.bin(bin, loss = loss, degree = degree,
ncv = ncv , warn = warn, ...)$h
return(locpol(bin, h = hopt, degree = degree, drv = drv, hat.bin = hat.bin))
}
#····································································
#' @rdname np.svar
#' @param lp local polynomial estimate of the trend function (object of class
#' \code{\link{locpol.bin}}).
#' @param tol convergence tolerance. The algorithm stops if the average of the
#' relative squared diferences is less than \code{tol}. Defaults to 0.04.
#' @param max.iter maximum number of iterations. Defaults to 10.
#' @param plot logical; if \code{TRUE}, the estimates obtained at each iteration
#' are plotted.
#' @param verbose logical; if \code{TRUE}, the errors (averages of the
#' relative squared differences) at each iteration are printed.
#' @param ylim y-limits of the plot (if \code{plot == TRUE}).
## @param col colors for lines and points if \code{plot == TRUE}.
#' @details
#' \code{np.svariso.corr} computes a bias-corrected nonparametric semivariogram
#' estimate using an iterative algorithm similar to that described in
#' Fernandez-Casal and Francisco-Fernandez (2014). This procedure tries to correct
#' the bias due to the direct use of residuals (obtained in this case from a
#' nonparametric estimation of the trend function) in semivariogram estimation.
# (additional arguments \code{...} are passed to \code{plot}).
#' @references
#' Fernandez-Casal R. and Francisco-Fernandez M. (2014)
#' Nonparametric bias-corrected variogram estimation under non-constant trend,
#' \emph{Stoch. Environ. Res. Ris. Assess}, \bold{28}, 1247-1259.
#' @export
np.svariso.corr <- function(lp, x = lp$data$x, h = NULL, maxlag = NULL, nlags = NULL,
minlag = maxlag/nlags, degree = 1, drv = FALSE, hat.bin = TRUE,
tol = 0.05, max.iter = 10, plot = FALSE, verbose = plot,
ylim = c(0,2*max(svar$biny, na.rm = TRUE))) {
#····································································
if (!inherits(lp, "locpol.bin"))
stop("function only works for objects of class (or extending) 'locpol.bin'")
if (is.null(lp$locpol$hat))
stop("'lp' must have a '$locpol$hat' component")
ny <- nrow(x)
nd <- ncol(x) # number of dimensions
if (is.null(maxlag))
maxlag <- 0.55*sqrt(sum(diff(apply(x, 2, range))^2)) # 55% of largest lag
if (is.null(nlags)) nlags <- 101 # dimension of the binning grid
if(is.null(h)) {
stop("argument 'h' (bandwith) must be provided")
# h <- 1 # bandwith matrix PENDENTE
} else if (!is.numeric(h) || length(h)!= 1L)
stop("bandwith 'h' is not a numeric value")
lpdat <- predict(lp, hat.data = TRUE)
lp.resid <- lp$data$y - lpdat$y.est
hat.trend <- lpdat$y.hat
svar <- np.svariso(x, lp.resid, h = h, maxlag = maxlag, nlags = nlags,
degree = degree, drv = FALSE, hat.bin = FALSE)
sv.lags <- coords(svar)
if(plot) {
col <- rainbow(max.iter)
plot(sv.lags, svar$biny, xlab = "distance", ylab = "semivariance",
ylim = ylim, col = col[1])
lines(sv.lags, svar$est, col = col[1])
}
svar.biased <- svar$biny
svarold <- 0
dists <- c(0, dist(x)) # 0 + lower triangle of the distance matrix
for (iter in 2:max.iter) {
# iter <- 1; iter <- iter +1
# cov.est <- varcov(svar, coords = x, sill = 2*max(svar$est))
cov.est <- varcov(svar, coords = x)
cov.bias.est <- hat.trend %*% cov.est
cov.bias.est <- cov.bias.est %*% t(hat.trend) - cov.bias.est - t(cov.bias.est)
cov.bias.diag <- diag(cov.bias.est) / 2
svar.bias <- cov.bias.diag + rep(cov.bias.diag, each = ny) - cov.bias.est
svar.bias <- svar.bias[lower.tri(svar.bias)]
tmp <- binning(dists, c(0, svar.bias), nbin = 2*svar$grid$n)
tmp <- approx(coords(tmp), tmp$biny, sv.lags)$y
svar$biny <- svar.biased - tmp
# if (hat.bin && !drv) svar$est <- svar$locpol$hat %*% svar$biny
svar <- locpol(svar, h = h,
degree = degree, drv = drv, hat.bin = hat.bin)
# COIDADO posible division por 0
# PENDENTE ponderar por num de saltos
error <- sqrt(mean((svarold/svar$est - 1)^2, na.rm = TRUE))
if(plot) {
lines(sv.lags, svar$est, col = col[iter])
points(sv.lags, svar$biny, col = col[iter])
}
if (verbose) cat('Iteration ', iter, ': ', error, '\n')
if (error < tol) break
svarold <- svar$est
} # for (iter in 2:sv.niter)
if (iter == max.iter)
warning('The maximun number of iterations has been reached.')
if(plot) {
lines(sv.lags, svar$est, lwd = 2)
points(sv.lags, svar$biny)
legend("topleft", paste("Iteration", seq(iter)), col = c(col[seq(iter-1)], "black"), lty = 1)
}
svar$svar$estimator <- "bias-corrected (residuals based)"
svar$svar$iter <- iter
svar$svar$error <- error
return(svar)
#····································································
} # np.svariso.corr
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.