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

Contents

• Tests using POISSRND (the stats toolbox Poisson sampler)
• Now using a more efficient Poisson sampler
• Theoretical speed limit

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)

Elapsed time is 237.690146 seconds.

ans =

   16.4879

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)

Elapsed time is 15.586183 seconds.

ans =

   16.4892

arrayfun

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

Elapsed time is 15.676805 seconds.

ans =

   16.4779

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)

Elapsed time is 3.794724 seconds.

ans =

   16.4860

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)

Elapsed time is 3.966500 seconds.

ans =

   16.4790

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)

Elapsed time is 3.753551 seconds.

ans =

   16.4920

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)

Elapsed time is 2.959969 seconds.

ans =

   16.4847

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)

Elapsed time is 2.813371 seconds.

ans =

   16.4752

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)

Elapsed time is 1.314490 seconds.

ans =

   16.5029

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)

Elapsed time is 1.156294 seconds.

ans =

   16.4855

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;

Elapsed time is 0.721762 seconds.

Published with MATLAB® 7.9