# How to estimate the mode using non-parametric methods of a 4-variate random vector drawn from a continuous multivariate distribution?

I have a sample of size 10,000 of a 4-variate random vector coming from a (unknown) continuous multivariate distribution.

How can I estimate the mode of this density using nonparametric methods?

I am currently estimating the corresponding density using the sample and a kernel estimator. Then I take the value in the grid that maximises this estimator (there are several theoretical results that justifies this procedure). The packages ks and np are very slow in this case, particularly in the calculation of the bandwidth matrix.

Do you have any suggestion?

• Hmm, I don't think I can help but perhaps if you posted your code there might be something about the way you are using packages that could be optimised? – Peter Ellis Mar 16 '12 at 18:26
• If the code below is the askers', there isn't much to optimize without digging into the code of the packages. – russellpierce Mar 17 '12 at 3:21
• Perhaps this should be on stack overflow? Depends if it is seen as a programming question or a "is there a better stats package/function" question. – Peter Ellis Mar 20 '12 at 4:29

This is not an answer. I am just posting a code that @PeterEllis suggested but the edit was not accepted for some reason.

library(mvtnorm)
library(misc3d)
library(ks)
library(MASS)

# simulated data, toy example, still hard
dat =  mvrnorm(n=10000, rep(0, 6),diag(6))

#Bandwidth matrix,  takes a looooooooong time
H.scv=Hscv(dat)

# Kernel density estimator, how to specify eval.points?
den = kde(dat, H=H.scv, eval.points=???????)

#Point that maximises the density estimator
ind=which(den$estimate==max(den$estimate), arr.ind=TRUE)

• Are you the original asker of this question? – whuber Mar 16 '12 at 19:51
• Yes, it's important: your attempt to edit the question was rejected because you are not the owner of the question! As far as anyone at SE can tell, the two accounts have nothing in common, not even an IP address. When logged on as Corsario, you have only as many privileges to alter material created by Demian as anyone else would with Corsario's reputation. Thus it's in your interest to merge the accounts, but first you would need to prove they have the same owners. – whuber Mar 16 '12 at 22:16

A quick glance suggests that this is over my head... but here I go anyway.

In terms of options the docs make it seem like if you had a prior for Hstart things might go a little faster. Things also look like they might go faster if binned, but I get the sense that binning is forced off for if ncol(dat) > 4. However at four (as in your text, but not the purported author's code) it looks like you could turn on binning. The consequences of this are entirely unknown to me.

Manually debugging the code... I saw that most of the time seems to be lost in Hscv (at least for k < 6 (where my computer will finish this month)) during calls in that function to Hpi, gamse.scv, and ultimately nlm. nlm is being passed around as an optimization function. No nasty loops or easily paralleled apply statements jumped out at me (which didn't/doesn't mean much). Both nlm and optim (the other choice mentioned in the docs) are .Internal, so I don't think we'll speed them up by much.

Trying Rprof: Poking way beyond my level of familiarity, I tried using Rprof in utils to profile the code. But, I got nowhere. Most of the crunching that is happening is related to matrix multiplication. I haven't used it myself yet, but I understand that R 2.14 has 'parallel' built in with parApply, parCapply, and parRapply functions. If you can find a high enough order non-serial loop or apply statement, perhaps replacing the existing calls to be parallel would help and not just load you down with overhead.

Perhaps someone with a beefer computer than mine, experience in R 2.14, and/or more experience profiling can tell you more.

Recently we have suggested a fast consistent mode estimator:

P.S. Ruzankin and A.V. Logachov (2019). A fast mode estimator in multidimensional space. Statistics & Probability Letters

Besides, there is a theme on this site devoted to fast mode estimators: Computationally efficient estimation of multivariate mode