Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: ksdensity
Replies: 3   Last Post: Apr 2, 2013 1:57 PM

 Messages: [ Previous | Next ]
 Peter Perkins Posts: 1,665 Registered: 12/7/04
Re: ksdensity
Posted: Jun 5, 2006 9:51 AM

michael wrote:
> Probability Density Estimation of two dimensional data
> i have set of points x,y .
> [1,4]
> [6,1]
> [1,5]
> ....

Below is a very rough version of 2D (Gaussian) kernel density estimation. It

- Peter Perkins
The MathWorks, Inc.

gridx1 = 0:.05:1;
gridx2 = 5:.1:10;
X = [0+.5*rand(20,1) 5+2.5*rand(20,1);
.75+.25*rand(10,1) 8.75+1.25*rand(10,1)];
ksdensity2d(X,gridx1,gridx2);

function f = ksdensity2d(x,gridx1,gridx2,bw)
% KSDENSITY2D Compute kernel density estimate in 2D.
% F = KSDENSITY2D(X,GRIDX,GRIDX2,BW) computes a nonparametric estimate
% of the probability density function of the sample in the N-by-2
% matrix X. F is the vector of density values evaluated at the points
% in the grid defined by the vectors GRIDX1 and GRIDX2. The estimate
% is based on a normal kernel function, using a window parameter
% (bandwidth) that is a function of the number of points in X.
[n,p] = size(x);
m1 = length(gridx1);
m2 = length(gridx2);

% Choose bandwidths optimally for Gaussian kernel
if nargin < 4 || isempty(bw)
sig1 = median(abs(gridx1-median(gridx1))) / 0.6745;
if sig1 <= 0, sig1 = max(gridx1) - min(gridx1); end
if sig1 > 0
bw(1) = sig1 * (1/n)^(1/6);
else
bw(1) = 1;
end
sig2 = median(abs(gridx2-median(gridx2))) / 0.6745;
if sig2 <= 0, sig2 = max(gridx2) - min(gridx2); end
if sig2 > 0
bw(2) = sig2 * (1/n)^(1/6);
else
bw(2) = 1;
end
end

% Compute the kernel density estimate
[gridx2,gridx1] = meshgrid(gridx2,gridx1);
x1 = repmat(gridx1, [1,1,n]);
x2 = repmat(gridx2, [1,1,n]);
mu1(1,1,:) = x(:,1); mu1 = repmat(mu1,[m1,m2,1]);
mu2(1,1,:) = x(:,2); mu2 = repmat(mu2,[m1,m2,1]);
f = sum(normpdf(x1,mu1,bw(1)) .* normpdf(x2,mu2,bw(2)), 3) / n;

% Plot the estimate
% colormap(repmat((256:-1:0)'./256,1,3));
% image(256*f./max(f(:)));
surf(gridx1,gridx2,f);
hold on;
plot3(x(:,1),x(:,2),zeros(n,1),'bo');
hold off;
view(-37.50,30);

Date Subject Author
6/1/06 michael
6/5/06 Peter Perkins
3/28/13 Mike
4/2/13 Mike