加速用于FDR估计的MATLAB代码

Question

我有2个输入变量:

• 带有N个元素的p值(p)向量(未排序)
• 和N × M矩阵,其具有通过随机排列(pr)以M次迭代获得的p值.N非常大,10K到100K或更多.M让我们说100.

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

我写与ARRAYFUN的功能,但它需要大量的时间用于大N(2 分钟为Ñ = 20K),媲美的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

首先,使用分析器进行分析.在尝试提高性能时,分析应始终是第一步.我们都可以猜测导致性能下降的原因,但确定并专注于正确部分的唯一方法是检查分析器报告.

我没有在你的代码上运行探查器,因为我不想生成测试数据这样做; 但我对于哪些工作是徒劳无功而有所了解.在你的功能中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;是主要的资源.这可以p通过利用排序的性质来类似地解决pr.

我测试了我的代码(相同的结果和时间消耗)与你的代码.因为N=20e3; M=100,我运行你的代码大约需要63秒,在主计算机上运行我需要43秒(在64位Arch Linux上的MATLAB 2011a,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毫秒.对于Andrey的参数,它需要104毫秒,所以在我的计算机上它变得有点慢,但我认为这个代码比复杂的循环更具可读性(就像我们的例子中的情况一样).

Answer 2

嗯,诀窍确实是对矢量进行排序.我赞扬了@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