Speed up RandomVariate on a custom probability distribution

Question

I have created a custom probability distribution:

aDistLogisticMakeham[α_, β_, s_, cMakeham_] := ProbabilityDistribution[(α Exp[β t]/(1 + (α s/\
        β) (Exp[β t] - 1)) + 
        cMakeham) ((β - α s + α s Exp[β t])/\
        β)^(-1/s) Exp[-cMakeham t], {t, 0, ∞}]

Which I then use to create a function which draws a single independent sample:

fRandomLogisticMakeham[α_, β_, s_, cMakeham_] :=  RandomVariate[
   aDistLogisticMakeham[α, β, s, cMakeham], {1}][[1]]

I then want to apply this sampler to a list of inputs, say of length 100:

{α, β, s, cMakeham} = Table[RandomVariate[ExponentialDistribution[1], {100}], {i, 1, 4, 1}];
inputs = Thread[{α, β, s, cMakeham}];

So, I do:

fRandomLogisticMakeham @@@ inputs // Timing

The timing for this is 14s on my machine. The problem is that I need to use this sampler on a list of length 100,000 or more.

I know I could use rejection sampling, or importance sampling, but wanting to avoid these, I was wondering whether there were any changes to my existing code/framework that would make it run faster?

I have tried to find the inverse CDF, using:

InverseCDF[aDistLogisticMakeham[α, β, s, cMakeham],t]

In order to do inverse transform sampling, but Mathematica gets stuck here (I tried adding positivity assumptions about the parameters to my distribution definition, but this didn't help either.)

Edit: I have now changed the density above, as per Wolfies' remark below - the density is now a proper probability distribution. However, the problem remains - the above is still very very slow.

Any help would be much appreciated!

Best,

Ben

Answer 1

A few thoughts:

Notation

Defining:

{α, β, s, cMakeham} = Table[ blah matrix]  

... looks like an invitation for trouble, since you are setting up $\alpha$, $\beta$ as vectors etc ... when those are the exact same symbols used in your pdf, for parameters. You may have avoided the issue by setting up your pdf as a function of the parameters, but it is still an invitation for problems. Moreover, for stack exchange, why use $\alpha$ and $beta$ when $a$ and $b$ are easier. To avoid all these issues, let us write your pdf as:

pdf = (c + (a*E^(b*t))/(1 + (a*(-1 + E^(b*t))*s)/b))/(E^(c*t)*((b-a*s+ a*E^(b*t)*s)/b)^s^(-1))

Why your problem takes time

You then consider say 100 realisations of the parameter vector $\{a,b,s,c\}$:

inputs = Table[RandomVariate[ExponentialDistribution[1], {100}], 4] // Transpose;

Importantly, this means that the Timing problems you are incurring are part of the design of your problem. Usually, when one wishes to generate pseudorandom drawings from a distribution, we fix the parameter values, lock everything down to one specific form of the distribution ... and then generate say 10,000 drawings from that distribution. That means we can spend time setting up a numerical inverse CDF, or find a similar shaped distribution that we know how to generate from and use the Rejection Method etc

But we can't do that for your example, because ON EACH AND EVERY RANDOM DRAWING, you want to change the parameter values of the parent pdf. Every time you draw a value, you do so from a different shaped pdf ... which means, every time, we have to start the whole problem from fresh. You don't want 100,000 drawings from 1 pdf ... you want 1 drawing from 100,000 different pdf's. As such, what is surprising, in my view, is not that Mma is slow (as you suggest), but quite the opposite, how fast the inbuilt generator is, given what you are asking it to do... namely, to solve the whole problem fresh for each iteration.

Speeding things up: Method 1 - Parallel

Given that you are running the problem from fresh on every drawing, one easy and obvious way to speed things up is to simply run the entire operation in parallel:

 ParallelMap[ RandomVariate[ProbabilityDistribution[
     pdf /. {a -> #[[1]], b -> #[[2]], s -> #[[3]], c -> #[[4]]}, 
        {t, 0, Infinity}]] & , 
     inputs] // AbsoluteTiming

This takes less than 1 second on my Mac, for the same 100 samples you are considering (which took 14 seconds on your machine). The same code does 10,000 samples in abut 85 seconds on my Mac ... but that's really not so bad, given that you are actually sampling from 10,000 different distributions.

REALLY Speeding things up: Method 2 - Symbolic CDF and numeric inversion

It is possible to find the exact symbolic cdf, say $F(t)$:

 F = 1 - 1/(E^(c*t) ((b + a (-1 + E^(b*t)) s)/b)^s^(-1))

Although we cannot invert it symbolically, we can invert it using numerical methods which ... subject to some checking ... should generally work fine. The idea is we want to solve:

$$t=F^{-1}(u)$$

where $u$ = RandomReal. To illustrate that this works, let us consider the parameter case:

Fparam = F /. {a -> .1, b -> 3, s -> .3, c -> .2};

Then, 50000 pseudorandom drawings will be given by:

  data = ParallelTable[  t /. FindRoot[ RandomReal[] == Fparam, {t, 1}],  {50000}];  

The following diagram compares:

• the theoretical pdf (red dashed curve) to the
• the empirical pdf of the generated data (squiggly blue)

when {a -> .1, b -> 3, s -> .3, c -> .2} ... and all seems fine.

All that remains to do is adjust this for the problem at hand, with the OP's varying parameter inputs on each drawing. Our pseudorandom data is then:

data1 = ParallelMap[ t /. FindRoot[ 
         RandomReal[] == F /. {a->#[[1]], b->#[[2]], s->#[[3]], c->#[[4]]}, {t, 1}]&, 
           inputs]; // AbsoluteTiming

{0.031362, Null}

This is more than 100 times faster than the OP's code.

It is possible that, for some parameter combinations, that funny numerical results will be obtained (e.g. with small complex components). We can just drop these funny values and make sure all our data is positive:

 cleandata = Select[data1, Positive];

... and all should be resolved.

Answer 2

The trouble is, as you note the CDF can not be computed analytically so RandomVariate needs to numerically crunch out the inverse cdf for each value.

This does so directly, faster by a factor of 6..

randvariate[inputs_] := Module[{target = RandomReal[]},
  res = c /. 
    First@NDSolve[ { 
       c'[t] == PDF[aDistLogisticMakeham[Sequence @@ inputs]][t]  , 
       c[0] == 0,
       WhenEvent[c[t] == target, "StopIntegration"]},
      c, {t, 0, Infinity}, MaxStepSize -> .01];
  res["Domain"][[1, 2]]
  ]
(r2 = randvariate /@ inputs) // Timing // First

Note I have not carefully validated this, but histograms of the results look the same.