


HISTXY 2-Dimensional Density Histogram. COUNTS = HISTXY(X,Y,D), where X and Y are matrices of the same size, returns a D x D matrix COUNTS containing the numbers of points in a binned scatter plot of the entries of X vs the corresponding entries of Y. The elements of X and Y are each independently discretized into D evenly-spaced bin centers, offset such that the extrema are are bin centers. If D is not specified (or is empty) it defaults to 100. COUNTS = HISTXY(X,Y,[D1,D2]) produces a D2 x D1 matrix COUNTS, with the corresponding binning of X and Y. [COUNTS,X_INDS,Y_INDS] = HISTXY(X,Y,D) also returns the bin centers. The density is then visualized with IMAGESC(X_INDS,Y_INDS,COUNTS). [...] = HISTXT(X,Y,D,BANDWIDTH), for scalar BANDWIDTH, convolves COUNTS with an isotropic 2-D Gaussian kernel of standard deviation BANDWIDTH bins. If BANDWIDTH is a two element vector, an anisotropic Gaussian is used with column standard deviation BANDWIDTH(1) bins and row standard deviation BANDWIDTH(2) bins. [...] = HISTXY(...,'log') uses the log of the counts (0's yield -Inf). HISTXY(...) without output arguments produces an image of the density.


0001 function [counts,x_inds,y_inds] = histxy(x, y, varargin) 0002 %HISTXY 2-Dimensional Density Histogram. 0003 % COUNTS = HISTXY(X,Y,D), where X and Y are matrices of the same size, 0004 % returns a D x D matrix COUNTS containing the numbers of points in a 0005 % binned scatter plot of the entries of X vs the corresponding entries 0006 % of Y. The elements of X and Y are each independently discretized 0007 % into D evenly-spaced bin centers, offset such that the extrema are 0008 % are bin centers. If D is not specified (or is empty) it defaults 0009 % to 100. 0010 % 0011 % COUNTS = HISTXY(X,Y,[D1,D2]) produces a D2 x D1 matrix COUNTS, with 0012 % the corresponding binning of X and Y. 0013 % 0014 % [COUNTS,X_INDS,Y_INDS] = HISTXY(X,Y,D) also returns the bin centers. 0015 % The density is then visualized with IMAGESC(X_INDS,Y_INDS,COUNTS). 0016 % 0017 % [...] = HISTXT(X,Y,D,BANDWIDTH), for scalar BANDWIDTH, convolves 0018 % COUNTS with an isotropic 2-D Gaussian kernel of standard deviation 0019 % BANDWIDTH bins. If BANDWIDTH is a two element vector, an anisotropic 0020 % Gaussian is used with column standard deviation BANDWIDTH(1) bins and 0021 % row standard deviation BANDWIDTH(2) bins. 0022 % 0023 % [...] = HISTXY(...,'log') uses the log of the counts (0's yield -Inf). 0024 % 0025 % HISTXY(...) without output arguments produces an image of the density. 0026 0027 0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Parse Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0029 D = [100,100]; logflag = 0; blur = 0; % defaults 0030 0031 x = x(:); M = length(x); 0032 y = y(:); N = length(y); 0033 if (~isnumeric(x) || ~isnumeric(y) || (M~=N)) 0034 error('First 2 input arguments must be numeric matrices of equal size.'); 0035 end 0036 0037 if (length(varargin) > 0) 0038 tail = varargin{end}; 0039 if (ischar(tail) && strcmpi(tail,'log')) % If the last arg was 'log' ... 0040 varargin = varargin(1:(end-1)); % ... chomp it and set a flag. 0041 logflag = 1; 0042 end 0043 if (length(varargin) > 1) % If two args left, ... 0044 tail = varargin{end}; % ... try to chomp bandwidth. 0045 if (isnumeric(tail) && length(tail)<=2) 0046 varargin = varargin(1:(end-1)); 0047 if (length(tail) == 2), sigma = tail; 0048 elseif (length(tail) == 1), sigma = [tail, tail]; 0049 else, error('Bandwidth can not be empty.'); 0050 end 0051 blur = 1; 0052 sigma(sigma == 0) = 0.01; % equiv to 0 b/c 100 stds till next bin => below realmin 0053 end 0054 end 0055 if (length(varargin) > 0) 0056 tail = varargin{end}; 0057 if (isnumeric(tail) && length(tail)<=2) % If remaining arg was a len<=2 vector, ... 0058 varargin = varargin(1:end-1); % ... chomp it and set the bin count. 0059 if (length(tail) == 2), D = tail; 0060 elseif (length(tail) == 1), D = [tail, tail]; 0061 end % if tail was empty, use default ... 0062 end 0063 end 0064 if (length(varargin) > 0), error('Unknown syntax.'); end; 0065 end 0066 0067 if (~isequal(D,round(D))), error('Number of bins D must be integer-valued.'); end; 0068 0069 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Rescale Data %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0070 % Separately scale the x and y data 0071 [x,min1,max1] = rescale(x,1,D(1)); x = round(x); 0072 [y,min2,max2] = rescale(y,1,D(2)); y = round(y); 0073 0074 % Bin centers are equally spaced over the range of the unscaled data 0075 x_inds = linspace(min1,max1,D(1)); 0076 y_inds = linspace(min2,max2,D(2)); 0077 0078 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Histogram %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0079 counts = CORE_histxy(x,y,D(1),D(2)); 0080 0081 if (blur) % use separable kernels for convolution with least odd length >= num bins 0082 x_kernel = exp(-([-floor(D(1)/2):floor(D(1)/2)].^2)./(2*sigma(1).^2)); 0083 x_kernel = x_kernel ./ sum(x_kernel); 0084 y_kernel = exp(-([-floor(D(2)/2):floor(D(2)/2)].^2)./(2*sigma(2).^2))'; 0085 y_kernel = y_kernel ./ sum(y_kernel); 0086 0087 counts = conv2(y_kernel, x_kernel, counts, 'same'); 0088 end 0089 0090 if (logflag) 0091 old = warning('off'); counts = log(counts); warning(old); 0092 end 0093 0094 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Graphics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0095 if (nargout == 0) 0096 imagesc(x_inds, y_inds, counts); axis xy; 0097 clear counts x_inds y_inds % clear these so nothing is dumped to output 0098 end