#' @rdname CV
#' @export
LIK.density <- function(pp,hlim=NULL,hseq=NULL,resolution=64,edge=TRUE,auto.optim=TRUE,
type=c("fixed","adaptive"),seqres=30,parallelise=NULL,
zero.action=0,verbose=TRUE,...){
if(!inherits(pp,"ppp")) stop("data object 'pp' must be of class \"ppp\"")
W <- Window(pp)
if(is.null(hlim)){
ppu <- pp
marks(ppu) <- NULL
md <- min(nndist(unique(ppu)))
hlim <- c(md,max(md*50,min(diff(W$xrange),diff(W$yrange))/6))
} else {
hlim <- checkran(hlim,"'hlim'")
}
if(!zero.action%in%((-1):2)) stop("invalid 'zero.action'")
typ <- type[1]
if(typ=="fixed"){
if(auto.optim){
if(verbose) cat("Searching for optimal h in ",prange(hlim),"...",sep="")
result <- suppressWarnings(optimise(LIK.density.spatial.single,interval=hlim,pp=pp,res=resolution,edge=edge,za=zero.action,maximum=TRUE)$maximum)
if(verbose) cat("Done.\n")
} else {
if(is.null(hseq)) hseq <- seq(hlim[1],hlim[2],length=seqres)
hn <- length(hseq)
if(is.null(parallelise)){
lik.vec <- rep(NA,hn)
if(verbose) pb <- txtProgressBar(1,hn)
for(i in 1:hn){
lik.vec[i] <- LIK.density.spatial.single(hseq[i],pp,resolution,edge,za=zero.action)
if(verbose) setTxtProgressBar(pb,i)
}
if(verbose) close(pb)
} else {
ncores <- detectCores()
if(verbose) cat(paste("Evaluating criterion on",parallelise,"/",ncores,"cores..."))
if(parallelise>ncores) stop("cores requested exceeds available count")
registerDoParallel(cores=parallelise)
lik.vec <- foreach(i=1:hn,.packages="spatstat",.combine=c) %dopar% {
return(LIK.density.spatial.single(hseq[i],pp,resolution,edge,zero.action))
}
if(verbose) cat("Done.\n")
}
result <- cbind(hseq,lik.vec)
dimnames(result)[[2]] <- c("h","CV")
}
} else if(typ=="adaptive"){
ellip <- list(...)
if(is.null(ellip$hp)){
if(verbose) cat("Selecting pilot bandwidth...")
hp <- LSCV.density(pp,verbose=FALSE,zero.action=zero.action)
if(verbose) cat(paste("Done.\n [ Found hp =",hp,"]\n"))
} else {
hp <- ellip$hp
}
if(is.null(ellip$pilot.density)){
pilot.density <- pp
} else {
pilot.density <- ellip$pilot.density
}
if(is.null(ellip$gamma.scale)){
gamma.scale <- "geometric"
} else {
gamma.scale <- ellip$gamma.scale
}
if(is.null(ellip$trim)){
trim <- 5
} else {
trim <- ellip$trim
}
if(is.null(ellip$dimz)){
dimz <- resolution
} else {
dimz <- ellip$dimz
}
if(verbose) cat("Computing multi-scale estimate...")
hhash <- mean(hlim)
msobject <- multiscale.density(pp,h0=hhash,hp=hp,h0fac=hlim/hhash,edge=ifelse(edge,"uniform","none"),resolution=resolution,dimz=dimz,gamma.scale=gamma.scale,trim=trim,intensity=FALSE,pilot.density=pilot.density,verbose=FALSE)
if(verbose) cat("Done.\n")
h0range <- range(as.numeric(names(msobject$z)))
if(auto.optim){
if(verbose) cat("Searching for optimal h0 in ",prange(h0range),"...",sep="")
h0opt <- suppressWarnings(optimise(ms.loo.lik,interval=h0range,object=msobject,za=zero.action,maximum=TRUE)$maximum)
if(verbose) cat("Done.\n")
return(h0opt)
} else {
if(is.null(hseq)) hseq <- seq(h0range[1],h0range[2],length=seqres)
hn <- length(hseq)
lik.vec <- rep(NA,hn)
if(verbose) pb <- txtProgressBar(1,hn)
for(i in 1:hn){
lik.vec[i] <- ms.loo.lik(hseq[i],msobject,zero.action)
if(verbose) setTxtProgressBar(pb,i)
}
if(verbose) close(pb)
result <- cbind(hseq,lik.vec)
dimnames(result)[[2]] <- c("h0","CV")
}
} else stop("invalid 'type'")
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.