Matlab中三點相關函數的矢量化計算？

Question

我想通過合適的矢量元素的直方圖（長度爲system_size的num_samples採樣）和相應的聚類函數T2，T3來計算矢量樣本的二點和三點相關函數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（垃圾箱，垃圾箱......）或交叉表（箱，桶），幾乎可以做我想做的，這是尋找向量的元素的相關事件，但它不」噸。

---

實施例：

如果最外層循環內我的輸入是

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 

Answer 1

「R2 and R3`很容易：

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是一個行向量。

因爲代碼中有些東西還不清楚，所以您可能需要稍微玩一下，但這應該與您最終需要的東西非常接近。

EDIT澄清後：

對於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應該看起來仍然類似於我以前的答案。