加速用于 FDR 估计的 MATLAB 代码

我有 2 个输入变量：

• p 值向量 (p) with N元素（未排序）
• and N x M具有通过随机排列获得的 p 值的矩阵 (pr) with M迭代。N相当大，10K到100K甚至更多。M假设是 100。

我正在估计每个元素的错误发现率（FDR）p表示如果当前 p 值（来自p）将是阈值。

我用 ARRAYFUN 编写了该函数，但是对于大 N（2min for N=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);.

嗯，诀窍确实是对向量进行排序。我对此表示感谢@EgonGeerardyn。另外，没有必要使用mean。您可以将所有内容除以M. When 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

and for N = 10000 ; M = 100 ;

已用时间为 39.730695 秒。
已用时间为 0.088870 秒。
2.220446049250313e-016