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