Help vectorising for loop for kernel density

조회 수: 1(최근 30일)
I am having a devil of a time figuring out how to vectorise the following for loopp
dens = NaN(length(y),1);
for i=1:length(y)
dens(i) = (1/(n*h))*sum(kernel((Y-y(i))/h));
end
Where Y has 1000 points and y has 50 points, i think perhaps it could be done by subtracting each element in y from Y such that i get a 1000 by 50 matrix. And then i can sum across columns to get the result perhaps ? but that would requre some fairly large matrices for changed dimensions and i also dont know how to do it
any help will be greatly appreciated
  댓글 수: 2
DoVile Last Name:
DoVile Last Name: 2012년 11월 24일
n = 1000 (number of simulated variables = length(Y), h is the bandwidth = 0.15 and the kernel function i have declared already so when the loop calls kernel(x) it should call that function, the function it self is:
K_Gauss =@(z) exp((-z.^2)./2)./sqrt(2*pi);

댓글을 달려면 로그인하십시오.

채택된 답변

Matt Fig
Matt Fig 2012년 11월 24일
편집: Matt Fig 2012년 11월 24일
With these values for Y and y:
Y = rand(10,1)*10;
y = rand(5,1)*10;
This gives the same result as your FOR loop:
D = (1/(n*h))*sum(kernel(bsxfun(@minus,Y.',y)/h),2);
Here is the code I used to check the equality of the two approaches, in case it helps you:
Y = rand(1000,1)*10;
y = rand(50,1)*10;
dens = NaN(length(y),1);
n = length(Y);
h = .15;
kernel = @(z) exp((-z.^2)./2)./sqrt(2*pi);
for i=1:length(y)
dens(i) = (1/(n*h))*sum(kernel((Y-y(i))/h));
end
D = (1/(n*h))*sum(kernel(bsxfun(@minus,Y.',y)/h),2);
isequal(dens,D) % Check for equality.
  댓글 수: 5
DoVile Last Name:
DoVile Last Name: 2012년 11월 25일
It unfortunately is not working, i accepted the answer by mistake :)
I tested the example code you posted and it works fine, however my main code still does not. I have zipped to a rar file here http://www.mediafire.com/?9o846f8f481q7sr.
Thanks alot for your help btw Matt, i am trying to get into the habbit of vectorizing my code.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by