加速MATLAB程序进行FDR估计

Question

我有两个输入变量：

• 具有N个元素的p值向量(p) (未排序)
• N×M矩阵，p-值由随机置换(pr，)和M次迭代得到.N相当大，10K到100 K或更多。I‘让我们说100。

我正在估计p的每个元素的错误发现率(FDR)，它表示如果当前p值(来自p)为阈值，随机排列将通过多少p值。

我用ARRAYFUN编写了这个函数，但是大N(N=20K为2分钟)需要大量时间，可以与for -循环相媲美。

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少。

编辑(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

首先，使用轮廓仪来分析这一点。在试图提高性能时，分析应该始终是第一步。我们都可以猜测是什么导致了性能下降，但是确保并专注于正确部分的唯一方法是检查分析器报告。

我没有在您的代码上运行分析器，因为我不想为此生成测试数据；但是，对于正在执行的工作，我有一些想法。在您的函数mean(sum(pr<=x))./sum(p<=x)中，您重复地在p<=x上求和。总之，一个调用包括N比较和N-1求和。因此，在N中，当计算出p的所有N值时，您的行为都是二次的。

如果您逐步执行排序版本的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;是主要的资源占优势。通过使用p的排序特性，可以处理类似于pr的问题。

我用您的代码测试了我的代码(结果相同，时间消耗也一样)。对于N=20e3; M=100，我有63秒的时间运行你的代码，43秒的时间在我的主计算机上运行我的代码(Matlab2011 a在64位Arch上，8 GiB内存，核心i7 860)。对于较小的M值，增益更大。但这一增加部分是由于并行化。

edit2:显然，我得到了与Andrey非常相似的结果，如果我采用同样的方法，我的结果会非常相似。

然而，我意识到有一些内置的函数或多或少地满足了你的需要，即非常类似于确定经验累积密度函数。这可以通过构造直方图来实现：

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毫秒。Andrey的参数需要104毫秒，所以在我的计算机上，它的速度要慢一些，但是我认为这段代码比复杂的循环要容易读得多(就像我们两个例子中的情况一样)。

Answer 3

在这个问题中我和Andrey进行了讨论之后，这个非常晚的答案只是为了向Andrey证明矢量化的解决方案仍然比JIT‘’ed循环更快，它们有时并不容易找到。

如果“任择议定书”认为这一答案不恰当，我非常愿意删除这一答复。

现在，关于业务，这里是最初的arrayfun，由Andrey编写的循环版本，以及Egon的矢量化版本：

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倍。