Given a pair of two variables (X,Y), how can you generate and plot the empirical joint PDF and CDF in vanilla MATLAB (no toolboxes)?
The data are:
X
, Y
.x_axis
, y_axis
. The edges must obviously be increasing, but need not be uniformly spaced.The resulting PDF and CDF are defined at the centers of the rectangles determined by the x and y edges.
To plot the results in 3D, use surf(...)
instead of imagesc(...)
.
clear all
%// Data (example):
X = randn(1,1e5); %// random variables.
Y = randn(1,1e5);
x_axis = -3:.2:3; %// Define edges of bins for x axis. Column vector
y_axis = -3:.2:3; %// Same for y axis
%// Compute corners of 2D-bins:
[x_mesh_upper,y_mesh_upper] = meshgrid(x_axis(2:end),y_axis(2:end));
[x_mesh_lower,y_mesh_lower] = meshgrid(x_axis(1:end-1),y_axis(1:end-1));
%// Compute centers of 1D-bins:
x_centers = (x_axis(2:end)+x_axis(1:end-1))/2;
y_centers = (y_axis(2:end)+y_axis(1:end-1))/2;
%// Compute pdf:
pdf = mean( bsxfun(@le, X(:), x_mesh_upper(:).') ...
& bsxfun(@gt, X(:), x_mesh_lower(:).') ...
& bsxfun(@le, Y(:), y_mesh_upper(:).') ...
& bsxfun(@gt, Y(:), y_mesh_lower(:).') );
pdf = reshape(pdf,length(x_axis)-1,length(y_axis)-1); %// pdf values at the
%// grid points defined by x_centers, y_centers
pdf = pdf ./ (y_mesh_upper-y_mesh_lower) ./ (x_mesh_upper-x_mesh_lower);
%// normalize pdf to unit integral
%// Compute cdf:
cdf = mean( bsxfun(@le, X(:), x_mesh_upper(:).') ...
& bsxfun(@le, Y(:), y_mesh_upper(:).') );
cdf = reshape(cdf,length(x_axis)-1,length(y_axis)-1);
%// Plot pdf
figure
imagesc(x_centers,y_centers,pdf)
axis xy
axis equal
colorbar
title 'pdf'
%// Plot cdf
figure
imagesc(x_centers,y_centers,cdf)
axis xy
axis equal
colorbar
title 'cdf'
Matlab R2015b includes an histogram2
function that does all the work. It automatically does the normalization to obtain the PDF (given the appropriate input flag), or even the CDF.
Using the same example as above,
clear all
%// Data (example):
X = randn(1,1e5); % random variables.
Y = randn(1,1e5);
x_axis = -3:.2:3; % Define edges of bins for x axis. Column vector
y_axis = -3:.2:3; % Same for y axis
%// Compute and plot pdf
figure
histogram2(X, Y, x_axis, y_axis, 'Normalization', 'pdf')
%// Compute and plot cdf
figure
histogram2(X, Y, x_axis, y_axis, 'Normalization', 'cdf')