Calculating a 2D joint probability distribution

Question 1

Improving on code in question

Your loops (and the nested dot product) can be eliminated with bsxfun and matrix multiplication as follows:

xcomps = bsxfun(@ge,xy(:,1),int_x);
ycomps = bsxfun(@ge,xy(:,2),int_y);
count_again = double(xcomps).'*double(ycomps); %' 143x143 = 143x1e5 * 1e5x143
count_again_fix = diff(diff(count_again')');

The multiplication step accomplishes the AND and summation done in sum(inds1 .* inds2), but without looping over the density matrix. EDIT: If you use single instead of double, execution time is nearly halved, but be sure to convert your answer to double or whatever is required for the rest of the code. On my computer this takes around 0.5 sec.

Note: With rot90(count_again/size(xy,1),2) you have a CDF, and in rot90(count_again_fix/size(xy,1),2) you have a PDF.

Using accumarray

Another approach is to use accumarray to make the joint histogram after we bin the data.

Starting with int_x, int_y, K1, xy, etc.:

% take (0,1) data onto [1 K1], following A.Dondas approach for easy comparison
ii = floor(xy(:,1)*(K1-eps))+1; ii(ii<1) = 1; ii(ii>K1) = K1;
jj = floor(xy(:,2)*(K1-eps))+1; jj(jj<1) = 1; jj(jj>K1) = K1;

% create the histogram and normalize
H = accumarray([ii jj],ones(1,size(ii,1)));
PDF = H / size(xy,1); % for probabilities summing to 1

On my computer, this takes around 0.01 sec.

The output is the same as A.Donda's converted from sparse to full (full(H)). Although, as he A.Donda pointed out, it is correct to have the dimensions be K1xK1, rather than the size of count_again_fix in the OPs code that was K1+1xK1+1.

To get the CDF, I believe you can simply apply cumsum to each axis of the PDF.

Question 2

I have written a simple mex function which works very well when N is large. Of course it's cheating but still ...

The function is

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    unsigned long int hh, ctrl;       /*  counters                       */
    unsigned long int N, m, n;        /*  size of matrices               */
    unsigned long int *xy;            /*  data                           */
    unsigned long int *count_cells;   /*  joint frequencies              */
    /*  matrices needed */
    mxArray *count_cellsArray;

/*  Now we need to get the data */
    if (nrhs == 3) {
        xy = (unsigned long int*) mxGetData(prhs[0]);
        N = (unsigned long int) mxGetM(prhs[0]);
        m = (unsigned long int) mxGetScalar(prhs[1]);
        n = (unsigned long int) mxGetScalar(prhs[2]);
    }

/*  Then build the matrices for the output */
    count_cellsArray = mxCreateNumericMatrix(m + 1, n + 1, mxUINT32_CLASS, mxREAL);
    count_cells = mxGetData(count_cellsArray);
    plhs[0] = count_cellsArray;

    hh = 0; /* counter for elements of xy */
    /* for all points from 1 to N */
    for(hh=0; hh<N; hh++) {
        ctrl = (m + 1) * xy[N + hh] + xy[hh];
        count_cells[ctrl] = count_cells[ctrl] + 1;
    }
}

It can be saved in a file "joint_dist_points_2D.c", then compiled:

mex joint_dist_points_2D.c

And check it out:

% Data
N = 1e7;    % number of points
xy = rand(N, 2);    % coordinates of points
xy(randi(2*N, 1000, 1)) = 0;    % add some points on one side
xy(randi(2*N, 1000, 1)) = 1;    % add some points on the other side
xy(randi(N, 1000, 1), :) = 0;    % add some points on one corner
xy(randi(N, 1000, 1), :) = 1;    % add some points on one corner
inds= unique(randi(N, 1000, 1)); xy(inds, :) = repmat([0 1], numel(inds), 1);    % add some points on one corner
inds= unique(randi(N, 1000, 1)); xy(inds, :) = repmat([1 0], numel(inds), 1);    % add some points on one corner

% Intervals for rectangles
K1 = ceil(sqrt(N/5));    % number of intervals along x
K2 = ceil(sqrt(N/7));    % number of intervals along y
int_x = [0:(1 / K1):1, 1+eps];    % intervals along x
int_y = [0:(1 / K2):1, 1+eps];    % intervals along y

% Use Statistics Toolbox: hist3
tic
count_cells_hist = hist3(xy, 'Edges', {int_x int_y});
count_cells_hist(end, :) = []; count_cells_hist(:, end) = [];
toc
% Elapsed time is 4.414768 seconds.

% Use mex function
tic
xy2 = uint32(floor(xy ./ repmat([1 / K1, 1 / K2], N, 1)));
count_cells = joint_dist_points_2D(xy2, uint32(K1), uint32(K2));
toc
% Elapsed time is 0.586855 seconds.

% Check: the two solutions are equivalent
all(count_cells_hist(:) == count_cells(:))

Question 3

chappjc's answer and using hist3 are all good, but since I happened to want to have something like this some time ago and for some reason didn't find hist3 I wrote it myself, and I thought I'd post it here as a bonus. It uses sparse to do the actual counting and returns the result as a sparse matrix, so it may be useful for dealing with a multimodal distribution where different modes are far apart – or for someone who doesn't have the Statistics Toolbox.

Application to francesco's data:

K1 = ceil(sqrt(N/5));
[H, xs, ys] = hist2d(xy(:, 1), xy(:, 2), [K1 K1], [0, 1 + eps, 0, 1 + eps]);

Called with output parameters the function just returns the result, without it makes a color plot.

Here's the function:

function [H, xs, ys] = hist2d(x, y, n, ax)

% plot 2d-histogram as an image
%
% hist2d(x, y, n, ax)
% [H, xs, ys] = hist2d(x, y, n, ax)
%
% x:    data for horizontal axis
% y:    data for vertical axis
% n:    how many bins to use for each axis, default is [100 100]
% ax:   axis limits for the plot, default is [min(x), max(x), min(y), max(y)]
% H:    2d-histogram as a sparse matrix, indices 1 & 2 correspond to x & y
% xs:   corresponding vector of x-values
% ys:   corresponding vector of y-values
%
% x and y have to be column vectors of the same size. Data points
% outside of the axis limits are allocated to the first or last bin,
% respectively. If output arguments are given, no plot is generated;
% it can be reproduced by "imagesc(ys, xs, H'); axis xy".


% defaults
if nargin < 3
    n = [100 100];
end
if nargin < 4
    ax = [min(x), max(x), min(y), max(y)];
end

% parameters
nx = n(1);
ny = n(2);
xl = ax(1 : 2);
yl = ax(3 : 4);

% generate histogram
i = floor((x - xl(1)) / diff(xl) * nx) + 1;
i(i < 1) = 1;
i(i > nx) = nx;
j = floor((y - yl(1)) / diff(yl) * ny) + 1;
j(j < 1) = 1;
j(j > ny) = ny;
H = sparse(i, j, ones(size(i)), nx, ny);

% generate axes
xs = (0.5 : nx) / nx * diff(xl) + xl(1);
ys = (0.5 : ny) / ny * diff(yl) + yl(1);

% possibly plot
if nargout == 0
    imagesc(ys, xs, H')
    axis xy
    clear H xs ys
end