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.
Abs
and use the form in a normalized distribution with support from 0 to Infinity, sample viaInverseCDF
, then multiply samples by a random choice from{-1,1}
, since the distribution is symmetric around 0. I ran a few thousand runs of sampling 100 at a time per your post this way, no hangs. $\endgroup$