Compound_Poisson_Bench.m

%% Monte Carlo speed test - Compound Poisson
% Benchmarks several different ways of generating Compound Poisson random 
% variables by Monte Carlo simulation.
%
% Compound Poisson random variables are of the form Y=X1+...+XN where the
% Xj are iid random variables and N is random with the Poisson
% distribution.
%
% Here we benchmark for the particular example where N has Poisson
% parameter 10 and the Xj are Lognormal with parameter mu=0 and sigma=1.
%
% The published version was run on a 2.6GHz dual core (MATLAB 2009b,
% Windows XP SP3).  Actual run times could vary by about 10% between runs.
%
% Requires POISSRND (Statistics Toolbox) and RANDRAW (FEX #7309).
%
% Author: Ben Petschel 10/12/2009

%%
n=1e6; % number of samples
lam=10; % Poisson lambda
mu=0; % lognormal mu parameter
sig=1; % lognormal sigma parameter


%% Tests using POISSRND (the stats toolbox Poisson sampler)
fn=@(m)poissrnd(lam,m,1); % fn(m) returns a Poisson sample of size m

fx=@(m)exp(mu+sig*randn(m,1)); % fx(m) returns a lognormal sample of size m


%%
% double for loop
tic;
y=zeros(n,1);
for i=1:numel(y),
  nr=fn(1);
  for j=1:nr,
    y(i)=y(i)+fx(1);
  end;
end;
toc;
mean(y)


%%
% loop over number of samples
tic;
nr=fn(n);
y=zeros(size(nr));
for i=1:numel(nr),
  y(i)=sum(fx(nr(i)));
end;
toc;
mean(y)


%%
% arrayfun
tic;
y=arrayfun(@(m)sum(fx(m)),fn(n));
toc;
mean(y)


%%
% vectorized - collapsing array
tic;
nr=fn(n);
mx=max(nr);
y=sum(reshape(fx(n*mx),n,mx).*bsxfun(@le,1:mx,nr),2);
toc;
mean(y)


%%
% vectorized - sample and discard
tic;
nr=fn(n);
y=zeros(size(nr));
for i=1:max(nr),
  y=y+fx(numel(y)).*(i<=nr);
end;
toc;
mean(y)


%%
% block-wise vectorized samples (discarding)
tic;
nr=fn(n);
y=zeros(size(nr));
nb=1000;
for i=1:ceil(numel(nr)/nb),
  imin=nb*(i-1)+1;
  imax=min(nb*i,numel(nr));
  x1=nr(imin:imax);
  y1=zeros(size(x1));
  for j=1:max(x1),
    y1=y1+fx(nb).*(j<=x1);
  end;y(imin:imax)=y1;
end;
toc;
mean(y)


%%
% vectorized - shrinking sample
tic;
nr=fn(n);
y=zeros(size(nr));
for i=1:max(nr),
  b=(i<=nr);
  y(b)=y(b)+fx(sum(b));
end;
toc;
mean(y)

%%
% vectorized - hybrid sample/discard and shrinking sample
tic;
nr=fn(n);
y=zeros(size(nr));
mx=floor(mean(nr));
for i=1:mx,
  y=y+fx(n).*(i<=nr);
end;
for i=(mx+1):max(nr),
  b=(i<=nr);
  y(b)=y(b)+fx(sum(b));
end;
toc;
mean(y)


%% Now using a more efficient Poisson sampler
fn=@(m)randraw('po',lam,m,1); % FEX #7309 (10x faster than poissrnd)


%%
% vectorized - shrinking sample
tic;
nr=fn(n);
y=zeros(size(nr));
for i=1:max(nr),
  b=(i<=nr);
  y(b)=y(b)+fx(sum(b));
end;
toc;
mean(y)


%%
% vectorized - hybrid sample/discard and shrinking sample
tic;
nr=fn(n);
y=zeros(size(nr));
mx=floor(mean(nr));
for i=1:mx,
  y=y+fx(n).*(i<=nr);
end;
for i=(mx+1):max(nr),
  b=(i<=nr);
  y(b)=y(b)+fx(sum(b));
end;
toc;
mean(y)


%% Theoretical speed limit
% time used by randn alone (assumes randn is optimal and ignores Poisson generator)
tic;
for i=1:ceil(lam)
  x=exp(randn(n,1));
end;
toc;


%%