Nothing
###################################################################################
#' Calculate negative log-likelihood
#'
#' Used to determine theta/lambda/p values for the Kriging model in \code{\link{modelKriging}}
#' with Maximum Likelihood Estimation (MLE).
#'
#' @param xt vector, containing parameters like theta, p and lambda.
#' @param D matrix (or list of multiple matrices) of distances between training samples. In case of multiple distance matrices, theta (part of xt) has to be a vector, giving a weighting parameter for each matrix.
#' @param y vector of observations at sample locations.
#' @param useLambda whether to use nugget effect, i.e., lambda (FALSE at default).
#' @param corr whether to use nugget effect, i.e., lambda (fcorrGauss at default).
#' @param inverter string specifying method for inversion of correlation matrix ("chol", cholesky decomposition at default, any other string leads to using the solve function).
#' @param indefiniteMethod The specific method used for correction: spectrum \code{"clip"}, spectrum \code{"flip"}, spectrum \code{"square"}, spectrum \code{"diffusion"}, feature embedding "feature", nearest definite matrix "near". Default is no correction: \code{"none"}. See Zaefferer and Bartz-Beielstein (2016).
#' @param indefiniteType The general type of correction for indefiniteness: \code{"NSD"},\code{"CNSD"} or the default \code{"PSD"}. See Zaefferer and Bartz-Beielstein (2016).
#' @param indefiniteRepair boolean, whether conditions of the distance matrix (in case of \code{"NSD"},\code{"CNSD"} correction type) or correlation matrix (in case of \code{"PSD"} correction type) are repaired.
#' @param returnLikelihoodOnly boolean, whether the function should return only the likelihood, or a else a list (see return information below).
#' @param inverter string, defining the inverter to use. default \code{"chol"} is inversion via \code{chol2inv}. A different string will lead to use of \code{solve}.
#' @param ntheta number of kernel parameters.
#'
#' @return the numeric Likelihood value (if \code{returnLikelihoodOnly} is TRUE) or a list with elements:
#' \describe{
#' \item{\code{NegLnLike}}{ concentrated log-likelihood *-1 for minimising }
#' \item{\code{Psi}}{ correlation matrix}
#' \item{\code{Psinv}}{ inverse of correlation matrix (to save computation time in forrRegPredictor)}
#' \item{\code{mu}}{ MLE of model parameter mu }
#' \item{\code{yMu}}{ vector of observations y minus mu}
#' \item{\code{SSQ}}{ MLE of model parameter sigma^2}
#' \item{\code{a}}{ transformation vector for eigenspectrum transformation, see Zaefferer and Bartz-Beielstein (2016)}
#' \item{\code{U}}{ Matrix of eigenvectors for eigenspectrum transformation, see Zaefferer and Bartz-Beielstein (2016)}
#' \item{\code{isIndefinite}}{ whether the uncorrected correlation (kernel) matrix is indefinite}
#' }
#'
#' @references Forrester, Alexander I.J.; Sobester, Andras; Keane, Andy J. (2008). Engineering Design via Surrogate Modelling - A Practical Guide. John Wiley & Sons.
#' @references Zaefferer, Martin; Stork, Joerg; Friese, Martina; Fischbach, Andreas; Naujoks, Boris; Bartz-Beielstein, Thomas. (2014). Efficient global optimization for combinatorial problems. In Proceedings of the 2014 conference on Genetic and evolutionary computation (GECCO '14). ACM, New York, NY, USA, 871-878. DOI=10.1145/2576768.2598282
#' @references Zaefferer, Martin; Stork, Joerg; Bartz-Beielstein, Thomas. (2014). Distance Measures for Permutations in Combinatorial Efficient Global Optimization. In Parallel Problem Solving from Nature - PPSN XIII (p. 373-383). Springer International Publishing.
#' @references Martin Zaefferer and Thomas Bartz-Beielstein. (2016). Efficient Global Optimization with Indefinite Kernels. Parallel Problem Solving from Nature-PPSN XIV. Accepted, in press. Springer.
#'
#' @seealso \code{\link{modelKriging}}
#' @keywords internal
###################################################################################
modelKrigingLikelihood <- function(xt,D,y,useLambda=FALSE,corr=fcorrGauss,
indefiniteMethod="none",indefiniteType="PSD",indefiniteRepair=FALSE,
returnLikelihoodOnly=TRUE,inverter = "chol",ntheta=1){
n <- dim(y)[1] #number of observations
defaultPenalty <- n * log(var(y)) + 1e4
isIndefinite <- NA
if(is.list(D)){ #in case of multiple distance matrices
distanceWeights <- 10^xt[ntheta+(1:length(D))]
D <- Reduce("+",mapply("*",D,distanceWeights,SIMPLIFY=FALSE)) #combine matrice by corresponding weight value, and compute sum of the matrices
origD <- D # Fix Definiteness (NSDness, CNSDness) of the provided distance matrix
ret <- correctionDistanceMatrix(D,indefiniteType,indefiniteMethod,indefiniteRepair)
D <- ret$mat
isCNSD <- ret$isCNSD
A <- ret$A
}
if(ntheta<1){ #corr function has no parameters
Psi <- corr(D)
}else{
theta <- xt[1:ntheta]
Psi <- corr(D,theta)
}
if(any(is.infinite(Psi))){ # this is required especially if distance matrices are forced to be CNSD/NSD and hence have zero distances
penalty <- defaultPenalty
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
U <- a <- NA
origPsi <- NA
unrepairedPsi <- NA
if(indefiniteType == "PSD"){
origPsi <- Psi
ret <- correctionKernelMatrix(Psi,indefiniteMethod,indefiniteRepair)
a <- ret$a
U <- ret$U
isIndefinite <- !ret$isPSD
Psi <- ret$mat
unrepairedPsi <- ret$matNoRep
#check whether indef-correction somehow yielded malformed values
if(any(is.na(Psi))){
#warning("NaN or NA values due to failed indefiniteness-correction in (in modelKrigingLikelihood). Returning penalty.")
penalty <- defaultPenalty
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
}
if(useLambda){
lambda <- 10^xt[length(xt)];
Psi <- Psi + diag(lambda,n)
}
if(inverter=="chol"){
## cholesky decomposition
cholPsi <- try(chol(Psi), TRUE)
## give penalty if fail
if(class(cholPsi)[1] == "try-error"){
#warning("Correlation matrix is not positive semi-definite (in modelKrigingLikelihood). Returning penalty.")
penalty <- defaultPenalty - min(eigen(Psi,symmetric=TRUE,only.values=TRUE)$values) #the minimal eigenvalue should push the search towards positive eigenvalues
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
#calculate natural log of the determinant of Psi (numerically more reliable and also faster than using det or determinant)
LnDetPsi <- 2*sum(log(abs(diag(cholPsi))))
#inverse with cholesky decomposed Psi
Psinv <- try(chol2inv(cholPsi),TRUE)
if(class(Psinv)[1] == "try-error"){
#warning("Correlation matrix is not positive semi-definite (in modelKrigingLikelihood). Returning penalty.")
penalty <- defaultPenalty - min(eigen(Psi,symmetric=TRUE,only.values=TRUE)$values) #the minimal eigenvalue should push the search towards positive eigenvalues
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
}else{
Psinv <- try(solve(Psi),TRUE) #inverse with LU decomposition
if(class(Psinv)[1] == "try-error"){
#warning("Correlation matrix is not positive semi-definite (in modelKrigingLikelihood). Returning penalty.")
penalty <- defaultPenalty - min(eigen(Psi,symmetric=TRUE,only.values=TRUE)$values) #the minimal eigenvalue should push the search towards larger eigenvalues
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
#calculate natural log of the determinant of Psi
LnDetPsi <- determinant(Psi,logarithm=TRUE)$modulus
}
## Check whether Psi is ill-conditioned
kap <- 1 / ( max(colSums(abs(Psi))) * max(colSums(abs(Psinv)))) # == rcond(Psi)
if(is.na(kap))
kap <- 0
if(kap < 1e-10){
#warning("Correlation matrix is ill-conditioned (in modelKrigingLikelihood). Returning penalty.")
penalty <- defaultPenalty
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
psisum <- sum(Psinv) #this sum of all matrix elements may sometimes become zero, which may be caused by inaccuracies. then, the following may help
if(psisum==0){
psisum <- as.numeric(rep(1,n) %*% Psinv %*% rep(1,n))
if(psisum==0){ #if it is still zero, return penalty
#warning("Sum of elements in inverse correlation matrix is zero (in modelKrigingLikelihood). Returning penalty.")
penalty <- defaultPenalty
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
}
mu <- sum(Psinv%*%y)/psisum
if(is.infinite(mu)|is.na(mu)){
#warning("MLE estimate of mu is infinite or NaN. Returning penalty.")
penalty <- defaultPenalty
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
yMu <- y-mu
#SigmaSqr <- (t(yMu)%*%Psinv%*%yMu)/n
SigmaSqr <- crossprod(yMu,Psinv)%*%yMu/n
if(SigmaSqr < 0){
#warning("Maximum Likelihood Estimate of model parameter sigma^2 is negative (in modelKrigingLikelihood). Returning penalty. ")
penalty <- as.numeric(defaultPenalty-SigmaSqr)
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
NegLnLike <- n*log(SigmaSqr) + LnDetPsi
if(is.na(NegLnLike)|is.infinite(NegLnLike)){#this may happen eg if all y are 0
penalty <- defaultPenalty
if(returnLikelihoodOnly){
return(penalty)
}
return(list(NegLnLike=penalty,origPsi=NA,Psi=NA,unrepairedPsi=NA,Psinv=NA,mu=NA,yMu=NA,SSQ=NA,a=NA,U=NA,isIndefinite=isIndefinite))
}
if(returnLikelihoodOnly){
return(as.numeric(NegLnLike))
}
ret <- list(NegLnLike=NegLnLike,origPsi=origPsi,Psi=Psi,unrepairedPsi=unrepairedPsi,Psinv=Psinv,mu=mu,yMu=yMu,SSQ=SigmaSqr,a=a,U=U,isIndefinite=isIndefinite)
if(exists("origD")){
ret$D <- D
ret$origD <- origD
ret$isCNSD <- isCNSD
ret$A <- A
}
ret
}
###################################################################################
#' Calculate negative log-likelihood
#'
#' This is a wrapper for the Kriging likelihood function \code{\link{modelKrigingLikelihood}}.
#' It is intended for the case where parameters of the distance function are also optimized
#' during maximum likelihood estimation. Thus, the wrapper receives the data, computes the
#' parameterized distance matrix and passes it to the standard likelihood function.
#'
#' @param xt vector, containing parameters like theta, p and lambda.
#' @param xs training samples, which are the input for the distance function. Should be in list format.
#' @param ys vector of observations at training sample locations.
#' @param useLambda whether to use nugget effect, i.e., lambda (FALSE at default).
#' @param corr whether to use nugget effect, i.e., lambda (fcorrGauss at default).
#' @param inverter string specifying method for inversion of correlation matrix ("chol", cholesky decomposition at default, any other string leads to using the solve function).
#' @param indefiniteMethod The specific method used for correction: spectrum \code{"clip"}, spectrum \code{"flip"}, spectrum \code{"square"}, spectrum \code{"diffusion"}, feature embedding "feature", nearest definite matrix "near". Default is no correction: \code{"none"}. See Zaefferer and Bartz-Beielstein (2016).
#' @param indefiniteType The general type of correction for indefiniteness: \code{"NSD"},\code{"CNSD"} or the default \code{"PSD"}. See Zaefferer and Bartz-Beielstein (2016).
#' @param indefiniteRepair boolean, whether conditions of the distance matrix (in case of \code{"NSD"},\code{"CNSD"} correction type) or correlation matrix (in case of \code{"PSD"} correction type) are repaired.
#' @param returnLikelihoodOnly boolean, whether the function should return only the likelihood, or a else a list (see return information below).
#' @param distanceFunction the distance function.
#' @param combineDistances boolean, whether to combine several distances provided as a list of distance functions.
#' @param distanceParametersLower lower boundary for the distance function(s) parameters. A vector in case of one distance, a list of vectors in case of several functions. The parameters are passed as a vector to each respective distance function.
#' @param ntheta number of kernel parameters.
#' @param scaling boolean, whether to scale the distance matrix.
#'
#' @return the numeric Likelihood value (if \code{returnLikelihoodOnly} is TRUE) or a list with elements:
#' \describe{
#' \item{\code{NegLnLike}}{ concentrated log-likelihood *-1 for minimising}
#' \item{\code{Psi}}{ correlation matrix}
#' \item{\code{Psinv}}{ inverse of correlation matrix (to save computation time in forrRegPredictor)}
#' \item{\code{mu}}{ MLE of model parameter mu }
#' \item{\code{yMu}}{ vector of observations y minus mu}
#' \item{\code{SSQ}}{ MLE of model parameter sigma^2}
#' \item{\code{a}}{ transformation vector for eigenspectrum transformation, see Zaefferer and Bartz-Beielstein (2016)}
#' \item{\code{U}}{ Matrix of eigenvectors for eigenspectrum transformation, see Zaefferer and Bartz-Beielstein (2016)}
#' \item{\code{isIndefinite}}{ whether the uncorrected correlation (kernel) matrix is indefinite}
#' }
#'
#' @seealso \code{\link{modelKrigingLikelihood}}
#' @keywords internal
###################################################################################
modelKrigingParameterizedLikelihood <- function(xt,xs,ys,useLambda=FALSE,corr=fcorrGauss,
indefiniteMethod="none",indefiniteType="PSD",indefiniteRepair=FALSE,returnLikelihoodOnly=TRUE,inverter = "chol",
distanceFunction,combineDistances,distanceParametersLower,ntheta,scaling){
#######
nd <- length(distanceFunction) # number of distance functions
if(combineDistances & nd>1){
nweights=nd #number of weight parameters for combining distances
}else{
nweights=0 #number of weight parameters for combining distances
}
# parameters of the distance function(s)
if(ntheta > 0 | nweights > 0 | useLambda)
distanceParameters <- xt[-(1:(ntheta+nweights+useLambda))]
else
distanceParameters <- xt
#calculate distance matrix
ret <- modelKrigingDistanceCalculation(xs,distanceFunction=distanceFunction,parameters=distanceParameters,
distances=NULL,scaling=scaling,combineDistances=combineDistances,indefiniteMethod=indefiniteMethod,
indefiniteType=indefiniteType,indefiniteRepair=indefiniteRepair,distanceParametersLower)
#Call ordinary likelihood function
ret2 <- modelKrigingLikelihood(xt=xt[1:(ntheta+nweights+useLambda)],ret$D,ys,useLambda,corr,
indefiniteMethod,indefiniteType,indefiniteRepair,returnLikelihoodOnly,inverter,ntheta)
if(returnLikelihoodOnly){
return(ret2)
}
ret2[names(ret)] <- ret
ret2
}
#todo: testing!
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.