ENH: kernel density estimation with asymmetric kernels, on R+, unit interval

[josef-pkt]

part of #7338

beta, gamma, invgauss and recipinvgauss kernels can be obtained through sccipy's distributions, with appropriate parameterization. Birnbaum-Saunders (fatiguelife) should also be possible but I haven't tried yet.

I don't know if MultivariateKDE is setup for asymmetric kernels. In contrast to symmetric kernels, the kernel does not depend on a the distance t - x, but it's a lonlinear function K(t, x) where x becomes a shape parameter and bandwidth is a scale parameter.

easy:

• pdf: just uses pdf of distribution averaged over sample points
• cdf: AFAICS, this just uses the sf of the distribution averaged over sample points

maybe easy

• add weights so we can use it for binned, histogram data
• rvs: I haven't tried yet, but I guess we can use rvs of distribution in a mixture, (I have something similar for bernstein beta histogram smoothing).

extras

• bandwidth selection
  • IMSE optimal bandwidth, explicit formulas in articles, requires density and first two derivatives of density
  • reference rules of thumb: it's possible to compute optimal bandwidth for a reference distribution on R+ like log-normal or gamma. One reference stated that log-normal reference bw is too small and not recommended in practise. Other articles estimate a gamma distribution and estimate the optimal bandwidth using it as reference distribution
  • cross-validation, .... I didn't look at those, available for MultivariateKDE
• var, mse at a point: also explicit formulas involving density and derivatives

The last part needs kernel specific formulas, i.e. we need kernel classes with extras. Currently, I'm working with just simple functions that compute kernel-pdf or kernel-cdf

other targets

• hazard function
• density derivatives

multivariate

• The univariate kernels can be easily extended to multivariate case using product kernel. as in MultivariateKDE
• There is a literature on bias reduction for multivariate product kernels, and for using non-product multivariate kernel. I guess that's not a priority.

performance, speedups, not clear to me yet

• currently I'm vectorizing in sample points, but loop over evaluation/prediction points
• histogram/binning: I guess for large samples we can switch to binning as in fft kde, except we cannot use convolution as speedup

kernels

• unit interval, cube: beta
• R+: gamma and Birnbaum-Saunders seems to be good, others less so gamma allows unbound prob at zero, invgauss and recipinvgauss force prob(x=0) = 0 designer distributions: using power BS or similar transformation we can get different kernels. I don't know if any of them have a large advantage.

tails

I haven't seen any references yet. But the kernels should imply different tail behaviors, e.g. can we get heavy tails? I guess we might want to choose kernels depending on behavior around x=0 boundary and on the behavior in tails.

status Currently I have mainly the scipy distribution parameterization of the kernels, which gives pdf and cdf (and maybe rvs) They work if I choose the bandwidth by visual inspection.

I would like to park those functions and leave the rest for another year.

(list of references coming later)

[josef-pkt]

found an R package with BS, Gamma, Erlang and LN kernels while looking for Birnbaum-Saunders kernel articles https://cran.r-project.org/web/packages/DELTD/index.html

I haven't done a systematic search for functions in R yet.

[josef-pkt]

browsing current nonparametric kde code

MultivariateKDE does not assume symmetric, distance kernels K(x - xi), aitchison_aitken for categorical does not use simple distance The univariate kde and kernels defined in sandbox uses distance measure |x - xi| and so will not work for asymmetric kernels.

That means that it should be possible to add asymmetric kernels and additional data types to MultivariateKDE (eg. "u" for unit interval and "p" for R+)

[josef-pkt]

binning again

I ran my notebook with all kernels using histogram binning. With 50 bins, several kernels show spikes, gamma and beta look fine in the example. With 100 bins, all kde and kernel-cdf look good, weibull has some wiggles and might need larger bw than in my original example.

examples use nobs=1000, and the binning function, where rvs_ is the original random sample

def get_bins(rvs, bins=100):
    count, edges = np.histogram(rvs, bins=bins)
    center = edges[:-1] + np.diff(edges) / 2
    probs = count / count.sum()
    return center, probs

rvs, weights = get_bins(rvs_)
kde = kern.pdf_kernel_asym(x_plot, rvs, bw, "gamma2", weights=weights)
kce = kern.cdf_kernel_asym(x_plot, rvs, bw, "gamma2", weights=weights)

(this comment was supposed to be in the PR, but ok here)

[josef-pkt]

some kernels require density at zero is zero: f(0) = 0, Those kernels cannot estimate an f(0) > 0. I didn't keep track of which kernels require that (gamma, or gamma2 does not). I'm starting to add references as I see them again

kernels with f(0) = 0

log-normal : Igarashi 2016 with changes to kernel to allow f(0) > 0 (generalized) bs: mentioned in Igarashi 2016

Gaku Igarashi (2016): Weighted log-normal kernel density estimation, Communications in Statistics - Theory and Methods, DOI: 10.1080/03610926.2014.963623

[josef-pkt]

followup : R package evmix beta1 beta2 and gamma1, gamma2 kernel estimators, and several traditional boundary correction methods for symmetric kernels

Hu, Yang, and Carl Scarrott. 2018. “Evmix: An R Package for Extreme Value Mixture Modeling, Threshold Estimation and Boundary Corrected Kernel Density Estimation.” Journal of Statistical Software 84 (1): 1–27. https://doi.org/10.18637/jss.v084.i05.

kdensity also has gamma and beta, and a kernel based on gaussian copula by Jones and Henderson (I didn't read that article) https://cran.r-project.org/web/packages/kdensity/readme/README.html