Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » Software » comp.soft-sys.matlab

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

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Peter Perkins

Posts: 1,665
Registered: 12/7/04
Re: ksdensity
Posted: Jun 5, 2006 9:51 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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
may help you get started.

- 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
Read ksdensity
michael
6/5/06
Read Re: ksdensity
Peter Perkins
3/28/13
Read Re: ksdensity
Mike
4/2/13
Read Re: ksdensity
Mike

Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.