kernel smoothing given a set of points (produces ECDF, PMF, smoothed CDF and differentiated smoothed PDF) , ... and other useful R stats-y stuff

Rstats.R

smooth_kernels <- function(distances, tol = .Machine$double.eps^0.75){

  d_ecdf <- stats::ecdf(distances)
  d_knots <- stats::knots(d_ecdf)

  d_pmf <- vapply(d_knots, function(x,tol){
    d_ecdf(x+.Machine$double.eps^0.75) - d_ecdf(x-.Machine$double.eps^0.75)
  }, numeric(1), tol = tol)
  
  # might want to check into isotonic regression here to force monotonic increasing fn for smooth CDF
  d_cdf <- lokern::glkerns(d_knots,d_ecdf(d_knots),deriv = 0,korder = 4)
  d_pdf <- lokern::glkerns(d_knots,d_ecdf(d_knots),deriv = 1,korder = 3)

  return(list(
    ecdf=d_ecdf,
    knots=d_knots,
    pmf=d_pmf,
    cdf=d_cdf,
    pdf=d_pdf
  ))
}

# calculate QSD (quasi-stationary distribution) for transition matrix P over transient set (names of absorbing state(s) given by argument abs)
# see eqn 3: Darroch, J. N.; Seneta, E. (1965). "On Quasi-Stationary Distributions in Absorbing Discrete-Time Finite Markov Chains". Journal of Applied Probability. 2 (1): 88–100. doi:10.2307/3211876
qsd_dtmc <- function(M,abs = c("D")){
  P <- M[-which(names(M) %in% abs),-which(names(M) %in% abs)]
  pi <- c(F=1,B=0,R=0,L=0,O=0)
  e <- rep(1,5)
  num <- (pi %*% solve(diag(nrow(P)) - P))
  num / as.numeric((num %*% e))
}