3

正方形の中に多くの点があります。正方形を多くの小さな長方形に分割し、各長方形に含まれるポイントの数を確認したい、つまり、ポイントの同時確率分布を計算したい。私は、ループを使用し、あまり効率的ではない、いくつかの常識的なアプローチを報告しています。

% Data
N = 1e5;    % number of points
xy = rand(N, 2);    % coordinates of points
xy(randi(2*N, 100, 1)) = 0;    % add some points on one side
xy(randi(2*N, 100, 1)) = 1;    % add some points on the other side
xy(randi(N, 100, 1), :) = 0;    % add some points on one corner
xy(randi(N, 100, 1), :) = 1;    % add some points on one corner
inds= unique(randi(N, 100, 1)); xy(inds, :) = repmat([0 1], numel(inds), 1);    % add some points on one corner
inds= unique(randi(N, 100, 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 = K1;    % 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

% First approach
tic
count_cells = zeros(K1 + 1, K2 + 1);
for k1 = 1:K1+1
  inds1 = (xy(:, 1) >= int_x(k1)) & (xy(:, 1) < int_x(k1 + 1));
  for k2 = 1:K2+1
    inds2 = (xy(:, 2) >= int_y(k2)) & (xy(:, 2) < int_y(k2 + 1));
    count_cells(k1, k2) = sum(inds1 .* inds2);
  end
end
toc
% Elapsed time is 46.090677 seconds.

% Second approach
tic
count_again = zeros(K1 + 2, K2 + 2);
for k1 = 1:K1+1
  inds1 = (xy(:, 1) >= int_x(k1));
  for k2 = 1:K2+1
    inds2 = (xy(:, 2) >= int_y(k2));
    count_again(k1, k2) = sum(inds1 .* inds2);
  end
end
count_again_fix = diff(diff(count_again')');
toc
% Elapsed time is 22.903767 seconds.

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

時間、メモリ、およびおそらくループを回避するという点で、より効率的に行うにはどうすればよいですか?

編集 -->私もこれを見つけました。これは、これまでに見つかった最良の解決策です。

tic
count_cells_hist = hist3(xy, 'Edges', {int_x int_y});
count_cells_hist(end, :) = []; count_cells_hist(:, end) = [];
toc
all(count_cells(:) == count_cells_hist(:))
% Elapsed time is 0.245298 seconds.

ただし、Statistics Toolbox が必要です。

編集 --> chappjc によって提案されたテスト ソリューション

tic
xcomps = single(bsxfun(@ge,xy(:,1),int_x));
ycomps = single(bsxfun(@ge,xy(:,2),int_y));
count_again = xcomps.' * ycomps; %' 143x143 = 143x1e5 * 1e5x143
count_again_fix = diff(diff(count_again')');
toc
% Elapsed time is 0.737546 seconds.
all(count_cells(:) == count_again_fix(:))
4

3 に答える 3

2

N が大きい場合に非常にうまく機能する単純な mex 関数を作成しました。もちろん詐欺ですが…

機能は

#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;
    }
}

ファイル「joint_dist_points_2D.c」に保存してからコンパイルできます。

mex joint_dist_points_2D.c

そしてそれをチェックしてください：

% 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(:))
于 2013-11-03T12:19:08.373 に答える
0

chappjcの答えと使用法hist3はすべて良いですが、たまたまこのようなものが欲しかったので、何らかの理由hist3で自分で書いたので、おまけとしてここに投稿しようと思いました。実際のカウントを行うために使用sparseし、疎行列として結果を返すため、さまざまなモードが遠く離れているマルチモーダル分布を処理する場合や、Statistics Toolbox を持っていない人にとって役立つ場合があります。

francesco のデータへの適用:

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

出力パラメーターを指定して呼び出された関数は、カラー プロットを作成せずに結果を返すだけです。

関数は次のとおりです。

関数 [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
于 2013-11-03T12:01:18.600 に答える
0

問題のコードの改善

ループ (およびネストされた内積) はbsxfun、次のように行列の乗算で削除できます。

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')');

乗算ステップは、 で行われた AND と合計を実行しsum(inds1 .* inds2)ますが、密度行列をループすることはありません。 EDIT :singleの代わりにを使用するとdouble、実行時間はほぼ半分になりますが、回答を に変換するdoubleか、残りのコードに必要なものに変換してください。私のコンピュータでは、これには約0.5 秒かかります。

注: Withrot90(count_again/size(xy,1),2)には CDF があり、In にrot90(count_again_fix/size(xy,1),2)は PDF があります。

accumarray の使用

もう 1 つの方法はaccumarray、データをビン化した後に結合ヒストグラムを作成するために使用することです。

int_x、int_y、K1、などで始まるxy:

% 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

私のコンピューターでは、これには約0.01 秒かかります。

出力は、スパースからフルに変換された A.Donda のものと同じです ( full(H))。ただし、A.Donda が指摘したように、サイズが x であった OP コードではなく、寸法が x であることは正しいことK1です。K1count_again_fixK1+1K1+1

CDF を取得するにcumsumは、PDF の各軸に適用するだけでよいと思います。

于 2013-11-02T19:49:51.290 に答える