2012-09-15 49 views
1

默認情況下,所有用於計算相關性或協方差的內置函數返回一個矩陣。我試圖編寫一個有效的函數來計算種子區域和其他各個區域之間的相關性,但我不需要其他區域之間的相關性。我假設計算全相關矩陣因此效率不高。MATLAB:與種子區域的相關

我可以改爲計算每個區域和種子區域之間的相關矩陣,選擇一個非對角點並存儲它,但我覺得在這種情況下循環也是低效的。

更具體地說,我的三維空間中的每個點都有一個時間維度。我試圖計算給定點與給定半徑內所有空間點之間的平均相關性。我想重複這個過程數十萬次,對於許多不同的半徑長度,等等,所以我希望這樣做盡可能高效。

那麼,計算單個向量和其他幾個向量之間的相關性的最佳方法是什麼,而不計算我將忽略的相關性?

謝謝 克里斯

編輯:這裏是我的代碼現在...

function [corrMap] = TIME_meanCorrMap(A,radius) 
% Even though the variable is "radius", we work with cubes for simplicity... 
% So, the radius is the distance (in voxels) from the center of the cube an edge. 
denom = ((radius*2)^3)-1; 
dim = size(A); 
corrMap = zeros(dim(1:3)); 
for x = radius+1:dim(1)-radius 
    rx = [x-radius : x+radius]; 
    for y = radius+1:dim(2)-radius 
     ry = [y-radius : y+radius]; 
     for z = radius+1:dim(3)-radius 
      rz = [z-radius : z+radius]; 
      corrCoefs = zeros(1,denom); 
      seed = A(x,y,z,:); 
      i=0; 
      for xx = rx 
       for yy = ry 
        for zz = rz 
         if ~all([x y z] == [xx yy zz]) 
          i = i + 1; 
          temp = corrcoef(seed,A(xx,yy,zz,:)); 
          corrCoeffs(i) = temp(1,2); 
         end 
        end 
       end 
      end 
      corrMap = mean(corrCoeffs); 
     end 
    end 
end 

編輯:這裏有一些更多的時間來補充公認的答案。 使用bsxfun()做歸一化,和矩陣乘法來計算相關性:

tic; for i=1:10000                 
    x=rand(100); 
    xz = bsxfun(@rdivide,bsxfun(@minus,x,mean(x)),std(x)); 
    cc = xz(:,2:end)' * xz(:,1) ./ 99; 
end; toc 
Elapsed time is 6.928251 seconds. 

使用zscore()來歸一化,矩陣乘法來計算相關性:

tic; for i=1:10000          
    x=rand(100);           
    xz = zscore(x);          
    cc = xz(:,2:end)' * xz(:,1) ./ 99; 
end; toc 
Elapsed time is 7.040677 seconds. 

使用bsxfun()來歸一化,和corr()來計算相關性。

tic; for i=1:10000          
    x=rand(100); 
    xz = bsxfun(@rdivide,bsxfun(@minus,x,mean(x)),std(x)); 
    cc = corr(x(:,1),x(:,2:end)); 
end; toc 
Elapsed time is 11.385707 seconds. 

回答

4

肯定可以改進您當前使用的for循環。如果有足夠的RAM,則可以使用矩陣乘法並行化相關計算。但是,它需要你將你的4維數據矩陣A打開成不同的形狀。很可能你正在處理三維立體fMRI數據,在這種情況下,你將不得不從[x y z時間]矩陣重塑爲[索引時間]矩陣。我會假設你可以處理重塑。一旦你有你的seed timecourse [時間1]和你的target timecourses [NumTargets的時間]準備就緒,你可以執行一些更有效的計算。

有效計算所需相關性的快速方法是使用MATLAB中的corr函數。該函數將接受2個矩陣參數,並且它將非常有效地計算參數1的列與參數2的列之間的所有成對相關性,例如,

T = 200; %time samples 
N = 20; %number of other voxels 

seed = randn(T,1);  %data from seed voxel 
targets = randn(T,N); %data from target voxels 

%here is the for loop method 
tic 
for n = 1:N 
    tmp = corrcoef(seed, targets(:,n)); 
    tmpcc = tmp(1,2); 
end 
looptime = toc; 

%here is the parallel method 
tic 
cc = corr(seed, targets); 
matrixtime = toc; 

在我的機器,在corr並行操作比正比於T * N的因子環路方法更快。

如果你願意自己完成底層矩陣操作,那麼有可能比corr函數快一點,並且在任何情況下都值得知道它們是什麼。兩個向量之間的關係基本上是一個歸一化點的產品,所以在使用的約定上面可以計算按以下方式

zseed = zscore(seed); %normalize the seed timecourse by z-scoring 
ztargets= zscore(targets); %normalize the target timecourses by z-scoring 
ztargets = ztargets';  %flip columns and rows for convenience 
cc2 = ztargets*zseed./(T-1); %compute many dot products with one matrix multiplication 

上面的代碼的相關性基本上是什麼corr功能會做這就是爲什麼它是非常比循環更快。請注意,大多數操作時間位於zscore操作中,如果使用bsxfun命令有效計算zscore,則可以提高corr函數的性能。就目前而言,我希望這給你一些指導,說明如何計算種子時間過程與許多目標時間過程之間的相關性,而不必分別進行循環和計算。

+0

這非常有幫助,謝謝!事實上,你從我的描述和代碼中推導出我的數據的性質是一個很好的添加觸摸。 –

+0

嗯,我過去曾經研究過類似的問題,所以猜測並不太可能。樂意效勞! – cjh