0

I have a code to find the pdf's approximation of a vector based on the formula for kernel estimation: enter image description here
I implemented this formula in the code below (see previous question). However, that code takes long time to run (two loops are used). Could you see the below code and help me to make it faster?

This is the code:

function pdf_est=KDE()
close all;
%%Random values of 20 pixels, range=[1 256]
data=randi([1 256],1,20)-1; %// changed: "-1"

%% Estimate histogram%%%%% 
pdf_est=zeros(1,256);
z=256;

for i=0:z-1 %// changed_ subtracted 1 
    for j=1:length(data)
        pdf_est(i+1)=pdf_est(i+1)+Gaussian(i-data(j)); %// changed: "+1" (twice)
    end
end
%% Plot real histogram 1 to 256; binsize=1;
hold on
plot(0:255, imhist(uint8(data))./length(data),'r'); %// changed: explicit x axis
%% Plot histogram estimation
plot(0:255, pdf_est./length(data),'b'); %// changed: explicit x axis
hold off
function K=Gaussian(x)
   sigma=1; %// change? Set as desired
   K=1./(sqrt(2*pi)*sigma)*exp(-x^2./(2*sigma^2));
Community
  • 1
  • 1
user3051460
  • 1,455
  • 22
  • 58
  • possible duplicate of [Estimate pdf of a vector using Gaussian Kernel](http://stackoverflow.com/questions/28541384/estimate-pdf-of-a-vector-using-gaussian-kernel) – knedlsepp Feb 16 '15 at 17:03
  • You should probably use a filter instead... – knedlsepp Feb 16 '15 at 17:04
  • Don't you have a *real* mathematical formula for this problem? The above has serious problems... – knedlsepp Feb 16 '15 at 17:07
  • @knedlsepp: Yes. But the previous answer is high computational time that I mentioned in my question. The formual is very clear. You can find in the kernel based estimation key world – user3051460 Feb 16 '15 at 17:28
  • You keep stating your formulations are very clear, yet your integral is neither a well defined continuous integral, nor are `x` or `z` defined. – knedlsepp Feb 16 '15 at 17:30
  • 1
    You might want to look into a clean version of this instead: http://de.mathworks.com/help/stats/ksdensity.html – knedlsepp Feb 16 '15 at 17:37
  • Thank knedlsepp. The original formula is continuous integral. However, we know that to easily implement, we can use discrete formula. I just convert it to discrete to make easy implementation. Have some mistake in my formula? – user3051460 Feb 16 '15 at 17:37
  • @knedlsepp: Let see the formula (15) in the paper https://www.dropbox.com/s/0tmr81rcv7z8afz/1-s2.0-S1047320313000473-main.pdf?dl=0 . I think normalize for x. z is number of bin. For image z is 255 – user3051460 Feb 16 '15 at 17:45
  • From what I can tell the paper is talking about 2D-integrals, whereas I can't see this in your 'formula'. So well. Could be correct, but who knows... – knedlsepp Feb 16 '15 at 17:58
  • 3
    @knedlsepp It's not really a _duplicate_ but rather a _follow-up_ to the question you linked. That question asked about some problems with the results the OP was obtaining. After those issues had been solved in the accepted answer, the OP said they wanted to vectorize the code. Since that's a different problem, and the original question had been solved, I suggested that they posted a new question (this one) – Luis Mendo Feb 16 '15 at 18:05
  • @knedlsepp: Let see figure 6. I think they draw a 1D histogram. Actually, the above formual is similar the formula (15). my notation x is I(y) in the (15). They multiple with H function to get pixels inside and outside region and using pixels to find pdf. That is my understand – user3051460 Feb 16 '15 at 18:15

1 Answers1

1

You can get rid of both of those nasty nested loops and then use the hardcoded sigma to have a mega-reduced vectorized solution -

pdf_est = sum(1./(sqrt(2*pi))*exp(-bsxfun(@minus,[0:z-1]',data).^2/2),2)

Or if you would like to have the flexibility to have sigma around, use this -

sum(1./(sqrt(2*pi)*sigma)*exp(-bsxfun(@minus,[0:z-1]',data).^2/(2*sigma^2)),2)

That's all really there is!

Quick tests put this to speedup the original code by 10x!

Divakar
  • 218,885
  • 19
  • 262
  • 358