MATLab Bootstrap без цикла for

вчера я реализовал свой первый бутстрап в MATLab. (и да, я знаю, что циклы — это зло.):

%data is an mxn matrix where the data should be sampled per column but there 
can be a NaNs Elements
%from the array (a column of data) n values are sampled nReps times

function result = bootstrap_std(data, n, nReps,quantil)

result = zeros(1,size(data,2));

for i=1:size(data,2)
    bootstrap_data = zeros(n,nReps);
    values = find(~isnan(data(:,i)));
    if isempty(values)
       bootstrap_data(:,:) = NaN;
    else
        for k=1:nReps
            bootstrap_data(:,k) = datasample(data(values,i),n);
        end
    end

    stat = zeros(1,nReps);

    for k=1:nReps
        stat(k) = nanstd(bootstrap_data(:,k));
    end

    sort(stat);
    result(i) = quantile(stat,quantil);      
end
end

Как видите, эта версия работает по столбцам. Алгоритм делает то, что должен, но очень медленно, когда размер данных увеличивается. Теперь у меня вопрос: можно ли реализовать эту логику без использования циклов for? Моя проблема в том, что я не смог найти версию datasample, которая выполняет выборку по столбцам. Или есть более удобная функция?

Я рад любому намеку или идее, как я могу ускорить эту реализацию.

Спасибо и всего наилучшего!

Стефан


matlab statistics-bootstrap
person Stephan Claus    schedule 08.08.2014    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Узкие места в вашей реализации:

  1. Функция проводит много времени внутри nanstd, что не нужно, поскольку вы все равно исключаете значения NaN из своей выборки.
  2. Есть много функций, которые работают со столбцами, но вы тратите время на перебор столбцов и многократный их вызов.
  3. Вы делаете много вызовов datasample, что является относительно медленной функцией. Гораздо быстрее создать случайный вектор индексов с помощью randi и использовать его вместо этого.

Вот как я бы написал функцию (на самом деле я, вероятно, не стал бы оставлять столько комментариев и не использовал бы так много временных переменных, но я делаю это сейчас, чтобы вы могли видеть, каковы все этапы вычисления). ).

function result = bootstrap_std_new(data, n, nRep, quantil)

    result = zeros(1, size(data,2));

    for i = 1:size(data,2)
        isbad = isnan(data(:,i));                   %// Vector of NaN values
        if all(isbad)
            result(i) = NaN;
        else
            data0 = data(~isbad, i);                %// Temp copy of this column for indexing
            index = randi(size(data0,1), n, nRep);  %// Create the indexing vector
            bootstrapdata = data0(index);           %// Sample the data
            stdevs = std(bootstrapdata);            %// Stdev of sampled data
            result(i) = quantile(stdevs, quantil);  %// Find the correct quantile
        end
    end

end

Вот некоторые тайминги

>> data = randn(100,10);
>> data(randi(1000, 50, 1)) = NaN;
>> tic, bootstrap_std(data, 50, 1000, 0.5); toc
Elapsed time is 1.359529 seconds.
>> tic, bootstrap_std_new(data, 50, 1000, 0.5); toc
Elapsed time is 0.038558 seconds.

Таким образом, это дает вам примерно 35-кратное ускорение.

person Chris Taylor    schedule 08.08.2014
comment
Здравствуйте, я думаю, что nanstd не решит проблему, так как я бы включил значения NaN в свой образец данных, который не нужен. Или я что-то здесь упускаю? Спасибо и с наилучшими пожеланиями Стефан - person Stephan Claus; 08.08.2014
comment
А, кажется, я понимаю — вы хотите выкинуть все значения NaN из каждого столбца, чтобы даже не пробовать их? - person Chris Taylor; 08.08.2014
comment
точно. Хотя использовать nanstd здесь бесполезно. Я думаю, что это вводило в заблуждение. - person Stephan Claus; 08.08.2014
comment
Здравствуйте, Крис, спасибо за отредактированную версию - я думаю, вы правы насчет использования образца данных - я думаю, что это и внутренние циклы for являются основными узкими местами. Одно небольшое замечание: ваш пример слишком быстрый, поскольку нет значений NaN для фильтрации. Но даже с фильтрацией это ускорение в 10 раз. Спасибо за совет! - person Stephan Claus; 08.08.2014

arrow_upward
0
arrow_downward

Ваша основная проблема, по-видимому, заключается в том, что у вас могут быть разные числа/позиции NaN в каждом столбце, поэтому вы не можете работать с полной матрицей, если вы не согласны также с выборкой NaN. Однако некоторые внутренние циклы можно упростить.

for k=1:nReps
    bootstrap_data(:,k) = datasample(data(values,i),n);
end

Поскольку вы выполняете выборку с заменой, вы должны просто сделать:

bootstrap_data = datasample(data(values,i), n*nReps);
bootstrap_data = reshape(bootstrap_data, [n nReps]);

Также nanstd может работать с полной матрицей, поэтому нет необходимости зацикливаться:

stat = nanstd(bootstrap_data); % or nanstd(x,0,2) to change dimension

Также стоило бы просто просмотреть ваш код с profile, чтобы увидеть, где находятся узкие места.

person nkjt    schedule 08.08.2014
comment
Привет @nkjt, спасибо за ваш ответ - функция изменения формы хороша, но, как отметил Крис, образец данных очень медленный, поэтому лучше использовать версию randi. С уважением Стефан - person Stephan Claus; 08.08.2014
comment
Да, его ответ лучше. :) - person nkjt; 08.08.2014