Можно ли векторизовать этот цикл случайной выборки для оптимизации?

Question

Я пытаюсь найти способ векторизации цикла FOR в коде ниже:

h=load('water-column');                 % load file 
perm=5;                         % make 10000 permutation
n_1=5;                            % number of random sample
dm_ale=zeros(1,perm);               % create a vector 
sz=length(h);                     % count size of matrix data    
for k=1:perm                      % making loop for the permutation
    i_1=randsample(sz,n_1);      
    x_3=h(i_1);            
    x_4=h(setdiff(1:sz,i_1));    
    dm_ale(k)=abs(mean(x_3)-mean(x_4)); % calculate difference of mean for each permutation
end

Что касается ввода файла, у меня есть что-то вроде этого (просто пример, реальный файл содержит больше данных):

   3792.615000000000
   3792.625000000000
   3792.634000000000
   3792.640000000000
   3792.647000000000
   3792.654000000000
   3792.662000000000
   3792.668000000000
   3792.673000000000

Я не могу понять, куда я могу поместить приращение в векторизованном выражении. Можно ли векторизовать это?

Когда код, предложенный Крисом Луенго (извините, я не могу понять, как пометить пользователя), натолкнулся на ошибку:

error: randsample: The input k must be a non-negative integer. Sampling without replacement needs k <= n.
error: called from
    randsample at line 46 column 5
    random_sampling at line 8 column 5

где random_sampling - название кода.

Изначально мне нужно иметь perm = 10000 (чтобы иметь надежный тест случайной выборки) и n_1 = 600 (количество населения, необходимое для того, чтобы мой тест мог работать). Код выше, кажется, не работает, даже если я подчиняюсь условию: n_1^2 << <code>perm. Я предполагаю, что ошибка связана с n_1, который все еще достаточно велик по отношению к perm. Любое другое руководство? Я думаю об увеличении perm.

Answer 1

Вы не можете использовать randsample для одновременной генерации нескольких случайных выборок (или так кажется из чтения документации).Если h достаточно велико, а perm и n_1 достаточно мало (sz >> perm*n_1), то вы можете создать случайную выборку с элементами perm*n_1, а затем разделить ее на perm комплекты.Это может быть примерно нормально, но не совсем то, что вы делаете сейчас.

Ваш код будет выглядеть так (с использованием упрощения, предложенного Джеффри Брентом в комментарии ):

h = load('col-deau');
perm = 5;
n_1 = 5;
sz = numel(h);  % numel is always better than length if you use one index h(i_1) rather than two h(i_1,1)
sum_h = sum(h)
i_1 = randsample(sz, n_1 * perm);
i_1 = reshape(i_1, n_1, perm);
x_3 = h(i_1);                    % x_3 has the same size as i_1
x_3 = sum(x_3, 1);               % sum over columns, x_3 has perm elements now
x_4 = sum_h - x_3;
dm_ale = abs(x_3 / n_1 - x_4 / (sz-n_1));

Если perm также большой (как указано в комментарии), но n_1 все еще мал, вы можете аппроксимировать это, используя случайную выборку с заменой (с маленькой n_1вероятность того, что у вас есть повторяющийся элемент в одном наборе, мала):

i_1 = randsample(sz, n_1 * perm, true);