Nothing
###################################################################################
#' Kriging Prediction
#'
#' Predict with a model fit resulting from \code{\link{modelKriging}}.
#'
#' @param object fit of the Kriging model (settings and parameters), of class \code{modelKriging}.
#' @param x list of samples to be predicted
#' @param ... further arguments, not used
#'
#' @return Returned value depends on the setting of \code{object$predAll}\cr
#' TRUE: list with function value (mean) \code{object$y} and uncertainty estimate \code{object$s} (standard deviation)\cr
#' FALSE:\code{object$y}only
#'
#' @seealso \code{\link{modelKriging}}
#' @seealso \code{\link{simulate.modelKriging}}
#' @export
###################################################################################
predict.modelKriging <- function(object,x,...){
ret <- modelKrigingInternalPredictor(object,x)
psi <- ret$psi
## return value:
res <- list(y=ret$y)
##########################################################################
if (object$predAll){
Psinv <- object$Psinv
lambda <- object$lambda
SigmaSqr <- object$SSQ
if(object$indefiniteType=="PSD" & any(object$indefiniteMethod==c("clip","flip","square","diffusion"))){
if(object$isIndefinite){
if(!object$indefiniteRepair){
Psinv <- object$PsinvA
}
}
}
if(object$reinterpolate & lambda > 0){
SigmaSqr <- object$SSQReint
Psinv <- object$PsinvReint
lambda <- 0
}
# Psinv / PsinvReint has t(A)%*%Psi%*%A included already, if necessary. else, transformation has already been done performed for psi
SSqr <- SigmaSqr*(1+lambda-diag(psi%*%Psinv%*%t(psi)))
s <- sqrt(abs(SSqr))
res$s <- as.numeric(s) #return value
}
res
}
###################################################################################
#' Kriging Simulation
#'
#' (Conditional) Simulate at given locations, with a model fit resulting from \code{\link{modelKriging}}.
#' In contrast to prediction or estimation, the goal is to reproduce the covariance
#' structure, rather than the data itself. Note, that the conditional simulation
#' also reproduces the training data, but
#' has a two times larger error than the Kriging predictor.
#'
#' @param object fit of the Kriging model (settings and parameters), of class \code{modelKriging}.
#' @param nsim number of simulations
#' @param seed random number generator seed. Defaults to NA, in which case no seed is set
#' @param xsim list of samples in input space, to be simulated
#' @param conditionalSimulation logical, if set to TRUE (default), the simulation is conditioned with the training data of the Kriging model.
#' Else, the simulation is non-conditional.
#' @param returnAll if set to TRUE, a list with the simulated values (y) and the corresponding covariance matrix (covar)
#' of the simulated samples is returned.
#' @param ... further arguments, not used
#'
#' @return Returned value depends on the setting of \code{object$simulationReturnAll}
#'
#' @references N. A. Cressie. Statistics for Spatial Data. JOHN WILEY & SONS INC, 1993.
#' @references C. Lantuejoul. Geostatistical Simulation - Models and Algorithms. Springer-Verlag Berlin Heidelberg, 2002.
#'
#' @seealso \code{\link{modelKriging}}, \code{\link{predict.modelKriging}}
#' @export
###################################################################################
simulate.modelKriging <- function(object,nsim=1,seed=NA,xsim,conditionalSimulation=TRUE,returnAll=FALSE,...){
if (!is.na(seed)){
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
#
len <- length(xsim) #number of simulated samples (points
noise <- matrix(rnorm(len*nsim),len, nsim)
#
res <- computeCorrelationMatrix(object,xsim)
covar <- res$psi
#
if(conditionalSimulation){
ret <- modelKrigingInternalPredictor(object,xsim)
y <- ret$y
psi <- ret$psi
covarDifference <- covar - psi %*% object$Psinv %*% t(psi)
eigv <- eigen(object$SSQ *covarDifference,symmetric=T) #eigen decomposition
covarDecomposed <- eigv$vectors %*% diag(sqrt(abs(eigv$values))) %*% eigv$vectors
ysim <- covarDecomposed %*% noise
#and the following adds the simulation part to the predictor
y <- matrix(y,len,nsim) + ysim
}else{
eigv <- eigen(object$SSQ *covar,symmetric=T) #eigen decomposition
covarDecomposed <- eigv$vectors %*% diag(sqrt(abs(eigv$values))) %*% eigv$vectors
y <- object$mu + covarDecomposed %*% noise
}
res$y <- y
if(returnAll)
return(res)
else
return(y)
}
###################################################################################
#' Kriging Prediction (internal)
#'
#' Predict with a model fit resulting from \code{\link{modelKriging}}.
#'
#' @param object fit of the Kriging model (settings and parameters), of class \code{modelKriging}.
#' @param x list of samples to be predicted
#'
#' @return returns a list with:
#' \describe{
#' \item{\code{y}}{predicted values}
#' \item{\code{psi}}{correlations between x and training data}
#' }
#'
#' @seealso \code{\link{simulate.modelKriging}}
#' @seealso \code{\link{predict.modelKriging}}
#' @keywords internal
###################################################################################
modelKrigingInternalPredictor <- function(object,x){
if(!is.list(x))x<-list(x)
xo <- object$x
Psinv <- object$Psinv
n <- length(xo)
#one <- rep(1,n)
mu <- object$mu
yMu <- object$yMu
psi <- matrix(1,length(x),n)
fundist <- object$distanceFunction
if(is.list(fundist)){ # multiple distance functions to be combined
psi <- replicate(length(fundist),psi,simplify=FALSE)
if(object$useDistanceParameters){
indices <- rep(1:length(fundist),sapply(object$distanceParametersLower,length)) #indices assigning each parameter to a distances function
}
for(j in 1:length(fundist)){
for (i in 1:n){
if(!object$useDistanceParameters){
psi[[j]][,i] <- distanceVector(xo[[i]],x,fundist[[j]])
}else{
psi[[j]][,i] <- distanceVector(xo[[i]],x,fundist[[j]],object$distanceParameters[indices==j])
}
}
if(object$scaling){
psi[[j]] <- psi[[j]]/object$maximumDistance[[j]]
}
}
psi <- Reduce("+",mapply("*",psi,object$distanceWeights,SIMPLIFY=FALSE)) #combine result by weighted sum
}else{ #only one distance function
for (i in 1:n){
if(!object$useDistanceParameters){
psi[,i] <- distanceVector(xo[[i]],x,fundist)
}else{
psi[,i] <- distanceVector(xo[[i]],x,fundist,object$distanceParameters)
}
}
if(object$scaling){
psi <- psi/object$maximumDistance
}
}
#
##
#
if(object$indefiniteRepair == 3){ #Required distance data for Weighted Sum Repair.
rd <- 1/psi #reciprocal distances
}else if(object$indefiniteRepair == 4){ #required indices of nearest neighbours for NN repair
rd1 <- apply(psi,1,function(x) which(x==min(x))) #which are the nearest neighbors to the predicted points
}
#
##
#
if(any(object$indefiniteMethod==c("clip","flip","near","square","diffusion"))){ #Distance matrix was not CNSD, and a suitable correction method was chosen.
if(object$indefiniteType=="NSD" & !object$indefiniteRepair){ #no repair, NSD-correction: transformation can be used directly
if(!object$isCNSD)
psi <- psi %*% t(object$A)
}
if(object$indefiniteType=="CNSD" | (object$indefiniteType=="NSD" & object$indefiniteRepair)){
if(!object$isCNSD){
if(object$indefiniteRepair == 1 | object$indefiniteType=="CNSD"){ #1 is augmented repair
if(object$indefiniteMethod!="near") #this could be removed, but is currently computationally too expensive
psi <- correctionAugmentedDistanceVector(psi,object,x) #in case of repair, CNSD-correction: retransform with augmented distance matrix
}else if(object$indefiniteRepair == 2){#nystroem. will not work for CNSD, lacks linear transform.
#nystroem approximation instead of correctionAugmented*
#so: first compute transformed (not repaired) psi, then inverse of the actual, transformed (but not repaired) Psi matrix,
# (this may require saving that matrix separately, because A may already be in it.)
#then use product to get self-similarities. use as divider for transformed psi.
psi <- psi %*% t(object$A)
add <- diag(psi %*% ginv(object$matNoRep) %*% t(psi)) #nystroem. only diagonal needed as divider
add <- matrix(add,nrow(psi),nrow(object$matNoRep))
add2 <- matrix(diag(object$matNoRep),nrow(psi),nrow(object$matNoRep),byrow=T)
psi <- 2*psi - add - add2 #repair
#correct lower bound add by
rdif <- rowSums(psi<0)>0 # rows with at least one negative value
if(any(rdif))
psi[rdif,] <- psi[rdif,,drop=FALSE] - apply(psi[rdif,,drop=FALSE],1,min)
}else if(object$indefiniteRepair == 3){ #weighted sum repair. will not work for CNSD, lacks linear transform.
rd1 <- rd / rowSums(rd) #scale by row sum, such that rows sum to 1
if(any(is.na(rd1))){ #inf values produce NAs that should be 1
nas <- is.na(rd1)
naIndex <- which(nas)
rd1[naIndex] <- (nas / rowSums(nas))[naIndex]
}
psi <- psi %*% t(object$A)
#rd is a weight vecor. need weighted sum of diagonal of unrepaired distance matrix (after correction, before repair.)
add <- rd1 %*% diag(object$matNoRep)
add <- matrix(add,nrow(psi),nrow(object$matNoRep))
add2 <- matrix(diag(object$matNoRep),nrow(psi),nrow(object$matNoRep),byrow=T)
#lower bound for add (the estimated repair factor): 2*psi - add2 (add2 is the known repair factor of the observations)
psi <- 2*psi - add2 - add #repair procedure
#correct lower bound add by
rdif <- rowSums(psi<0)>0 # rows with at least one negative value
if(any(rdif))
psi[rdif,] <- psi[rdif,,drop=FALSE] - apply(psi[rdif,,drop=FALSE],1,min)
}else if(object$indefiniteRepair == 4){ #nearest neighbour repair. will not work for CNSD, lacks linear transform.
psi <- psi %*% t(object$A) #transform psi
#need mean of diagonal values indexed by rd1, of unrepaired distance matrix (after correction, before repair.)
diagNN <- diag(object$matNoRep)
meanidx <- function(x){
return(mean(diagNN[x]))
}
meanDiagNN <- sapply(rd1,meanidx) #get the mean diagonal elements of nearest neighbors
add <- matrix(meanDiagNN,nrow(psi),nrow(object$matNoRep)) #the NN elements are the first added factor
add2 <- matrix(diag(object$matNoRep),nrow(psi),nrow(object$matNoRep),byrow=T) #the elements of the matNoRep are the second element.
psi <- 2*psi - add2 - add #repair procedure
##
#lower bound for add (the estimated repair factor): 2*psi - add2 (add2 is the known repair factor of the observations)
#correct lower bound add by
rdif <- rowSums(psi<0)>0 # rows with at least one negative value
if(any(rdif))
psi[rdif,] <- psi[rdif,,drop=FALSE] - apply(psi[rdif,,drop=FALSE],1,min)
}
}
}
}
if((object$indefiniteType=="CNSD" | object$indefiniteType=="NSD") & object$indefiniteMethod=="feature"){ #distances as features
if(!object$isCNSD){
tempx <- split(psi,seq(nrow(psi)))
for (i in 1:n)
psi[,i] <- distanceVector(object$origD[i,],tempx,distanceRealEuclidean) #todo choice of distance
}
}
#
##
#
if(is.null(object$theta)) #corr function has no parameters
psi <- object$corr(psi)
else
psi <- object$corr(psi,object$theta)
if(object$indefiniteType=="PSD" & any(object$indefiniteMethod==c("clip","flip","near","square","diffusion"))){
if(object$isIndefinite){
if(!object$indefiniteRepair){
#psi <- psi %*% t(object$A) #This is already included in Psinv. do nothing.
}else{
if(object$indefiniteRepair == 1){
psi <- correctionAugmentedKernelVector(psi,object,x)
}else{
if(object$indefiniteRepair == 2){ # Nystroem repair.
#todo: note: this case is not working. zero variances.
#psi <- psi %*% t(object$A) #not needed, contained in unrepairedAPsinvA / unrepairedADivier
div <- diag(psi %*% object$unrepairedAPsinvA %*% t(psi)) #nystroem. only diagonal needed as divider
div <- diag(1/sqrt(div),nrow(psi))
psi <- div %*% psi %*% object$ADividedSqrtDiagPsi #repair
}else if(object$indefiniteRepair == 3){ #Weighted Sum Repair.
rd1 <- rd / rowSums(rd) #scale by row sum, such that rows sum to 1
if(any(is.na(rd1))){ #inf values produce NAs that should be 1
nas <- is.na(rd1)
naIndex <- which(nas)
rd1[naIndex] <- (nas / rowSums(nas))[naIndex]
}
#rd is a weight vecor. need weighted sum of diagonal of unrepaired Psi (after correction, before repair.)
rd2 <- rd1 %*% object$diagUnrepairedPsi
div <- diag(1/sqrt(as.numeric(rd2)),nrow(psi))
psi <- div %*% psi %*% object$ADividedSqrtDiagPsi #repair
}else if(object$indefiniteRepair == 4){ #nearest neighbour repair. will not work for CNSD, lacks linear transform.
#need mean of diagonal values indexed by rd1, of unrepaired distance matrix (after correction, before repair.)
diagNN <- object$diagUnrepairedPsi
meanidx <- function(x){
return(mean(diagNN[x]))
}
meanDiagNN <- sapply(rd1,meanidx) #get the mean diagonal elements of nearest neighbors
div <- diag(1/sqrt(as.numeric(meanDiagNN)),nrow(psi))
psi <- div %*% psi %*% object$ADividedSqrtDiagPsi #repair
}
psi[psi > 1] <- 1 #make sure that psi stays in bounds. to avoid numerical issues. see "NOTE" above for distance case.
psi[psi < -1] <- -1 #no guarantees that approximated repair keeps bounds intact for new data?
}
}
}
}
y <- as.numeric(psi%*%Psinv%*%yMu)+mu #todo: Psinv%*%yMu can be precomputed for speedup
list(y=y,psi=psi)
}
###################################################################################
#' Compute Correlation Matrix
#'
#' Compute the correlation matrix of samples x, given the model object.
#'
#' @param object fit of the Kriging model (settings and parameters), of class \code{modelKriging}.
#' @param x list of samples / data
#'
#' @return the correlation matrix
#'
#' @seealso \code{\link{simulate.modelKriging}}
#' @seealso \code{\link{predict.modelKriging}}
#' @keywords internal
###################################################################################
computeCorrelationMatrix <- function(object,x){
if(!is.list(x))x<-list(x)
if(is.null(object$distanceParameters))
object$distanceParameters <- NA
ret <- modelKrigingDistanceCalculation(x,object$distanceFunction,parameters=object$distanceParameters,
NULL,object$scaling,object$combineDistances,object$indefiniteMethod,object$indefiniteType,object$indefiniteRepair,object$distanceParametersLower)
psi <- ret$D
#
##
#
if(is.null(object$theta)) #corr function has no parameters
psi <- object$corr(psi)
else
psi <- object$corr(psi,object$theta)
#
##
#
ret$U <- NA
ret$a <- NA
ret$isIndefinite <- NA
ret$origPsi <- NA
if(object$indefiniteType=="PSD" & any(object$indefiniteMethod==c("clip","flip","near","square","diffusion"))){
#psi <- correctionKernelMatrix(psi,object$indefiniteMethod,object$indefiniteRepair)$mat
ret$origPsi <- psi
ret2 <- correctionKernelMatrix(psi,object$indefiniteMethod,object$indefiniteRepair)
ret$a <- ret2$a
ret$U <- ret2$U
ret$A <- ret2$A
ret$isIndefinite <- !ret2$isPSD
psi <- ret2$mat
}
#
##
#
if(object$useLambda){
psi <- psi + diag(object$lambda,length(x))
}
#
##
#
ret$psi <- psi
ret
}
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.