0

ベクトルの要素 (長さ system_size の num_samples サンプル) と対応するクラスター関数 T2、T3 の適切なヒストグラム処理によって、ベクトルのサンプルの 2 点と 3 点の相関関数 R2、R3 を計算したいと思います。簡単にするために、均一なビン全体のヒストグラムを検討しています。

次のコードをベクトル化および/または高速化する良い方法は何ですか?

n = length(mesh);
R2 = zeros(n, n);
R3 = zeros(n, n, n);
for sample_id=1:num_samples 
    s = samples(:, sample_id);
    d = mesh(2) - mesh(1);
    % Which bin does the ith sample s belong to?
    bins = ceil((s - mesh(1))/d);

    % Compute two-point correlation function
    for i = 1:system_size
        for j = 1:system_size
            if i ~= j
                R2(bins(i), bins(j))=R2(bins(i), bins(j))+1;
            end
        end
    end

    % Compute three-point correlation function
    for i = 1:system_size
        for j = 1:system_size
            if i ~= j
                for k = 1:system_size
                    if k ~= j && k ~= i
                        R3(bins(i), bins(j), bins(k))=R3(bins(i), bins(j), bins(k))+1;
                        T3(x1, x2, x3) = R3(x1,x2,x3)-R1(x1)*R2(x2,x3)-R1(x2)*R2(x1,x3)...
                             -R1(x3)*R2(x1,x2)+2*R1(x1)*R1(x2)*R1(x3);
                    end
                end
            end
        end
    end
end
R2 = R2/sum(R2(:));
R3 = R3/sum(R3(:));

T3 = zeros(n, n, n);
% Compute three-point cluster function
for i = 1:n
    for j = 1:n
        if i ~= j
            for k = 1:n
                if k ~= j && k ~= i
                    T3(x1, x2, x3) = R3(x1,x2,x3)-R1(x1)*R2(x2,x3)-R1(x2)*R2(x1,x3)...
                         -R1(x3)*R2(x1,x2)+2*R1(x1)*R1(x2)*R1(x3);
                end
            end
        end
    end
end

素朴に hist3(bins, bins...) または crosstab(bins, bins) は、ベクトルの要素の相関する出現を探すという、私が望むことをほとんど行うと思っていましたが、そうではありません。

例：

最も外側のループ内の入力が

s = [1.2 3.1 4.6 4.7 5.1]
mesh = 0:0.5:6

量子化されたデータは

bins = [3 7 10 10 11]

R2 は

>> R2

R2 =

     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     0     0     2     1     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0     2     1     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     2     0     0     0     2     0     0     2     2     0
     0     0     1     0     0     0     1     0     0     2     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
4

1 に答える 1

0

「R2 andR3」は簡単です:

R2 = R2 + 1 - diag(ones(size(R2, 1), 1); % you can replace the loop with this
eye3 = zeros(n, n, n);
eye3(linspace(1, numel(eye3), n)) = 1;
R3 = R3 + 1 - eye3; % can move R3 computation outside the loop

の場合T3:

temp = repmat(R2, [1 1 n]).*permute(repmat(R1, [n, 1, n]), [1, 3, 2]);
T3 = R3 - temp - permute(temp, [2 3 1]) - permute(temp, [3 1 2]);
temp2 = repmat(R1'*R1, [1 1 n]).*permute(repmat(R1, [n, 1, n]), [1, 3, 2]);
T3 = T3 + temp2;

R1は行ベクトルであると仮定します。

コードからはまだ不明な点がいくつかあるため、これを少しいじる必要があるかもしれませんが、これは最終的に必要になるものにかなり近いはずです。

明確化後に編集：

の場合R2:

ubins = unique(bins);
bincounts = histc(bins, ubins);
for i=1:max(bincounts)
    indices = find(bincounts == i);
    R2(indices, indices) = R2(indices, indices) + i
end

これは、大きなベクトルと配列の場合にのみ役立ちます。実際には、行列全体ではなく、行列のチャンクの計算をベクトル化しています ( で繰り返しが発生する可能性があるためbins)。

に似たようなものを書くことができますR3。T3まだ私の以前の答えに似ているはずです。

于 2012-08-10T06:37:46.670 に答える