2012-02-24 5 views
0

행렬의 각 요소가 요소가 선택 될 확률을 나타내는 매우 큰 정사각형 행렬 M (i, j)가 있다고 가정합니다. 가중치 임의 선택. 나는 (i, j) 인덱스에 의해 행렬로부터 n 개의 원소를 대체 할 필요가있다. 메인 루프의 반복마다 가중치가 변경됩니다. 때문에 for 루프에 MATLAB의 매우 큰 배열의 인덱스로 n 개의 가중치를 선택하십시오

for m = 1:M_size 
    xMean(m) = mean(M(:, m)); 
end 

[~, j_list] = histc(rand(n, 1), cumsum([0; xMean'./sum(xMean)'])); 
for c = 1:n 
    [~, i_list(c)] = ... 
     histc(rand(1, 1), cumsum([0;, M(:, j_list(c))./sum(M(:, j_list(c)))])); 
end 

을하지만이 또한 매우 긴 시간이 소요 다소 투박한 방법을 것 같다 :

현재, 나는 다음과 같은 것을 사용하고 있습니다. 더 효율적인 방법이 있습니까? 아마 내가 어떤 식으로 행렬을 벡터화한다면?

* 나는 통계 도구 상자에 미리

많은 감사를 액세스 할 수 없습니다 언급해야 편집.

+0

하기 위해 행렬을 벡터화 실험 후 실제로없이 (histc 한 번 cumsum, 그렇게하는 것은 더 오래 걸립니다 ... – Huggzorx

+0

당신이 M 행렬의 합을 알고 있는가 것으로 보인다 전화 그것을 요약)? – Pursuit

+0

확률이 임의로 평가되므로 합계가없는 합계를 알 수 없습니다. 이는 행렬의 일부 요소 값만 늘리는 보강 학습 알고리즘의 일부로 사용되기 때문입니다. 단순히 합계를 한 번 계산 한 다음 각 값에 새 값을 추가 할 수 있습니다. – Huggzorx

답변

1

randsample (docs)은 여기에 친구입니다. 당신은 적절한 크기를 얻을 수 M에 몇 전치을해야 할 수도 있습니다

selected_indexes = randsample(1:numel(M), n, true, M(:)); 
[sub_i, sub_j] = ind2sub(size(M), selected_indexes); 

: 나는 첨자의 뒤쪽으로 인덱스로 변환 다음 방법을 사용합니다.

+0

답장을 보내 주셔서 감사합니다.하지만 유감스럽게도 randsample은 통계 도구 상자에 포함되어 있습니다 (기관은 cheapskates입니다). – Huggzorx

+0

오, 죄송합니다. 나는 조심해야 했어. – AE426082

0
% M is ixj 
xMean = transpose(mean(M,1)); 
%xMean is jx1, so i hope n == j 
[~, j_list] = histc(rand(n, 1), cumsum([0; xMean./sum(xMean)])); 
% j_list is not used? but is j x 1 
cumsumvals = cumsum([zeros(1,jj);, M(:,j_list(1:n))./kron(sum(M(:,j_list(1:n))),ones(ii,1))],1),1) 
% cumsumvals is i+1 x j, so looks like it should work 
% but histc won't work with a matrix valued edge parameter 
% you'll need to look into hist3 for that 
for c = 1:n 
    [~, i_list(c)] = ... 
     histc(rand(1, 1), cumsumvals(:,c)); 
end 

그럼 더 가까워 지지만, 완전히 벡터화하려면 hist3이 필요합니다.

+0

죄송합니다. y_list는 내 코드에서 j_list로 읽어야합니다. – Huggzorx

0

사실 저는 이것을 비 벡터 라이징으로 해결할 것이라고 생각합니다. 즉, 사전 정의 된 배열과 간단한 조작 만 사용하여 높은 수준의 호출과 값 비싼 작업을 모두 제거하고 필수 요소까지 제거하십시오.

은 알고리즘의 핵심은 다음과 같습니다

  1. 은 선택 n은 0과 가중치의 합 사이의 임의의 숫자는,으로 정렬 가중치
  2. 의 합을 결정합니다.

  3. 수동으로 누누 루프를 구현하십시오. 그러나 모든 누적 합계를 저장하는 대신 cumsum이 현재의 난수보다 작은 곳에서 현재의 난수 이상으로 점프하는 인덱스를 저장하십시오. 코드에서

은 (타이밍 장비의 비트와 함께), 그 다음과 같습니다

tic 
for ixTiming = 1:1000 
    M = abs(randn(50)); 
    M_size = size(M, 2); 
    n = 8; 

    for m = 1:M_size 
     xMean(m) = mean(M(:, m)); 
    end 

    [~, j_list] = histc(rand(n, 1), cumsum([0; xMean'./sum(xMean)'])); 
    for c = 1:n 
     [~, i_list(c)] = ... 
      histc(rand(1, 1), cumsum([0;, M(:, j_list(c))./sum(M(:, j_list(c)))])); 
    end 
end 
toc; %1.10 sec on my computer 

:

tic 
for ixTiming = 1:1000 

    M = abs(randn(50)); 
    M_size = size(M, 2); 
    n = 8; 
    total = sum(M(:)); 

    randIndexes = sort(rand(n,1) * total); 

    list = zeros(n,1); 
    ixM = 1; 
    ixNextList = 1; 
    curSum = 0; 
    while ixNextList<=n && ixM<numel(M) 
     while curSum<randIndexes(ixNextList) && ixM<=numel(M) 
      curSum = curSum+M(ixM); 
      ixM = ixM + 1; 
     end 
     list(ixNextList) = ixM; 
     ixNextList = ixNextList+1; 
    end 
    [i_list, j_list] = ind2sub(size(M),list); 

end 
toc; %0.216 sec. on my computer 

이 원래의 질문에 코드이 비교

주의 사항 및 최적화.

  • 광범위하게 테스트하지 않았습니다. 무작위 수 연산은 적절한 임의 동작을하기가 어렵습니다. 많은 몬테카를로 세트에 대해 몇 가지 테스트 케이스를 실행하여 동작이 예상대로 작동하는지 확인하십시오. 특히 off-by-one 유형 오류가 있는지 조심하십시오.

  • 프로필을 참조한 다음 느린 단계에서 추가 개선 사항을 찾으십시오. 어떤 가능성. 당신이 M을 변경하면

    • total 값을 유지, 그래서 당신은 돈, t는 다시 계산해야합니다.

    • randIndexes의 가장 낮은 값과 가장 높은 값을 0total에 대해 확인하십시오. 이 경우 randIndexes(1) is larger than total-randIndexes (종료) , then increment ixM from numel (M) numel (M)`입니다.