Ускорение MATLAB код для оценки FDR

Question

У меня 2 входные переменные:

• вектор р-значений (р) с N элементов (несортированных)
• и Н х M матрица с p-значениями, полученными случайными перестановками (pr) с M итераций. N довольно большой, от 10 до 100 тыс. Или более. M скажет 100.

Я оценивающее ошибочный Discovery Rate (FDR) для каждого элемента, представляющего p сколько р-значение от случайных перестановок будут проходить, если текущий р-значение (от p) будет порогом.

Я написал функцию с ARRAYFUN, но это занимает много времени для больших N (2 мин для N = 20K), что сравнимо с для цикла.

function pfdr = fdr_from_random_permutations(p, pr) 
%# ... skipping arguments checks 
pfdr = arrayfun(@(x) mean(sum(pr<=x))./sum(p<=x), p); 

Любые идеи, как сделать это быстрее?

Комментарии о статистических проблемах здесь также приветствуются.

Данные испытаний могут быть сгенерированы как p = rand(N,1); pr = rand(N,M);.

Answer 1

Ну, трюк действительно сортировал векторы. Я отдаю должное @EgonGeerardyn за это. Кроме того, нет необходимости использовать mean. Вы можете просто все разделить на M. Когда p сортируется, поиск количества значений, которое меньше текущего x, является просто индексом работы. pr - более интересный случай - я использовал индекс бега с именем place, чтобы узнать, сколько элементов меньше x.

Edit (2): Вот это самый быстрый вариант я придумал:

function Speedup2() 
    N = 10000/4 ; 
    M = 100/4 ; 
    p = rand(N,1); pr = rand(N,M); 

    tic 
    pfdr = arrayfun(@(x) mean(sum(pr<=x))./sum(p<=x), p); 
    toc 

    tic 
    out = zeros(numel(p),1); 
    [p,sortIndex] = sort(p); 
    pr = sort(pr(:)); 
    pr(end+1) = Inf; 
    place = 1; 
    N = numel(pr); 
    for i=1:numel(p) 
     x = p(i); 
     while pr(place)<=x 
      place = place+1; 
     end 
     exp1a = place-1; 
     exp2 = i; 
     out(i) = exp1a/exp2; 
    end 
    out(sortIndex) = out/ M; 
    toc 
    disp(max(abs(pfdr-out))); 

end 

И результаты тестов для N = 10000/4 ; M = 100/4:

истекшее время 0.898689 секунд.
Истекшее время составляет 0.007697 секунд.
2.220446049250313e-016

и N = 10000 ; M = 100;

Истекшее время 39.730695 секунд.
Истекшее время: 0.088870 секунд.
2.220446049250313e-016

Answer 2

Прежде всего, tr проанализировать это, используя profiler. Профилирование должно ВСЕГДА быть первым шагом при попытке повысить производительность. Мы все можем догадываться о том, что приводит к снижению производительности, но единственный способ быть уверенным и сосредоточиться на правильной части - проверить отчет профилировщика.

Я не запускал профилировщик вашего кода, так как я не хочу генерировать тестовые данные для этого; но у меня есть некоторые идеи о том, какая работа выполняется напрасно. В вашей функции mean(sum(pr<=x))./sum(p<=x), вы повторно суммируете более p<=x. Всего один звонок включает в себя N сравнения и N-1 суммирования. Итак, для обоих, у вас есть поведение, которое квадратично в N, когда вычисляются все значения N значений p.

Если вы выполните отсортированную версию p, вам нужно меньше вычислений и сравнений, так как вы можете отслеживать текущую сумму (то есть поведение, которое является линейным в N). Я предполагаю, что аналогичный метод может быть применен к другой части расчета.

редактировать: Реализация моей идеи, выраженная выше:

function pfdr = fdr(p,pr) 
[N, M] = size(pr); 
[p, idxP] = sort(p); 
[pr] = sort(pr(:)); 

pfdr = NaN(N,1); 

parfor iP = 1:N 
    x = p(iP); 
    m = sum(pr<=x)/M; 
    pfdr(iP) = m/iP; 
end 

pfdr(idxP) = pfdr; 

Если у вас есть доступ к параллельной вычислительной инструментарии, цикл parfor позволит вам получить некоторое представление. Я использовал две основные идеи: mean(sum(pr<=x)) фактически равен sum(pr(:)<=x)/M.С другой стороны, поскольку p сортируется, это позволяет просто взять индекс как количество элементов (в предположении, что каждый элемент уникален, иначе вам придется работать с unique, чтобы выполнить полный тщательный анализ).

Как вы уже должны хорошо знать, запустив профайлер самостоятельно, линия m = sum(pr<=x)/M; является основным ресурсом hog. Это можно решить аналогично p, используя отсортированную природу pr.

Я проверил свой код (как для идентичных результатов, так и для потребления времени) против вашего. Для N=20e3; M=100 я получаю около 63 секунд для запуска вашего кода и 43 секунды для запуска моего на моем основном компьютере (MATLAB 2011a на 64-битной Arch Linux, 8 GiB RAM, Core i7 860). При меньших значениях M коэффициент усиления больше. Но этот выигрыш частично объясняется распараллеливанием.

edit2: По-видимому, я пришел к очень похожим результатам, как Андрей, мой результат был бы очень схожим, если бы я придерживался того же подхода.

Однако я понял, что есть некоторые встроенные функции, которые делают больше или меньше того, что вам нужно, то есть очень похоже на определение эмпирической функции кумулятивной плотности. И это можно сделать с помощью построения гистограммы:

function pfdr = fdr(p,pr) 
[N, M] = size(pr); 
[p, idxP] = sort(p); 

count = histc(pr(:), [0; p]); 
count = cumsum(count(1:N)); 

pfdr = count./(1:N).'; 

pfdr(idxP) = pfdr/M; 

Для того же M и N, как указано выше, этот код занимает 228 миллисекунд на моем компьютере. Для параметров Андрея требуется 104 миллисекунды, поэтому на моем компьютере это происходит немного медленнее, но я думаю, что этот код гораздо читабельнее, чем запутанный для циклов (как это было в обоих наших примерах).

Answer 3

После обсуждения между мной и Андреем в this question, это очень поздно ответ, чтобы доказать, что Андрей векторизованные решения еще быстрее, чем JIT'ed петли, иногда они просто не так легко найти.

Я более чем охотно удаляю этот ответ, если его ОП считают неуместным.

Теперь к делу, вот оригинальный arrayfun, петельные версии Андрея, и векторизованную версия Эгона:

function test 

    clc 

    N = 10000/4 ; 
    M = 100/4 ; 
    p = rand(N,1); 
    pr = rand(N,M); 

    %% first option 

    tic 

    pfdr = arrayfun(@(x) mean(sum(pr<=x))./sum(p<=x), p); 

    toc 


    %% second option 

    tic 

    out = zeros(numel(p),1); 
    [p2,sortIndex] = sort(p); 
    pr2 = sort(pr(:)); 
    pr2(end+1) = Inf; 
    place = 1;  
    for i=1:numel(p2) 
     x = p2(i); 
     while pr2(place)<=x 
      place = place+1; 
     end 
     exp1a = place-1; 
     exp2 = i; 
     out(i) = exp1a/exp2; 
    end 
    out(sortIndex) = out/ M; 

    toc 

    %% third option 

    tic 
    [p2,sortIndex] = sort(p); 

    count = histc(pr2(:), [0; p2]); 
    count = cumsum(count(1:N)); 

    out = count./(1:N).'; 

    out(sortIndex) = out/M; 

    toc 

end 

Результаты на моем ноутбуке:

Elapsed time is 0.916196 seconds. 
Elapsed time is 0.011429 seconds. 
Elapsed time is 0.007328 seconds. 

и для N=1000; M = 100;:

Elapsed time is 38.082718 seconds. 
Elapsed time is 0.127052 seconds. 
Elapsed time is 0.042686 seconds. 

Так: векторизованный в 2-3 раза быстрее.