2013-07-02 6 views
3

매트랩 1000 5 × 5 행렬의 공분산 계산 :I 이런 1000 5 × 5 행렬 (</strong><strong>내지 Xm)가

enter image description here

각각 $ (x_ij) 미터을 $로부터 인출 점 추정치 분포. 나는 각각의 $ x {ij} $의 공분산 cov을 계산하고 싶습니다. 여기서 i = 1..n, j = 1..n은 빨간색 화살표 방향입니다.

예를 들어, $ X_m $의 분산은 'var (X, 0,3)'이며, 분산의 5x5 행렬을 제공합니다. 같은 방식으로 공분산을 계산할 수 있습니까? 대답에서

시도가

지금까지이 작업을 완료했습니다 당신이 (명령 창에서 edit cov) cov 보면

for m=1:1000 
Xm_new(m,:)=reshape(Xm(:,:,m)',25,1); 
end 

cov(Xm_new) 
spy(Xm_new) gives me this unusual looking sparse matrix: 

enter image description here

답변

5

당신은 이유를 볼 수 있습니다 다차원 배열을 지원하지 않습니다. 이것은 입력 행렬의 전치와 행렬 곱셈을 수행합니다 : xc' * xc. 두 작업 모두 다차원 배열을 지원하지 않으며 함수를 작성한 사람이 일반화 작업을하지 않기로 결정한 것 같습니다. (여전히 Mathworks 및 make a feature request에 문의하는 것이 좋습니다.) 우리가 cov에서 기본 코드를 가지고 몇 가지 가정을 경우

귀하의 경우에는

는, 우리는 공분산 기능 M-파일에게 지원하는 3 차원 배열을 작성할 수 있습니다

function x = cov3d(x) 
% Based on Matlab's cov, version 5.16.4.10 

[m,n,p] = size(x); 
if m == 1 
    x = zeros(n,n,p,class(x)); 
else 
    x = bsxfun(@minus,x,sum(x,1)/m); 
    for i = 1:p 
     xi = x(:,:,i); 
     x(:,:,i) = xi'*xi; 
    end 
    x = x/(m-1); 
end 

참고가이 간단한 코드 x은 3 차원을 따라 겹쳐진 일련의 2D 행렬이라고 가정합니다. 그리고 정규화 플래그는 0이고 기본값은 cov입니다. 약간의 작업을 통해 var과 같은 여러 차원으로 확장 될 수 있습니다. 내 타이밍에서는 for 루프에서 cov(x(:,:,i))을 호출하는 함수보다 10 배 이상 빠릅니다.

예, 저는 for 루프를 사용했습니다. faster ways to do this 일 수도 있고 아닐 수도 있지만,이 경우 for 루프는 faster than most schemes이 될 것입니다. 특히 어레이의 크기가 선험적으로 알려지지 않은 경우에는 더욱 그렇습니다.

+0

아 매우 좋습니다! 고맙습니다! 그러나 각 변수 x_ij의 1000 회의 관측 (3 차원에서 서로 겹쳐진 온 탑)에 기반한 5x5 공분산 행렬을 찾고 있습니다. 예를 들어 Xm (1,1, :)은 같은 변수에 대한 관찰 값입니다.이 함수가 무엇을하고 있습니까? – HCAI

+0

5-by-5 ​​입력에 대해 'cov' 함수를 1000 번 호출하고 5x5 출력을 저장하는 것과 같습니다 :'Xm_out (:, :, 1) = cov (Xm (:, :, 1)), Xm_out (:, :, 2) = cov (Xm (:, :, 2))', ...,'Xm_out (:, :, 999))','Xm_out (:, :, 1000) = cov (Xm (:, :, 1000))'로 시작한다. – horchler

+0

아 아 ...각 x_ij에 대한 공분산을 찾으려고합니다. 즉, x12는 열 벡터 x1과 x2 사이의 공분산이고, 빨간색 화살표의 방향입니다. 올바른 질문을하지 않았습니까? – HCAI

2

또한 아래의 대답은 xi=x(:,:,i)

function xy = cov3d(x) 

[m,n,p] = size(x); 
if m == 1 
    x = zeros(n,n,p,class(x)); 
else 
    xc = bsxfun(@minus,x,sum(x,1)/m); 
    for i = 1:p 
     xci = xc(:,:,i); 
     xy(:,:,i) = xci'*xci; 
    end 
    xy = xy/(m-1); 
end 

내 대답은 horchler 매우 유사하다 직사각형 행렬 작동하지만 horchler의 코드는 (그 차원 xi'*xi 차원에서 다른) 직사각형 행렬 xi 작동하지 않습니다.

관련 문제