从计数向量中随机选取元素

Question

我目前正试图通过改变算法来优化一些MATLAB/Octave代码，但不知道如何处理一些随机性。假设我有一个整数向量V，每个元素代表一些东西的计数，在我的例子中是光子。现在，我想随机选取一些“东西”，创建一个大小相同的新向量，但要调整计数。

以下是我目前的做法：

function W = photonfilter(V, eff)
% W = photonfilter(V, eff)
% Randomly takes photons from V according to the given efficiency.
%
% Args:
%  V: Input vector containing the number of emitted photons in each
%     timeslot (one element is one timeslot). The elements are rounded
%     to integers before processing.
%  eff: Filter efficiency. On the average, every 1/eff photon will be
%       taken. This value must be in the range 0 < eff <= 1.
%  W: Output row vector with the same length as V and containing the number
%     of received photons in each timeslot.
%
% WARNING: This function operates on a photon-by-photon basis in that it
% constructs a vector with one element per photon. The storage requirements
% therefore directly depend on sum(V), not only on the length of V.

% Round V and make it flat.
Ntot = length(V);
V = round(V);
V = V(:);

% Initialize the photon-based vector, so that each element contains
% the original index of the photon.
idxV = zeros(1, sum(V), 'uint32');
iout = 1;
for i = 1:Ntot
  N = V(i);
  idxV(iout:iout+N-1) = i;
  iout = iout + N;
end;

% Take random photons.
idxV = idxV(randperm(length(idxV)));
idxV = idxV(1:round(length(idxV)*eff));

% Generate the output vector by placing the remaining photons back
% into their timeslots.
[W, trash] = hist(idxV, 1:Ntot);

这是上面描述的一个相当简单的实现。但是它有一个明显的性能缺陷：函数创建一个向量（idxV），每个光子包含一个元素。所以如果我的V只有1000个元素，但是每个元素的平均计数是10000，内部向量将有1000万元素，使得函数缓慢而沉重。

我现在想要实现的不是直接优化这个代码，而是使用其他类型的算法，可以立即计算新的计数，而不给每个光子某种“身份”。这一定是有可能的，但我就是不知道怎么做。

要求：

• 输出向量W的元素数必须与输入向量V的元素数相同。
• W（i）必须是一个整数，且以0为界

有没有办法实现这样的东西？理想的解决方案是只使用一个随机向量，然后使用一些概率和舍入的技巧，但到目前为止，我还没有在这方面取得任何成功。

谢谢向你问好，菲利普

Answer 1

根据Alexander Solovets的回答，代码现在是这样的：

function W = photonfilter(V, eff, impl=1)

Ntot = length(V);
V = V(:);

if impl == 0
  % Original "straightforward" solution.
  V = round(V);
  idxV = zeros(1, sum(V), 'uint32');
  iout = 1;
  for i = 1:Ntot
    N = V(i);
    idxV(iout:iout+N-1) = i;
    iout = iout + N;
  end;
  idxV = idxV(randperm(length(idxV)));
  idxV = idxV(1:round(length(idxV)*eff));
  [W, trash] = hist(idxV, 1:Ntot);

else
  % Monte Carlo approach.
  Nphot = sum(V);
  P = cumsum(V / Nphot);
  W = hist(lookup(P, rand(1, round(Nphot * eff))), 0:Ntot-1);

end;

只要eff不太接近1（eff=1，原始解产生W=V，而蒙特卡罗方法仍具有一定的随机性，从而违反了上限约束），结果就相当具有可比性。

在交互式倍频程外壳中测试：

octave:1> T=linspace(0,10*pi,10000);
octave:2> V=100*(1+sin(T));
octave:3> W1=photonfilter(V, 0.1, 0);
octave:4> W2=photonfilter(V, 0.1, 1);
octave:5> plot(T,V,T,W1,T,W2);
octave:6> legend('V','Random picking','Monte Carlo')
octave:7> sum(W1)
ans =  100000
octave:8> sum(W2)
ans =  100000

Answer 2

你用来计算W的方法叫做蒙特卡罗方法。而且确实可以进行一些优化。一旦这样做不是计算光子指数，让我们想象一组箱子。每个箱子都有一定的概率，所有箱子的概率之和加起来等于1。我们将段[0，1]分成若干部分，其长度与料仓的概率成正比。现在，对于我们生成的[0,1]中的每个随机数，我们可以快速找到它所属的箱子。最后，我们计算箱子中的数字以获得最终结果。下面的代码说明了这个想法。

% Population size (number of photons).
N = 1000000;
% Sample size, size of V and W as well.
% For convenience of plotting, V and W are of the same size, but
% the algorithm doesn't enforce this constraint.
M = 10000;
% Number of Monte Carlo iterations, greater numbers give better quality.
K = 100000;

% Generate population of counts, use gaussian distribution to test the method.
% If implemented correctly histograms should have the same shape eventually.
V = hist(randn(1, N), M);
P = cumsum(V / sum(V));
% For every generated random value find its bin and then count the bins.
% Finally we normalize counts by the ration of N / K.
W = hist(lookup(P, rand(1, K)), M) * N / K;
% Compare distribution plots, they should be the same.
hold on;
plot(W, '+r');
plot(V, '*b');
pause