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