Code covered by the BSD License

# pairwise/joint/bivariate histograms for many variables

### Shalin Mehta (view profile)

Explore correlations between variables with pairwise histograms plotted over grid.

varargout=hist3matrix(observations,binedges)
```function varargout=hist3matrix(observations,binedges)
% HIST3MATRIX Display joint histograms between multiple quantities on a
% square grid.
%   HIST3MATRIX(OBSERVATIONS,BINEDGES) displays joint histograms between
%   variables on a square grid. Joint histogram is a count of how often a
%   pair of variables assume a certain pair of values.
%
%   INPUTS: Column vectors of the matrix OBSERVATIONS are assumed to be
%   variables and row vectors are assumed to be measured values of the
%   variable. Thus, Ns=size(OBSERVATIONS,1) is the number of samples in the
%   measurement and Nv=size(OBSERVATIONS,2) is the number of variables.
%   BINEDGES is a cell array of vectors that specify edges to be used for
%   binning corresponding variables. The function computes joint histogram
%   between all possible pairs of variables using hist3 function.
%
%   If no output arguments are passed, the joint histograms are shown over
%   a matrix grid with self-histograms appearing as bar-plots across the
%   diagonal. EXAMPLE:
%       x=rand(1000,1); xedges=linspace(0,1,11); y=2*rand(1000,1);
%       yedges=linspace(0,2,11); z=x.^2; zedges=linspace(0,1,11);
%       hist3matrix([x y z],{xedges yedges zedges});
%
%   [hists histaxes]=HIST3MATRIX(observations,binedges) returns square cell
%   matrices hists and histaxes, whose each element (including diagonals)
%   are joint histograms and vectors of bin centers. Display of figure is
%   suppressed.
%

% Copyright: Shalin B. Mehta 2012. Revision history: Shalin B. Mehta.
% April 09, 2012, first version using edge-based binning.
% April 16, 2012, Added the linear correlation coefficient and P-value
% estimates as title on each joint histogram when displaying results.

if(~ismatrix(observations))
error('Input should be a 2D matrix whose rows are assumed observations and columns variables.');
end

[Ns Nv]=size(observations);
% Ns: number of samples
% Nv: number of observations

hists=cell(Nv,Nv);
histaxes=cell(Nv,Nv);

for i1=1:Nv
for i2=1:Nv
% Use hist3 to obtain joint histogram.
[jhist jhistaxes]=hist3([observations(:,i1) observations(:,i2)],'edges',{binedges{i1} binedges{i2}});
hists{i1,i2}=(1/Ns)*jhist;
% Normalize so that the area under histogram is 1, i.e., it is an
% experimental joint distribution function.

histaxes{i1,i2}=jhistaxes;

end
end

if(nargout == 0)
[C P]=corrcoef(observations,'rows','complete');
for i1=1:Nv
for i2=1:Nv
subplot(Nv,Nv,Nv*(i1-1)+i2); %Subplot indices are row-wise.
% Plot histograms at each grid position and mention the
% correlation coefficient and P-value.
if(i1==i2)
bar(histaxes{i1,i2}{1},sum(hists{i1,i2},2));
xlim([binedges{i2}(1) binedges{i2}(end)]);
else
imagesc(histaxes{i1,i2}{2},histaxes{i1,i2}{1},hists{i1,i2});
axis xy; axis tight; %colorbar;
title(['C=' num2str(C(i1,i2),2) ',P=' num2str(P(i1,i2),2)]);

xlim([binedges{i2}(1) binedges{i2}(end)]);
ylim([binedges{i1}(1) binedges{i1}(end)]);
end
end
end
elseif(nargout == 1)
varargout{1}=hists;
else
varargout{1}=hists;
varargout{2}=histaxes;
end

end```