Fast multivariate unimodal density estimator I have a sample $\boldsymbol{x}_i$ for $i$ in $1,\dots, n$, from a $d$ dimensional density $f(\boldsymbol{x})$ and I would like to estimate this unknown density. In addition I know that $f(\boldsymbol{x})$ is unimodal, but it can be skewed or fat tailed.
Given that $d$ is around 20, I think that using a standard kernel density estimator (KDE) is not an option (I can use a sample size $n$ of around $10^4$). In addition given that I know that $f(\boldsymbol{x})$ is unimodal, I don't think that all the flexibility of KDEs is needed in my case.
Hence I was looking for a parametric density estimator, that can be fitted reasonably fast. I've read a bit about multivariate skewed normal or student-t distributions but I would like to know if there are other options out there (maybe constrained KDEs or semi-parametric density estimators?).
 A: Another potential approach is the 'nonparanormal' density estimator, which transforms each marginal distribution to standard Normal then estimates a correlation matrix. This is making a Gaussian copula assumption, which is relatively strong, but that's part of the high dimension/finite sample size tradeoff.
A: Recently we have published a paper with a fast mode estimator for multidimensional unimodal distributions. https://doi.org/10.1016/j.spl.2019.108670 That is not the answer to this question, since we estimated the mode only, not the whole density. Besides, our algorithm needs some simple corrections to be applicable for the dimensionality $d=20$. However that shows that there do exist fast consistent estimators for the mode.
For semiparametric density estimators, you may see the paper by Hsu and Wu https://doi.org/10.1016/j.csda.2013.01.018 They use multivariate Box-Cox transform to fit the observed empirical distribution to a normal one. The main problem with the approach is that consistency is proved for (samples from a) lognormal dstribution only. Besides, the authors do not write out the complexity and consider samples of no more than 4096 points.
A: I have just stumbled across this question. Perhaps you'll find my Field Estimator for Arbitrary Spaces (FiEstAS) useful.
The algorithm estimates the continuous probability density field underlying a given discrete data sample in multiple, non-commensurate dimensions, by recursive binary splitting. It is fully described in two publications:

*

*Ascasibar & Binney (2005), MNRAS 356, 872

*Ascasibar (2010), CoPhC, 181, 1438
and I would be grateful if you cited them if you find the software useful. Source code is available in the Github repository linked above: https://github.com/paranoya/FiEstAS
In a nutshell, the idea is to divide the domain, one dimension at a time, until each data point ends up in an individual cell. The densities can be crudely estimated from the cell volumes, or they can be used to seed an adaptive kernel density estimator. I tested the algorithm with a few million particles in 6 dimensions, and it works pretty decently, but I wonder whether it will still provide meaningful results in 20 dimensions...
