cpda: Probability Density Approximation

#' Use Traditional KDE or Silverman's KDE-FFT Method to Calculate Likelihood
#' 
#' Our KDE-FFT method uses the following steps to calculate likelihood:
#' \enumerate{
#' \item draws Monte Carlo simulations to construct a simulated histogram,
#' \item uses FFT to transform the histogram to spectral domain,
#' \item applies standard Gaussian kernels to smooth it,
#' \item uses inverse FFT to get simulated PDF, and
#' \item interpolates data linearly on the sPDF to obtain likelihood.
#' }
#' 
#' Although \code{stats::density} also uses KDE-FFT method, our method includes
#' simulation and interpolation steps. 
#'   
#' @param xtilde vector of simulations.
#' @param x vector of data.
#' @param h the bandwidth for KDE. If not given, the function uses Silverman's 
#' \sQuote{rule of thumb}, Silverman (1986, page 48, eqn (3.31)).
#' @param m a multiplier to adjust bandwidth proportationally. Default is 0.8.
#' @param p a precision parameter defining the number of grid as power of 2. 
#' The default grid size is 1024. That p is set at 10 (i.e., 2^10).
#' @param N_s number of simulations. This is to adapt the simulation when one
#' uses PDA to calculate likelihood in a defective probability density function,
#' such as the LBA model or the DDM.  The function will count the number of 
#' simulations from \code{xtilde}.  If \code{xtilde} stores simulations from a 
#' a defective distribution, one has to supply N_s. 
#' @param pw whether to calculate gaussian kernel directly or using FFT 
#' algorithm. \eqn{pw = FALSE} uses KDE-FFT method.
#' @return A vector.
#' @references Holmes, W. (2015). A practical guide to the Probability Density
#' Approximation (PDA) with improved implementation and error characterization.
#' \emph{Journal of Mathematical Psychology}, \bold{68-69}, 13--24,
#' \url{https://dx.doi.org/10.1016/j.jmp.2015.08.006}. \cr
#' 
#' Silverman, B. W. (1982). Algorithm as 176: Kernel density estimation using
#' the fast Fourier transform. Journal of the Royal Statistical Society. Series
#' C (Applied Statistics), 31(1), 93-99. \url{https://dx.doi.org/10.2307/2347084}. 
#' 
#' @seealso
#' \code{\link{density}}, \code{\link[MASS]{bandwidth.nrd}}, 
#' \code{\link{bw.nrd}}, \code{\link{bandwidth.nrd}}.
#' @examples
#' ###################
#' ## Example 1     ##
#' ###################
#' x   <- seq(-3, 3, length.out=100) ## Data
#' sam <- rnorm(1e5)              ## Monte Carlo simulation
#' xlabel <- "Observations"
#' ylabel <- "Density"
#'
#' sam  <- rnorm(n)
#' den1 <- dnorm(x)
#' den2 <- likelihood(x, sam, pw=TRUE)
#' den3 <- likelihood(x, sam, pw=FALSE)
#'
#' par(mar=c(4, 5.3, 0.82, 1))
#' plot(x, den1, type="l", lty="dotted", xlab=xlabel, ylab=ylabel, cex.lab=3,
#'      cex.axis=1.5, lwd=3)
#' 
#' lines(x, den2, col="blue", lty="dashed", lwd=2)
#' lines(x, den3, col="red", lwd=2)
#' 
#' ## Define bandwidth using bw.nrd0
#' bandwidth <- 0.8*bw.nrd0(sam)          
#' den1 <- likelihood(x, sam, h=bandwidth, pw=TRUE)
#' den2 <- likelihood(x, sam, h=bandwidth, pw=FALSE)
#' 
#' plot(x, den1, type="l", lty="dotted", xlab=xlabel, ylab=ylabel, cex.lab=3,
#'      cex.axis=1.5, lwd=3)
#' lines(x, den2, col="blue", lty="dashed", lwd=2)
#' lines(x, den3, col="red", lwd=2)
#' 
#' 
#' ###################
#' ## Example 2     ##
#' ###################
#' rm(list=ls())
#' ## Assuming that this is empirical data
#' y <- rtdists::rLBA(1e3, A=.5, b=1, t0=.25, mean_v=c(2.4, 1.2), sd_v=c(1,1))
#' rt1  <- y[y$response==1,"rt"]
#' rt2  <- y[y$response==2,"rt"]
#' srt1 <- sort(rt1)
#' srt2 <- sort(rt2)
#' 
#' xlabel <- "RT"; ylabel <- "Density"
#' par(mfrow=c(1,2))
#' hist(rt1, "fd", freq=F, xlim=c(0.1,1.5), main="Choice 1", 
#'     xlab=xlabel, ylab=ylabel)
#' hist(rt2, "fd", freq=F, xlim=c(0.1,1.5), main="Choice 2",
#'     xlab=xlabel, ylab=ylabel)
#' par(mfrow=c(1,1))
#' 
#' ## Now draw simulations from another rLBA
#' n <- 1e5
#' samp <- rtdists::rLBA(n, A=.5, b=1, t0=.25, mean_v=c(2.4, 1.2), sd_v=c(1,1))
#' samp1 <- samp[samp[,2]==1, 1]
#' samp2 <- samp[samp[,2]==2, 1]
#' fft1  <- likelihood(srt1, samp1, n=n)  
#' fft2  <- likelihood(srt2, samp2, n=n)
#' 
#' ## Calculate theoretical densities
#' den0 <- rtdists::dLBA(y$rt, y$response, A=.5, b=1, t0=.25, 
#'    mean_v=c(2.4, 1.2), sd_v=c(1, 1))
#'    
#' df0 <- cbind(y, den0)
#' df1 <- df0[df0[,2]==1,]
#' df2 <- df0[df0[,2]==2,]
#' den1 <- df1[order(df1[,1]),3]
#' den2 <- df2[order(df2[,1]),3]
#' plot(srt1,  den1, type="l")
#' lines(srt2, den1)
#' lines(srt1, fft1, col="red")
#' lines(srt2, fft2, col="red")
#' @export
likelihood <- function(xtilde, x, h = 0, m = .8, p = 10, N_s = 0, pw = FALSE) {
  if (pw) {
    out <- lik_pw(xtilde, x, h, m, N_s)[,1]
  } else {
    out <- lik_fft(xtilde, x, h, m, p, N_s)[,1]
  }
  return(out)
}