Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
2  
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
1  
"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
1  
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
show 8 more comments

3 Answers

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}]]

histogram and PDF of custom distribution

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

Mathematica graphics

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

share|improve this answer
1  
Since OP is sampling from a discrete distribution, EmpiricalDistribution[] might be more appropriate... –  J. M. Nov 1 '12 at 16:47
add comment

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]

share|improve this answer
    
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
1  
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
add comment

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]

freqplot

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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