# How to generate a RandomVariate of a custom distribution?

I'm trying to generate a pseudorandom variate out of a custom distribution. Suppose I want define a custom distribution, and for the sake of simplicity I define a Poisson distribution (the distribution I want is a derived mixture of two Poisson distributions):

custom[a_] := ProbabilityDistribution[(a^k Exp[-a])/k!, {k, 0, ∞, 1}];


When I now try to generate a RandomVariate[] of that distribution,

RandomVariate[custom[2], 10]


I obtain the following error:

RandomVariate::noimp: Sampling from ProbabilityDistribution[2^x/(E^2 x!),{x,0,∞,1}] is not implemented.

For the record, when the distribution is defined as a continuous distribution:

custom[a_] := ProbabilityDistribution[(a^k Exp[-a])/k!, {k, 0, ∞}];


there is no problem with the use of RandomVariate.

I'd appreciate some insight into this problem, thanks!

-
It would work if for instance your new distribution is derived from a known one via TransformedDistribution as in RandomVariate[ TransformedDistribution[x - 1, x \[Distributed] NormalDistribution[]] // Release, 15] – chris Nov 1 '12 at 16:28
"the distribution I want is a derived mixture of two Poisson distributions" - could you maybe specify how the two Poisson distributions are mixed? I presume you've tried MixtureDistribution[] already? – J. M. Nov 1 '12 at 16:32
FWIW: it looks to me (after some experimentation) that RandomVariate[] is not equipped to handle arbitrary discrete probability distributions; in general, you might have to roll your own algorithm. – J. M. Nov 1 '12 at 16:41
@J.M. why does RandomVariate[\[ScriptCapitalD], 15] from reference.wolfram.com/mathematica/ref/… works though? – chris Nov 1 '12 at 19:35
@J.M. and more importantly why is it wrong?? Show[Plot[PDF[\[ScriptCapitalD], x], {x, -5, 5}], RandomVariate[\[ScriptCapitalD], 1500]//Histogram[#, Automatic, "Probability"]&] – chris Nov 1 '12 at 19:40

The short answer is: I don't think it is built to always work for general distributions (though see Sjoerd C. de Vries's answer and continuous example below).
You might also want to have a look at the nice tutorial how to Create Your Own Distribution with Oleksandr Pavlyk.

Lets see what are your options

Mixture

For mixture distributions, following J. M., we can define a mixture as:

DD[μ_, ν_, a_, b_] =
MixtureDistribution[{a, b}, {PoissonDistribution[μ], PoissonDistribution[ν]}]


so that

RandomVariate[DD[4, 3, 1, 2], 16]

(* {5, 5, 6, 3, 5, 7, 9, 1, 3, 7, 5, 5, 3, 5, 2, 3} *)


Transformation

One can also use TransformedDistribution to define a new distribution from a known one, and draw from it:

RandomVariate[TransformedDistribution[x - 4,
x \[Distributed] PoissonDistribution[10]] // Release, 15]


Draw via Cumulative Distrubution

For a set of distributions for which the cumulative distribution exists, within which your example falls (though it probably differs from what you really want to do), the inverse of the cumulative distribution can be used to define a draw:

 CDF[custom[a], x]

(* Gamma[1 + Floor[x], a]/Gamma[1 + Floor[x]] *)


So for instance,

  g[x_] = CDF[custom[10], x];
h = Table[{g[x], x}, {x, 0, 40}] // N // Interpolation[#, InterpolationOrder -> 0] &;


which allows us to define the inverse of the CDF:

  h /@ (RandomReal[{0, 1}, {12}])

(*  {12., 15., 6., 17., 12., 8., 6., 9., 8., 13., 13., 12.} *)

Show[h /@ (RandomReal[{0, 1}, {125000}]) // Histogram[#, Automatic, "Probability"] &,
Plot[PDF[custom[10], x - 1/2], {x, 0, 20}]]


Empirical Distribution

For a more general class of distribution, there is a twisted route within the current set of existing Mathematica functions: if you can draw samples which satisfy the distribution, then you can derive an EmpiricalDistribution (thanks J. M.) from your sample and then draw through it. It might be useful if it's costly to draw in the first place.

Continuous distribution where ProbabilityDistribution works with Random Variate

norm = Integrate[Exp[-x^2]/(c^2 + x^2), {x, -Infinity, Infinity},Assumptions -> c > 0];
DD[c_] = ProbabilityDistribution[1/norm Exp[-x^2]/(c^2 + x^2),
{x, -Infinity,Infinity},Assumptions -> c > 0];


Let's check that the random variate sample is consistent

 data =RandomVariate[DD[4], 1500];

DDE = EmpiricalDistribution[data];
{Plot[CDF[DD[4], x] // Release, {x, -4, 4},PlotStyle -> Directive[Red]],
Plot[CDF[DDE, x], {x, -4, 4}]} // Show


Note that it can involve more than one parameter:

 norm = Integrate[Exp[-β x^2]/(c^2 + x^2),{x,-Infinity,Infinity},Assumptions->{c > 0, β>0}];
DD[c_, β_] = ProbabilityDistribution[1/norm Exp[-β x^2]/(c^2 + x^2),
{x, -Infinity, Infinity}, Assumptions -> (c > 0 && β > 0)];

RandomVariate[DD[1, 1], 4]

(* {0.149204,-0.644442,-0.632239,-0.266628} *)


Final Note

For 2D generalization see RandomVariate from 2-dimensional probability distribution

-
Since OP is sampling from a discrete distribution, EmpiricalDistribution[] might be more appropriate... – J. M. Nov 1 '12 at 16:47

I have a hunch that the problem lies in the use of Infinity in combination with a finite step. So, it's not an issue of continuous vs discrete. Replace $\infty$ with a sufficiently large number (with the factorial involved 100 should be OK) and it works:

custom[a_] := ProbabilityDistribution[(a^k Exp[-a])/k!, {k, 0, 100, 1}];

RandomVariate[custom[2], 10]


{2, 2, 1, 4, 2, 2, 1, 4, 1, 1}

Let's test whether a set of 10 million random draws from the custom distribution are described by a Poisson distribution with the same parameter:

EstimatedDistribution[RandomVariate[custom[2], 10000000], PoissonDistribution[λ]]


PoissonDistribution[1.9999276]

EstimatedDistribution[RandomVariate[custom[8.5], 10000000], PoissonDistribution[λ]]


PoissonDistribution[8.5006323]

-
so it would be a kind of bug that it does not work? – chris Nov 2 '12 at 8:40
@chris One man's bug is another man's feauture ;-) Anyway, there are various Mathematica functions like Sum that handle infinite series and can deal with them in a symbolic way. It just seems that ProbabilityDistribution isn't one of them. – Sjoerd C. de Vries Nov 2 '12 at 8:57
I see your point; nonetheless its a bit unfortunate because quite a few distributions have infinities in their domain. – chris Nov 2 '12 at 11:00
ProbabilityDistribution is not able to sample custom univariate discrete distribution because it does not know where the majority of the probability is concentrated. Sjoerd's workaround is enabling, provided truncated density is properly re-normalized, because he manually truncated the density to where most density is concentrated. – Sasha Nov 2 '12 at 16:11
This is actually very helpful, thanks! – Rafal Dec 6 '12 at 17:13

This is trivially easy with mathStatica for Mathematica, for ANY arbitrary discrete distribution. For the original poster's problem:

f = Exp[-a] a^x / x!;
domain[f] = {x, 0, Infinity} &&  {a > 0}  &&  {Discrete};


Here are 10 pseudo-random drawings:

RandomNumber[10, f /. a -> 2]
{2, 2, 1, 4, 3, 0, 1, 1, 1, 7}


And here are 100,000 more:

data = RandomNumber[100000, f /. a -> 2];


Next, compare a frequency plot of the empirical sample data we have just generated (the red dots) to the parent symbolic pmf (in blue):

FrequencyPlotDiscrete[data, f /. a -> 2]


Splendido. Just as easy for anything else :)

We have a chapter on these sorts of issues... see Chapter 3 of: Rose and Smith (Nov. 2011), Mathematical Statistics with Mathematica.

-