8
$\begingroup$

Consider the following code where we first take 1000 independent samples from a Poisson distribution, and then take 1000 independent samples from 1000 different Poisson distributions:

list = RandomVariate[NormalDistribution[0.5, 0.1],1000];
RandomVariate[PoissonDistribution[0.5], 1000] // Timing // First
RandomVariate[PoissonDistribution[#]] & /@ list // Timing // First

which on my machine outputs something like 0.000343 and 0.137541, respectively. Note that it obviously doesn't matter if list contains different numbers or identical numbers.

Clearly RandomVariate has some kind of overhead. Is there a way to improve the efficiency of sampling from a list of distributions that vary only in their parameters?

$\endgroup$
6
  • 1
    $\begingroup$ RandomInteger[PoissonDistribution[#]] & /@ list is slightly faster $\endgroup$ Commented Oct 6, 2015 at 17:27
  • $\begingroup$ Poisson-specific possible duplicate: (56180) and general-distribution possible duplicate (75303). Note also that attempting to compile with a Poisson distribution won't work, as per (1124) $\endgroup$ Commented Oct 6, 2015 at 17:41
  • $\begingroup$ And in fact, another Poisson-based possible duplicate (35433) $\endgroup$ Commented Oct 6, 2015 at 17:44
  • 1
    $\begingroup$ Thanks. The last one you mentioned seems the most like a duplicate, and there are some good answers there, but I am most interested in a distribution-independent solution, which does not seem to be asked or answered elsewhere. $\endgroup$
    – Ian Hincks
    Commented Oct 6, 2015 at 18:03
  • $\begingroup$ Well (75303) is probably closest - for a distribution-independent solution, you can turn either to a C++ solution, or the clever method presented by @ciao $\endgroup$ Commented Oct 6, 2015 at 18:06

1 Answer 1

5
$\begingroup$

Per your example case:

(* setup (you'll pay time here, but once done, it's done) *)
dist = ParameterMixtureDistribution[PoissonDistribution[m],
                                     m \[Distributed] NormalDistribution[1/2, 1/10]];

pdf = PDF[dist, z];

(* pick some end-point with negligible probability *)
(* Here, using 30 truncates tail with total p of ~6.34*10^-40 *)
dc = Rule @@ Transpose@Table[{pdf /. z -> zz, zz}, {zz, 0, 30}];

(* use it and compare *)
result = RandomChoice[dc, 1000] // Timing // First

list = RandomVariate[NormalDistribution[0.5, 0.1], 1000];
resultold=RandomVariate[PoissonDistribution[#]] & /@ list // Timing // First

0.

0.842405

Timings on netbook...

$\endgroup$
3
  • $\begingroup$ I miss the cigar lounge $\endgroup$ Commented Oct 7, 2015 at 3:25
  • 4
    $\begingroup$ @belisariusisforth: there as I type - will take a puff in your honor ;-} $\endgroup$
    – ciao
    Commented Oct 7, 2015 at 5:33
  • $\begingroup$ PDF[dist, z] gives Undefined in 10.3.1. $\endgroup$
    – Kuba
    Commented Jan 31, 2016 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.