性能 – 加速用于FDR估计的MATLAB代码

我有2个输入变量：

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

嗯,诀窍确实是对矢量进行排序.我赞扬了@EgonGeerardyn.此外,没有必要使用平均值.之后你可以将所有内容除以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：

Elapsed time is 0.898689 seconds.
Elapsed time is 0.007697 seconds.
2.220446049250313e-016

并且对于N = 10000; M = 100;

Elapsed time is 39.730695 seconds. Elapsed time is 0.088870 seconds. 2.220446049250313e-016