


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 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Parse Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 D = [100,100]; logflag = 0; blur = 0; % defaults 0029 0030 x = x(:); M = length(x); 0031 y = y(:); N = length(y); 0032 if ((M~=N) || ~isnumeric(x) || ~isnumeric(y) || any(isinf(x(:))) || any(isinf(y(:)))) 0033 error('First 2 input arguments must be numeric matrices of equal size with no +/- Inf elements.'); 0034 end 0035 mask = (isnan(x) | isnan(y)); nanflag = any(mask(:)); 0036 if (nanflag), warning('HISTXY:NaNflag', 'NaN elements will be ignored.'); end 0037 0038 if (length(varargin) > 0) 0039 tail = varargin{end}; 0040 if (ischar(tail) && strcmpi(tail,'log')) % If the last arg was 'log' ... 0041 varargin = varargin(1:(end-1)); % ... chomp it and set a flag. 0042 logflag = 1; 0043 end 0044 if (length(varargin) > 1) % If two args left, ... 0045 tail = varargin{end}; % ... try to chomp bandwidth. 0046 if (isnumeric(tail) && length(tail)<=2) 0047 varargin = varargin(1:(end-1)); 0048 if (length(tail) == 2), sigma = tail; 0049 elseif (length(tail) == 1), sigma = [tail, tail]; 0050 else, error('Bandwidth can not be empty.'); 0051 end 0052 blur = 1; 0053 sigma(sigma == 0) = 0.01; % equiv to 0 b/c 100 stds till next bin => below realmin 0054 end 0055 end 0056 if (length(varargin) > 0) 0057 tail = varargin{end}; 0058 if (isnumeric(tail) && length(tail)<=2) % If remaining arg was a len<=2 vector, ... 0059 varargin = varargin(1:end-1); % ... chomp it and set the bin count. 0060 if (length(tail) == 2), D = tail; 0061 elseif (length(tail) == 1), D = [tail, tail]; 0062 end % if tail was empty, use default ... 0063 end 0064 end 0065 if (length(varargin) > 0), error('Unknown syntax.'); end; 0066 end 0067 0068 if (~isequal(D,round(D))), error('Number of bins D must be integer-valued.'); end; 0069 if (any(D==1)), error('Number of bins D must be greater than 1.'); end; 0070 0071 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Rescale Data %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0072 % Separately scale the x and y data 0073 [x,min1,max1] = rescale(x,1,D(1)); x = round(x); 0074 [y,min2,max2] = rescale(y,1,D(2)); y = round(y); 0075 0076 % Bin centers are equally spaced over the range of the unscaled data 0077 x_inds = linspace(min1,max1,D(1)); 0078 y_inds = linspace(min2,max2,D(2)); 0079 0080 % Mask NaNs 0081 if (nanflag), D=D+1; x(mask)=D(1); y(mask)=D(2); end; 0082 0083 0084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Histogram %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0085 counts = CORE_histxy(x,y,D(1),D(2)); 0086 0087 0088 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Clean Up %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0089 if (nanflag), counts(end,:) = []; counts(:,end) = []; D = D-1; end; 0090 0091 if (blur) % use separable kernels for convolution with least odd length >= num bins 0092 x_kernel = exp(-([-floor(D(1)/2):floor(D(1)/2)].^2)./(2*sigma(1).^2)); 0093 x_kernel = x_kernel ./ sum(x_kernel); 0094 y_kernel = exp(-([-floor(D(2)/2):floor(D(2)/2)].^2)./(2*sigma(2).^2))'; 0095 y_kernel = y_kernel ./ sum(y_kernel); 0096 0097 counts = conv2(y_kernel, x_kernel, counts, 'same'); 0098 end 0099 0100 if (logflag), o=warning('off'); counts=log(counts); warning(o); end; 0101 0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Graphics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0103 if (nargout > 0), return; end; 0104 imagesc(x_inds, y_inds, counts); axis xy; 0105 clear counts x_inds y_inds % clear these so nothing is dumped to output