RandomVariate drawn from a custom distribution function takes forever to finish

Question

I'm using the following custom-defined PDF:

dist = ProbabilityDistribution[(
   Cos[1/(4 Abs[x])] - 
    2 Cos[1/(4 Abs[x])] FresnelC[1/(Sqrt[2 \[Pi]] Sqrt[Abs[x]])] + 
    Sin[1/(4 Abs[x])] - 
    2 FresnelS[1/(Sqrt[2 \[Pi]] Sqrt[Abs[x]])] Sin[1/(4 Abs[x])])/(
   2 Sqrt[2 \[Pi]] Abs[x]^(3/2)), {x, -Infinity, Infinity}];

When I do something like (I'm playing around with statistics of a very special cooked distribution and sum of random deviates)

Total[RandomVariate[dist,100]]

I never get a result. I only found out why after I did

Monitor[Sum[RandomVariate[dist], {i, 1, 100}], i]

Sometimes it chokes at 5, sometimes it chokes at 10, or 19, it really is...eh, random. I also noticed that sometimes the routine RandomVariete takes longer, sometimes shorter (and sometimes it just chokes). I guess there is some internal problem Mathematica is having with generating random numbers from custom user-defined probability distributions? How can this task be accomplished? Does the distribution must have analytic expression for the inverse of the distribution function? I can use FindRoot together with the expression for CDF and use the uniform distribution on (0,1) to get a random variate from the PDF, but how come Mathematica struggles with this?

Answer 1

Yes, a custom ProbabilityDistribution can be tricky for RandomVariate. The easiest way to get a result is by using the undocumented MCMC (Markov Chain Monte Carlo) sampler like so:

samples = RandomVariate[dist, 10000, Method -> {"MCMC",
     "Thinning" -> 1, "InitialVariance" -> 1,
     "InitialPoint" -> 0.01, "BurnInPeriod" -> 100
    }
];
Total[samples]

-283976.

I included the sub-options you can supply. To get good results, you may need to tweak these options since MCMC is not a one-size-fits-all method for sampling distributions. The "Thinning" an "BurnInPeriod" options are probably the most important to play around with. A "BurnInPeriod" is necessary to make the MCMC chain reach equilibrium while the value of the "Thinning" option determines how many samples get skipped in the chain before the next one is drawn. Increasing this option is useful to reduce correlations between subsequent samples.

Alternatively, there is the possibility that you could write the distribution as a transformation of another one. In that case, it is recommended to use TransformedDistibution for sampling. For example, if you want to sample X^3 where X is distributed as NormalDistribution[], I recommend you use

RandomVariate[TransformedDistribution[x^3, x \[Distributed] NormalDistribution[]]]

instead of computing the PDF of X^3 and stuffing that into ProbabilityDistribution.

Extra information about MCMC sampling

If you really want to go into the nuts and bolts of the MCMC sampler that Mathematica has, you can use Statistics`MCMC`BuildMarkovChain and Statistics`MCMC`MarkovChainIterate to exert a little more control over the sampler.

First of all, this is how you can find the available methods:

In[20]:= Statistics`MCMC`MCMCData[]

Out[20]= {"Metropolis", {"Metropolis", "Log"}, "IndependentMetropolis", {"IndependentMetropolis", "Log"}, 
"TransformedMetropolis", {"TransformedMetropolis", "Log"}, "AdaptiveMetropolis", {"AdaptiveMetropolis", "Log"}, "Hamiltonian", "Gibbs"}

Once you've selected the method you want to use, you can find out more about how to use it with these two lines. I will use the adaptive Metropolis method in log-space since it works quite well for a large class of distributions. In general you'll want to use log-densities whenever possible when working with probability densities since you will run into fewer problems with machine number underflows that way.

Statistics`MCMC`MCMCData[{"AdaptiveMetropolis", "Log"}, "Usage"]
Statistics`MCMC`MCMCData[{"AdaptiveMetropolis", "Log"}, "Example"]

Here is an example of how to build a Markov chain object and sample from it:

obj = Statistics`MCMC`BuildMarkovChain[{"AdaptiveMetropolis", "Log"}][
   RandomVariate[NormalDistribution[]], (* Random starting point *)
   Function[x, Evaluate[FullSimplify[Log[PDF[dist, x]], x \[Element] Reals]]], (*Log-density to sample from *)
   {1, 10} (* for the first 10 steps used 1 SD for proposal, then switch to adaptive *)
];
Statistics`MCMC`MarkovChainIterate[obj, 1000]; (* Burn in for 1000 steps *)
samplesMCMC = Statistics`MCMC`MarkovChainIterate[obj, {1000, 10}]; (* Sample 1000 points; only keep every 10th point *)
obj["AcceptanceRate"] (* Check the acceptance rate *)

As a bonus, MCMC samplers can work with un-normalized probability densities so you don't have to go through the trouble of normalizing your PDF.

Finally, Statistics`MCMC`BuildMarkovChain also has a WorkingPrecision option that you can use to do the sampling in arbitrary precision arithmetic instead of machine precision. This can help if your PDF is tricky to evaluate numerically because of potential rounding errors.

Answer 2

Per my comment, a possible workaround to the hanging problem:

d = ProbabilityDistribution[(
  Cos[1/(4 x)] - 2 Cos[1/(4 x)] FresnelC[1/(Sqrt[2 \[Pi]] Sqrt[x])] + 
   Sin[1/(4 x)] - 
   2 FresnelS[1/(Sqrt[2 \[Pi]] Sqrt[x])] Sin[1/(4 x)])/(
  2 Sqrt[2 \[Pi]] x^(3/2)), {x, 0, Infinity}, Method -> "Normalize"];

numsamp = 20;

samps = InverseCDF[d, RandomReal[1, numsamp]]*RandomChoice[{-1, 1}, numsamp]

{0.812671,-4.5795,103751.,2.24011,-0.252999,77.4749, -11.0699,0.340867,0.0431401,7.19123,6967.5,-320.257, -0.538424,1.40892,0.0402734,1.24583,-1.20498,-0.949599,0.342374,-1.22637}