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.
%
% SEE ALSO: plotmatrix, hist3
% 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