I happened to need a reliable multivariete kernel density estimation method, and found it in this Matlab code, which turned out to work very well with my data. But I have trouble understanding its workings. Author provides no explanation (he links a paper, but unless I badly misunderstood it, it does not concern said method) nor have I been able to find it elsewhere.

I was hoping someone might recognize it / figure it out and explain it / direct me to resources explaining it.

What I do / don't understand

Once estimation parameters (set of kernel locations μ, weights w and squared bandwidths σ; for the sake of simplicity, I will stick with univariete case) are determined by iteratively evolving initial guess, estimate is provided simply as $$ \widehat{f}(x) = \sum_{i=1}^{\gamma} {\frac{w_i}{\sqrt{2 \pi \sigma_i}} e^{-\frac{(x - \mu_i )^2}{2 \sigma_i }}} $$ where $ \gamma $ is number of kernels used.

The procedure of updating parameters (univariete case, Matlab syntax) can be simplified without a loss of functionality as follows :

function [w,mu,Sig,bandwidth,ent]=regEM(w,mu,Sig,bandwidth,X)
p=zeros(n,gam); psig=p;
for i=1:gam
    xRinv = (X - mu(i)).^2 / Sig(i);
    xSig = (xRinv/s)+eps;
    p(:,i) = exp(-0.5 * xRinv - 0.5 * bandwidth^2 / Sig(i)) * w(i) / sqrt(2 * pi * Sig(i));
    psig(:,i) = p(:,i) .* xSig;
density = sum(p,2); psigd = sum(psig,2);
p = bsxfun(@rdivide, p, density); % normalize classification prob.
ent = sum(log(density)); w = sum(p,1);
for i=find(w>0)
    mu(i) = p(:,i)'*X/w(i); %compute mu's
    Sig(i) = p(:,i)'*((X - mu(i)).^2)/w(i)+bandwidth^2; % compute sigmas
curv = mean(psigd ./ density);
bandwidth = 1/(4*n*(4*pi)^(1/2)*curv)^(1/3);

I see how this resembles finding an optimal bandwidth that minimizes mean integrated square error in basic KDE (as described in this paper). But I don't get how exactly does this work here (deriving expression for MISE is much harder in this case). I don't know how does the 'bandwidth' variable come to play here, how to interpret the additional apperances of it, etc.


Could anyone help me understand motivation behind this method and explain it's workflow?

  • $\begingroup$ If you want to understand KDE, see if you can find accounts by Bernard Silverman, who seems to take special care to discuss the motivation for choices of kernels. The implementation of KDEs in R is the procedure density; I have found the default choice of kernels and bandwidth usually works well, but parameters of density permit one to try alternatives. R documentation for density has a brief, but useful, discussion along with references. $\endgroup$ – BruceET Oct 17 '20 at 21:58
  • $\begingroup$ @BruceET I understand standart KDE (sum over kernels with constant bandwith centered in all data points) fairly well, but I dont quite understand this method which uses fewer kernels, each with different bandwitdh. There's plenty of information about other similar approaches, but I haven't been able to find any about this one (which seems to be best suited for my purposes). $\endgroup$ – Regedin Oct 17 '20 at 22:06
  • $\begingroup$ I will look for references, but promise no helpful results. I have not looked into that particular issue before. $\endgroup$ – BruceET Oct 18 '20 at 0:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.