Skip to content

Instantly share code, notes, and snippets.

@zahlenteufel
Last active August 29, 2015 14:00
Show Gist options
  • Save zahlenteufel/0ea490a1b28b7cc62492 to your computer and use it in GitHub Desktop.
Save zahlenteufel/0ea490a1b28b7cc62492 to your computer and use it in GitHub Desktop.
Experimenting with Parzen Window in MATLAB for Density Estimation
% Generate the points
clf
subplot(2,1,1);
axis equal
mu1 = [5 4]';
sigma1 = [2 -0.8; -0.8 1];
mu2 = [0 1]';
sigma2 = [1 -0.06;-0.06 0.5];
chol(sigma1);
chol(sigma2); % verifica definida positiva
P1 = 0.3;
P2 = 1 - P1;
n = 500;
n1 = sum(rand(1, n) <= P1);
n2 = n - n1;
x1 = bsxfun(@plus, sigma1 * randn(2, n1), mu1);
x2 = bsxfun(@plus, sigma2 * randn(2, n2), mu2);
points = [x1 x2];
plot(x1(1,:), x1(2,:), 'r+');
hold on
plot(x2(1,:), x2(2,:), 'b+');
xlim([0 7]);
ylim([0 7]);
M = 50;
xs = linspace(0, 7, M);
ys = linspace(0, 7, M);
pXY = zeros(M);
[X Y] = meshgrid(xs,ys);
%
%h = 0.25;
hold on
for i = 1:M
for j= 1:M
pXY(i,j) = 0;
Hmax = 50;
Hmin = 0.001;
while Hmax - Hmin > 0.1
H = (Hmax+Hmin) / 2;
inside = sum(pointsdist(points, [xs(i) ys(j)]') < H);
if inside > sqrt(n)
Hmax = H;
else
Hmin = H;
end
end
h = H;
%circle(xs(i), ys(j), h);
for k=1:n
dist = norm([xs(i) ys(j)]' - points(:,k));
pXY(i,j) = pXY(i,j) + normpdf(dist/h,0,1);
end
end
end
subplot(2,1,2);
surf(X, Y, -pXY');
title(['h = ', num2str(h)]);
xlim([0 7]);
ylim([0 7]);
axis square
fprintf('n = %d, sqrt(n) = %f \n', n, sqrt(n));
subplot(2,1,1);
hold off
plot(x1(1,:), x1(2,:), 'r+');
hold on
plot(x2(1,:), x2(2,:), 'b+');
xlim([0 7]);
ylim([0 7]);
axis square
box off
while 1
[mx,my] = ginput(1);
for H = 0.01:0.2:10
inside = sum(pointsdist(points, [mx my]') < H);
if inside > sqrt(n)
h = H;
break;
end
end
circle(mx, my, h);
end
function d = pointsdist(points, x)
% calculates the distance from each point to x
difs = bsxfun(@minus, points, x);
d = sqrt(sum(difs.^2,1));
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment