Nothing
#····································································
# svarmod.sb.R (npsp package)
#····································································
# kappasb(x, dk)
# disc.sb(nx, dk, rmax)
# fitsvar.sb.iso(esv, dk, nx, rmax, min.contrib, method, iter, tol)
#
# (c) R. Fernandez-Casal
# Created: Apr 2013
#
# NOTE: Press Ctrl + Shift + O to show document outline in RStudio
#····································································
# PENDENTE:
# - documentacion
# - @examples
#····································································
#····································································
# kappasb(x, dk) ----
#····································································
#' Coefficients of an extended Shapiro-Botha variogram model
#'
#' Computes the coefficients of an extended Shapiro-Botha variogram model.
#' @details If \code{dk >= 1}, the coefficients are computed as:
#' \deqn{\kappa_d(x) = (2/x)^{(d-2)/2} \Gamma(d/2) J_{(d-2)/2}(x)}
#' where \eqn{J_p} is the Bessel function of order \eqn{p}.
#' \cr If \code{dk == 0}, the coefficients are computed as:
#' \deqn{\kappa _\infty(x) = e^{-x^2}}
#' (corresponding to a model valid in any spatial dimension).
#' \cr NOTE: some authors denote these functions as \eqn{\Omega_d}.
#'
#' @param x numeric vector (on which the kappa function will be evaluated).
#' @param dk dimension of the kappa function.
#' @return
#' A vector with the coefficients of an extended Shapiro-Botha variogram model.
#' @references
#' Shapiro, A. and Botha, J.D. (1991) Variogram fitting with a general class of
#' conditionally non-negative definite functions. \emph{Computational Statistics
#' and Data Analysis}, \bold{11}, 87-96.
#' @seealso
#' \code{\link{svarmod.sb.iso}}, \code{\link{besselJ}}.
#' @examples
#' kappasb(seq(0, 6*pi, len = 10), 2)
#'
#' curve(kappasb(x/5, 0), xlim = c(0, 6*pi), ylim = c(-1, 1), lty = 2)
#' for (i in 1:10) curve(kappasb(x, i), col = gray((i-1)/10), add = TRUE)
#' abline(h = 0, lty = 3)
#' @export
kappasb <- function(x, dk = 0) {
#····································································
if ( zeros <- any(index <- x < sqrt(.Machine$double.eps)) ) {
dx <- dim(x)
x <- x[!index]
# Alternativamente poderiase facer pmax(x, .Machine$double.eps^0.5)
}
res <- switch( min(dk + 1, 5),
exp(-x*x), # dk = 0
cos(x), # dk = 1
besselJ(x, nu = 0), # dk = 2
sin(x)/x, # dk = 3
gamma(dk/2) * (2/x)^((dk-2)/2) * besselJ(x, nu = (dk-2)/2) # dk >=4
) # switch
if (zeros) {
res[!index] <- res
res[index] <- 1
dim(res) <- dx
}
return(res)
}
#····································································
# disc.sb(nx, dk, rmax) ----
#····································································
#' Discretization nodes of a Shapiro-Botha variogram model
#'
#' Computes the discretization nodes of a `nonparametric' extended Shapiro-Botha
#' variogram model, following Gorsich and Genton (2004), as the scaled roots of
#' Bessel functions.
#' @param nx number of discretization nodes.
#' @param dk dimension of the kappa function (\code{dk >= 1}, see Details below).
#' @param rmax maximum lag considered.
#' @details
#' If \code{dk >= 1}, the nodes are computed as:
#' \deqn{x_i = q_i/rmax; i = 1,\ldots, nx,} where
#' \eqn{q_i} are the first \eqn{n} roots of \eqn{J_{(d-2)/2}}, \eqn{J_p}
#' is the Bessel function of order \eqn{p} and \eqn{rmax}
#' is the maximum lag considered. The computation of the zeros of the Bessel
#' function is done using the efficient algorithm developed by Ball (2000).
#'
#' If \code{dk == 0} (corresponding to a model valid in any spatial dimension),
#' the nodes are computed so the gaussian variogram models involved have
#' practical ranges:
# \deqn{r_i = ( 1 + 1.2(i-1))rmax/nx; i = 1,\ldots, nx.}
#' \deqn{r_i = 2 ( 1 + (i-1) ) rmax/nx; i = 1,\ldots, nx.}
#' @return
#' A vector with the discretization nodes.
#' @references
#' Ball, J.S. (2000) Automatic computation of zeros of Bessel functions and other
#' special functions. \emph{SIAM Journal on Scientific Computing}, \bold{21},
#' 1458-1464.
#'
#' Gorsich, D.J. and Genton, M.G. (2004) On the discretization of nonparametric
#' covariogram estimators. \emph{Statistics and Computing}, \bold{14}, 99-108.
#' @seealso
#' \code{\link{kappasb}}, \code{\link{fitsvar.sb.iso}}.
#' @examples
#' disc.sb( 12, 1, 1.0)
#'
#' nx <- 1
#' dk <- 0
#' x <- disc.sb(nx, dk, 1.0)
#' h <- seq(0, 1, length = 100)
#' plot(h, kappasb(x * h, 0), type="l", ylim = c(0, 1))
#' abline(h = 0.05, lty = 2)
#' @export
disc.sb <- function(nx, dk = 0, rmax = 1) {
#····································································
# NOTA: se llama primero a esta función
if (dk == 0)
# OLLO: equiv. modelos gausianos, poden aparecer inestabilidades
# Nodos de discretizacion "xeometricos"
# return( 1.732/(seq(1, by = -1/nx, length = nx) * rmax) )
# return( sqrt(3)/(seq(1/nx, by = 1.2/nx, length = nx) * rmax) )
# return( sqrt(3)/(seq(1/nx, 1, length = nx) * rmax) )
return( sqrt(3)/(seq(1/nx, 1, length = nx) * 2 * rmax) )
# dk > 1
# Let's go FORTRAN!
# subroutine disc_sbv(nx, x, dim, range)
ret <-.Fortran( "disc_sbv", nx = as.integer(nx), x = double(nx),
dk = as.integer(dk), as.double(rmax))
return(ret$x)
}
#····································································
# fitsvar.sb.iso(esv, dk, nx, rmax, min.contrib, method, iter, tol) ----
#····································································
#' Fit an isotropic Shapiro-Botha variogram model
#'
#' Fits a `nonparametric' isotropic Shapiro-Botha variogram model by WLS through
#' quadratic programming.
#' Following Gorsich and Genton (2004), the nodes are selected as the scaled
#' roots of Bessel functions (see \code{\link{disc.sb}}).
#' @aliases fitsvar
#' @aliases fitsvar-class
#' @param esv pilot semivariogram estimate, a \code{\link{np.svar}}-\code{\link{class}}
#' (or \code{\link{svar.bin}}) object. Typically an output of the function
#' \code{\link{np.svariso}}.
#' @param dk dimension of the kappa function (\code{dk == 0} corresponds to a model
#' valid in any dimension; if \code{dk > 0}, it should be greater than
#' or equal to the dimension of the spatial process \code{ncol(esv$data$x)}).
#' @param nx number of discretization nodes. Defaults to \code{min(nesv - 1, 50)},
#' where \code{nesv} is the number of semivariogram estimates.
#' @param rmax maximum lag considered in the discretization
#' (range of the fitted variogram on output).
#' @param min.contrib minimum number of (equivalent) contributing pairs
#' (pilot estimates with a lower number are ignored, with a warning).
#' @param method string indicating the WLS fitting method to be used
#' (e.g. \code{method = "cressie"}). See "Details" below.
#' @param iter maximum number of interations of the WLS algorithm (used only
#' if \code{method == "cressie"}).
#' @param tol absolute convergence tolerance (used only
#' if \code{method == "cressie"}).
#' @details
#' The fit is done using a (possibly iterated) weighted least squares criterion, minimizing:
#' \deqn{WLS(\theta) = \sum_i w_i[(\hat{\gamma}(h_i)) - \gamma(\theta; h_i)]^2.}
#' The different options for the argument \code{method} define the WLS algorithm used:
#' \describe{
#' \item{\code{"cressie"}}{The default method. The procedure is
#' iterative, with \eqn{w_i = 1} (OLS) used for the first step
#' and with the weights recalculated at each iteration,
#' following Cressie (1985), until convergence: \deqn{w_i =
#' N(h_i)/\gamma(\hat{\theta}; h_i)^2,} where \eqn{N(h_i)}
#' is the (equivalent) number of contributing pairs in the
#' estimation at lag \eqn{h_i}.}
#' \item{\code{"equal"}}{Ordinary least squares: \eqn{w_i = 1}.}
#' \item{\code{"npairs"}}{\eqn{w_i = N(h_i).}}
#' \item{\code{"linear"}}{\eqn{w_i = N(h_i)/h_i^2}
#' (default fitting method in \pkg{gstat} package).}
#' }
#' Function \code{\link[quadprog]{solve.QP}} of \pkg{quadprog} package is used
#' to solve a strictly convex quadratic program. To avoid problems, the Cholesky decomposition
#' of the matrix corresponding to the original problem is computed using \code{\link{chol}} with \code{pivot = TRUE}.
#' If this matrix is only positive semi-definite (non-strictly convex QP),
#' the number of discretization nodes will be less than \code{nx}.
#' @return
#' Returns the fitted variogram model, an object of \code{\link{class}} \code{fitsvar}.
#' A \code{\link{svarmod}} object
# (extending \code{\link{sb.iso}}: a \code{\link{svarmod}}) object
#' with additional components \code{esv} (pilot semivariogram estimate) and \code{fit} containing:
#' \item{u}{vector of lags/distances.}
#' \item{sv}{vector of pilot semivariogram estimates.}
#' \item{fitted.sv}{vector of fitted semivariances.}
#' \item{w}{vector of (least squares) weights.}
#' \item{wls}{value of the objective function.}
#' \item{method}{string indicating the WLS fitting method used.}
#' \item{iter}{number of WLS iterations (if \code{method == "cressie"}).}
#'
#' @references
#' Ball, J.S. (2000) Automatic computation of zeros of Bessel functions and other
#' special functions. \emph{SIAM Journal on Scientific Computing}, \bold{21},
#' 1458-1464.
#'
#' Cressie, N. (1985) Fitting variogram models by weighted least squares.
#' \emph{Mathematical Geology}, \bold{17}, 563-586.
#'
#' Cressie, N. (1993) \emph{Statistics for Spatial Data}. New York. Wiley.
#'
#' Fernandez Casal R., Gonzalez Manteiga W. and Febrero Bande M. (2003)
#' Flexible Spatio-Temporal Stationary Variogram Models,
#' \emph{Statistics and Computing}, \bold{13}, 127-136.
#'
#' Gorsich, D.J. and Genton, M.G. (2004) On the discretization of nonparametric
#' covariogram estimators. \emph{Statistics and Computing}, \bold{14}, 99-108.
#'
#' Shapiro, A. and Botha, J.D. (1991) Variogram fitting with a general class of
#' conditionally non-negative definite functions. \emph{Computational Statistics
#' and Data Analysis}, \bold{11}, 87-96.
#' @seealso
#' \code{\link{svarmod.sb.iso}}, \code{\link{disc.sb}}, \code{\link{plot.fitsvar}}.
#' @examples
#' # Trend estimation
#' lp <- locpol(aquifer[,1:2], aquifer$head, nbin = c(41,41),
#' h = diag(100, 2), hat.bin = TRUE)
#' # 'np.svariso.corr()' requires a 'lp$locpol$hat' component
#'
#' # Variogram estimation
#' esvar <- np.svariso.corr(lp, maxlag = 150, nlags = 60, h = 60, plot = FALSE)
#'
#' # Variogram fitting
#' svm2 <- fitsvar.sb.iso(esvar) # dk = 2
#' svm3 <- fitsvar.sb.iso(esvar, dk = 0) # To avoid negative covariances...
#' svm4 <- fitsvar.sb.iso(esvar, dk = 10) # To improve fit...
#'
#' plot(svm4, main = "Nonparametric bias-corrected semivariogram and fitted models", legend = FALSE)
#' plot(svm3, add = TRUE)
#' plot(svm2, add = TRUE, lty = 3)
#' legend("bottomright", legend = c("NP estimates", "fitted model (dk = 10)", "dk = 0", "dk = 2"),
#' lty = c(NA, 1, 1, 3), pch = c(1, NA, NA, NA), lwd = c(1, 2, 1, 1))
#' @export
#····································································
fitsvar.sb.iso <- function(esv, dk = 4*ncol(esv$data$x), nx = NULL, rmax = esv$grid$max,
min.contrib = 10, method = c("cressie", "equal", "npairs", "linear"),
iter = 10, tol = sqrt(.Machine$double.eps)) {
# PENDENTE:
# - Rounding errors? w <- w/sum(w)
# - rematar documentacion: details, examples, ...
# - Version preliminar, final 'fit.svar.sb' valida para modelos anisotropicos
# - verificar missing values
# - nodes un vector de ptos de discretizacion
#····································································
if (!inherits(esv, "svar.bin"))
stop("function only works for objects of class (or extending) 'svar.bin'.")
# if (esv$svar$type != "isotropic")
if (esv$grid$nd != 1)
stop("pilot variogram estimates 'esv' must be isotropic.")
method <- match.arg(method)
if (method != "cressie") iter <- 1
# if (!requireNamespace(quadprog)) stop("'quadprog' package is required.")
# Let's go...
u <- as.numeric(coords(esv))
if (inherits(esv, "np.svar")) {
# "np.svar" class
v <- esv$est
if (!is.null(esv$locpol$hat)) {
# Aproximacion varianza estilo Cressie (suponiendo independencia) para ajuste wls
# PENDIENTE: ESCRIBIR/REVISAR ESTAS CUENTAS
n <- 1 / with(esv, rowSums(locpol$hat^2 /
pmax(matrix(binw, nrow=grid$n, ncol=grid$n, byrow=TRUE), 1))) # num equivalente de aportacions
} else {
n <- esv$binw # num de aportacions
}
} else {
# "svar.bin" class
v <- esv$biny
n <- esv$binw # num de aportacions
}
if (!all(index <- n >= min.contrib)) {
warning("some pilot semivariogram estimates will be ignored (contrib < min.contrib)")
u <- u[index]
v <- v[index]
n <- n[index]
}
n.esv <- length(u)
# Discretization points
if (is.null(nx)) {
nx <- min(n.esv - 1, 50)
} else {
if (nx >= length(u)) {
warning("'nx' must be less than the number of variogram estimates (corrected)")
nx <- n.esv - 1
}
}
x <- disc.sb( nx, dk, rmax)
# x <- sqrt(3)/u[1:nx]
# x <- sqrt(3)/u[2:(nx + 1)]
# M(1:nesv,1:npar)
# M <- cbind(-outer(u, x, function(u, x) kappasb(u*x, dk)), 1)
M <- cbind( -kappasb(outer(u, x), dk), 1)
# Reescalar pesos
w <- switch(method,
npairs = n,
linear = n/pmax(u^2, sqrt(.Machine$double.eps)),
1 # default: "cressie", "equal"
)
w <- w/sum(w)
# (iterative) quadratic programming
n.par <- nx + 1
Amat <- diag(n.par)
Amat[1:nx, n.par] <- -1
# wls loop
i <- 0
conver <- FALSE
while( (i < iter) & !conver) {
i <- i + 1
# D = t(M) %*% diag(w) %*% M
Dmat <- suppressWarnings(chol(crossprod(sqrt(w) * M), pivot = TRUE))
n.par2 <- attr(Dmat, "rank")
pivot <- attr(Dmat, "pivot")[1:n.par2]
inu <- match(n.par, pivot, nomatch = 0)
if (!inu) stop('the algorithm has failed to achieve an optimal solution.')
# max(abs(crossprod(Dmat[1:n.par2, 1:n.par2]) - crossprod(sqrt(w) * M[, pivot])))
# Invert the upper triangular matrix (is there a more efficient way?)
Dmat <- backsolve(Dmat[1:n.par2, 1:n.par2], diag(n.par2))
# d = t(M) %*% diag(w) %*% v
dvec <- drop((w * v) %*% M[, pivot]) # crossprod?
# solve min(-d^T b + 1/2 b^T D b) with the constraints A^T b >= b_0
res <- quadprog::solve.QP(Dmat, dvec, Amat[pivot, pivot], factorized = TRUE) # bvec <- rep(0, n.par2)
fit <- drop(M[, pivot] %*% res$solution)
if (i > 1) {
# fitted values absolute difference convergence criteria
conver <- max(abs(fit - last.fit)) < tol
w <- n/pmax(fit^2, sqrt(.Machine$double.eps))
w <- w/sum(w)
}
wls <- sum(w*(v-fit)^2)
last.fit <- fit
}
if (method == "cressie") {
iter <- i
if (!conver) warning("the wls algorithm did not converge. \n",
" Maximum number of iterations exceeded; \n",
" the current values may be an approximate solution, \n",
" or 'tol' is too big or 'iter' is too small.")
}
# Return the fitted Shapiro-Botha model
nu <- res$solution[inu]
z <- res$solution[-inu]
x <- x[pivot[-inu]]
# Rounding errors in solve.QP, the constraints A^T b >= b_0 might not (exactly) hold
if (any(pivot <- z < 100*.Machine$double.eps)) {
if (!all(pivot)) {
z <- z[!pivot]
x <- x[!pivot]
} else {
z <- 0
x <- x[1]
}
}
nu <- max(nu, sum(z))
result <- svarmod.sb.iso( dk = dk, x = x, z = z, nu = nu, range = rmax)
# Add fitting details
result$fit <- list(u = u, sv = v, fitted.sv = fit, w = w, wls = wls,
method = method, iter = iter)
result$esv <- esv
oldClass(result) <- c("fitsvar", oldClass(result))
return(result)
#····································································
} # fitsvar.sb.iso
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.