#' @title Adaptive bw in one dimension
#'
#' @description Calculate adaptive bandwidths in one dimension
#'
#' @param events A numeric vector representing the moments of occurrence of events
#' @param w The weight of the events
#' @param bw A float, the bandiwdth to use
#' @param kernel_name The name of the kernel to use
#' @keywords internal
adaptive_bw_1d <- function(events, w, bw, kernel_name){
kern_func <- select_kernel(kernel_name)
dens <- sapply(events, function(s){
dists <- abs(s-events)
return(sum(kern_func(dists,bw)*w))
}) / bw
lambdaf <- exp(sum(log(1/sqrt(dens)))/length(events))
bws <- bw * (1/sqrt(dens)) * (1/lambdaf)
return(bws)
}
#' @title Temporal Kernel density estimate
#'
#' @description Calculate the Temporal kernel density estimate based on sampling points in
#' time and events
#' @param events A numeric vector representing the moments of occurrence of events
#' @param w The weight of the events
#' @param samples A numeric vector representing the moments to sample
#' @param bw A float, the bandwidth to use
#' @param kernel_name The name of the kernel to use
#' @param adaptive Boolean
#' @return A numeric vector with the density values at the requested timestamps
#' @export
#' @examples
#' data(bike_accidents)
#' bike_accidents$Date <- as.POSIXct(bike_accidents$Date, format = "%Y/%m/%d")
#' start <- min(bike_accidents$Date)
#' diff <- as.integer(difftime(bike_accidents$Date , start, units = "days"))
#' density <- tkde(diff, rep(1,length(diff)), seq(0,max(diff),1), 2, "quartic")
tkde <- function(events, w, samples, bw, kernel_name, adaptive = FALSE){
## selecting the kernel function
kern_func <- select_kernel(kernel_name)
## calculating the adaptive bw if required
if (adaptive){
bws <- adaptive_bw_1d(events, bw, kern_func)
}else{
bws <- rep(bw, length(events))
}
## calculating the kernel values with the bw
kernel_values <- sapply(samples, function(s){
dists <- abs(s-events)
return(sum(kern_func(dists,bws)*w))
})
return((1/bw)*kernel_values)
}
#' @title Bandwidth selection for Temporal Kernel density estimate by likelihood cross validation
#'
#' @description Calculate the likelihood cross validation score for several bandwidths for the
#' Temporal Kernel density
#' @param events A numeric vector representing the moments of occurrence of events
#' @param w The weight of the events
#' @param bws A numeric vector, the bandwidths to use
#' @param kernel_name The name of the kernel to use
#' @return A vector with the cross validation scores (the higher the better).
#' @export
#' @examples
#' data(bike_accidents)
#' bike_accidents$Date <- as.POSIXct(bike_accidents$Date, format = "%Y/%m/%d")
#' start <- min(bike_accidents$Date)
#' diff <- as.integer(difftime(bike_accidents$Date , start, units = "days"))
#' w <- rep(1,length(diff))
#' scores <- bw_cv_likelihood_calc_tkde(diff, w, seq(10,60,10), "quartic")
bw_cv_likelihood_calc_tkde <- function(events, w, bws, kernel_name){
kern_func <- select_kernel(kernel_name)
n <- length(events)
scores <- sapply(bws, function(bw){
vals <- sapply(1:length(events), function(i){
ei <- events[[i]]
wi <- w[[i]]
dists <- abs(ei-events)
dens <- kern_func(dists,bw)*w
tot_dens <- ((sum(dens) - (kern_func(0,bw)*wi)))
return(tot_dens)
})
score <- ((sum(log(vals)))/n) - log((n-1)*bw)
return(score)
})
return(scores)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.