A probability distribution can be created in Mathematica (I am using 8.0.1) with e.g.

distribution1 = ProbabilityDistribution[(Sqrt[2])/Pi*(1/((1+x^4))), 

Random variates from this distribution can be created with RandomVariate easily:

Histogram[dataDistribution1](*Just an optical control*)

How can I create random variates from a 2-dimensional (multivariate) probability distribution? Let's say my 2-dimensional distribution is of the following form (very similar to the previous one):

distribution2=ProbabilityDistribution[((((Sqrt[2] π^(3/2) Gamma[5/4])/

I thought it would be logical to try


But that does not work. I get the following message then:

RandomVariate::noimp: Sampling from ProbabilityDistribution[Gamma[3/4]…}] is not implemented.

I tried a lot of variations of this approach to create random variates from such a distribution but without any success. Either I am doing a lot of things wrong (very likely in my experience ;-) or Mathematica cannot deliver random variates from 2-dimensional probability distributions. But in the help of RandomVariate (under “Simulate a multivariate continuous distribution”) one can see that this should be possible.

Perhaps somebody of you can tell me, how I can generate random variates from 2 dimensional probability distributions? I would be very happy about any help!

  • $\begingroup$ Huh, from the docs I would expect the same thing. I'm looking forward to answers as well. $\endgroup$
    – tkott
    Mar 6, 2012 at 17:40
  • $\begingroup$ Looks like you might need to go with a hand-made random number generator: reference.wolfram.com/mathematica/tutorial/… $\endgroup$
    – tkott
    Mar 6, 2012 at 17:58
  • $\begingroup$ Hi! Thanks for commenting. Great to see that you share my expectation. I hope there will be more comments/answers too. $\endgroup$
    – partial81
    Mar 6, 2012 at 17:59
  • 1
    $\begingroup$ The obvious way is to integrate out one of the variables to obtain its marginal distribution. Sample from the marginal, then sample the other variable from the conditional distribution. This process is quite slow, though, and requires the marginal to have a closed form that Mathematica can find. You're probably better off using Gibbs sampling. $\endgroup$
    – whuber
    Mar 6, 2012 at 19:29
  • 2
    $\begingroup$ Is sampling from a multi-dimensional distribution now possible in Mathematica 10? $\endgroup$
    – Guillochon
    Aug 28, 2014 at 15:37

3 Answers 3


Mathematica v8 does not provide support for automated random number generation from multivariate distributions, specified in terms of its probability density functions, as you have already discovered it.

At the Wolfram Technology conference 2011, I gave a presentation "Create Your Own Distribution", where the issue of sampling from custom distribution is extensively discussed with many examples.

You can draw samples from the particular distribution at hand by several methods. Let

di = ProbabilityDistribution[
   Beta[3/4, 1/2]/(Sqrt[2] Pi^2) 1/(1 + x^4 + y^4), {x, -Infinity, 
    Infinity}, {y, -Infinity, Infinity}];

Conditional method

The idea here is to first generate the first component of the vector from a marginal, then a second one from a conditional distribution:

md1 = ProbabilityDistribution[
   PDF[MarginalDistribution[di, 1], x], {x, -Infinity, Infinity}];

cd2[a_] = 
   Simplify[PDF[di, {a, y}]/PDF[MarginalDistribution[di, 1], a], 
    a \[Element] Reals], {y, -Infinity, Infinity}, 
   Assumptions -> a \[Element] Reals];

Then the conditional method is easy to code:

diRNG[len_, prec_] := Module[{x1, x2},
  x1 = RandomVariate[md1, len, WorkingPrecision -> prec];
  x2 = Function[a, RandomVariate[cd2[a], WorkingPrecision -> prec]] /@
  Transpose[{x1, x2}]

You can not call it speedy:

In[196]:= AbsoluteTiming[sample1 = diRNG[10^3, MachinePrecision];]

Out[196]= {20.450045, Null}

But it works:

enter image description here

Transformed distribution method

This is somewhat a craft, but if such an approach pans out, it typically yields the best performing random number generation method. We start with a mathematical identity $$ \frac{1}{1+x^4+y^4} = \int_0^\infty \mathrm{e}^{-t(1+x^4+y^4)} \mathrm{d} t = \mathbb{E}_Z( \exp(-Z (x^4+y^4))) $$ where $Z \sim \mathcal{E}(1)$, i.e. $Z$ is exponential random variable with unit mean. Thus, for a random vector $(X,Y)$ with the distribution in question we have $$ \mathbb{E}_{X,Y}(f(X,Y)) = \mathbb{E}_{X,Y,Z}\left( f(X,Y) \exp(-Z X^4) \exp(-Z Y^4) \right) $$ This suggests to introduce $U = X Z^{1/4}$ and $V = Y Z^{1/4}$. It is easy to see, that the probability density function for $(Z, U, V)$ factors: $$ f_{Z,U,V}(t, u, v) = \frac{\operatorname{\mathrm{Beta}}(3/4,1/2)}{\sqrt{2} \pi^2} \cdot \frac{1}{\sqrt{t}} \mathrm{e}^{-t} \cdot \mathrm{e}^{-u^4} \cdot \mathrm{e}^{-v^4} $$ It is easy to generate $(W, U, V)$, since they are independent. Then $(X,Y) = (U, V) W^{-1/4}$, where $f_W(t) = \frac{1}{\sqrt{\pi}} \frac{1}{\sqrt{t}} \mathrm{e}^{-t}$, i.e. $W$ is $\Gamma(1/2)$ random variable.

This gives much more efficient algorithm:

  prec_] := (RandomVariate[NormalDistribution[], len, 
       WorkingPrecision -> prec]^2/2)^(-1/4) RandomVariate[
    1/(2 Gamma[5/4]) Exp[-x^4], {x, -Infinity, Infinity}], {len, 2}, 
   WorkingPrecision -> prec]

Noticing that $|W|$ is in fact a power of gamma random variable we can take it much further:

In[40]:= diRNG3[len_, prec_] := 
 Power[RandomVariate[GammaDistribution[1/4, 1], {len, 2}, 
     WorkingPrecision -> 
      prec]/(RandomVariate[NormalDistribution[], len, 
        WorkingPrecision -> prec]^2/2), 1/4] RandomChoice[
   Range[-1, 1, 2], {len, 2}]

In[42]:= AbsoluteTiming[sample3 = diRNG3[10^6, MachinePrecision];]

Out[42]= {0.7230723, Null}

Rejection method

Here the idea is to sample from a relatively simple to draw from hat distribution. It is again a craft to choose a good one. Once the one is chosen, we exercise the rejection sampling algorithm:

enter image description here

In the case at hand, a good hat is the bivariate T-distribution with 2 degrees of freedom, as it is easy to draw from, and it allows for easy computation of the scaling constant:

In[49]:= Maximize[(1/(1 + x^4 + y^4))/
  PDF[MultivariateTDistribution[{{1, 0}, {0, 1}}, 2], {x, y}], {x, y}]

Out[49]= {3 Pi, {x -> -(1/Sqrt[2]), y -> -(1/Sqrt[2])}}

This gives another algorithm:

diRNG4[len_, prec_] := Module[{dim = 0, bvs, u, res},
  res = Reap[While[dim < len,
      bvs = 
       RandomVariate[MultivariateTDistribution[{{1, 0}, {0, 1}}, 2], 
        len - dim, WorkingPrecision -> prec];
      u = RandomReal[3/2, len - dim, WorkingPrecision -> prec];
      bvs = 
        Sign[(Total[bvs^2, {2}]/2 + 1)^2 - u (1 + Total[bvs^4, {2}])],
      dim += Length[Sow[bvs]];
      ]][[2, 1]];
  Apply[Join, res]

This one proves to be quite efficient as well:

In[77]:= AbsoluteTiming[sample4 = diRNG4[10^6, MachinePrecision];]

Out[77]= {0.6910000, Null}
  • 1
    $\begingroup$ +1 Beautifully done! This will serve as a working reference for many people. $\endgroup$
    – whuber
    Mar 6, 2012 at 20:24
  • $\begingroup$ Hi @Sasha! Sorry for my late reply. It took me a while to adapt your solutions to my real problem. But it works! Thanks for this detailed answer! As whuber said, it will surely help a lot of people with similar problems! $\endgroup$
    – partial81
    Mar 12, 2012 at 23:21
  • $\begingroup$ it is so helpful, thanks so much. $\endgroup$
    – user6196
    Mar 4, 2013 at 22:15
  • $\begingroup$ Thanks for the detailed answer. I was just wondering if it is possible to extend the Conditional Method to construct a 3 dimensional distribution? $\endgroup$
    – user29165
    Jul 21, 2014 at 13:17
  • $\begingroup$ @user29165 Most certainly. You need to start with finding the marginal distribution of a suitably chosen component, say x1. Then, given the value of x1, you need a conditional distribution of x2 given x1. This is found, for continuous distribution, as $f_{X_1,X_2}(x_1,x_2)/f_{X_1}(x_1)$, where $f_{X_1}(x_1)$ is the 1D marginal pdf, and $f_{X_1, X_2}(x_1, x_2)$ is the 2D marginal pdf. To generate the last component you would need to find and sample $X_3 \mid X_1, X_2$ whose pdf is $f_{X_1, X_2,X_3}(x_1,x_2,x_3)/f_{X_1,X_2}(x_1,x_2)$. Try this out on a simple example of uniform in a 3D box. $\endgroup$
    – Sasha
    Jul 21, 2014 at 17:26

There is some hidden MCMC (Markov Chain Monte Carlo) functionality in Mathematica that can come in handy, though there are obviously no guarantees about how well it works since it's unexposed.

Here's a simple example:

dist = ProbabilityDistribution[
   Beta[3/4, 1/2]/(Sqrt[2] Pi^2) 1/(1 + x^4 + y^4), {x, -Infinity, 
    Infinity}, {y, -Infinity, Infinity}];
samples = RandomVariate[dist, 10^3, Method -> {"MCMC", "Thinning" -> 1, "InitialVariance" -> 1}]

I do not know exactly if there are other methods/options you can feed into this MCMC sampler.

General MCMC

It's also worth checking out the context:


This lists the available methods:

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

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

Usage examples and details of a given method (say, {"AdaptiveMetropolis", "Log"}) can be requested:

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

First, you will need to build an MCMC object that holds the state of the random walk. For this example, you need to define a single-argument pure function that returns the log density of the probability distribution:

chain = Statistics`MCMC`BuildMarkovChain[{"AdaptiveMetropolis", "Log"}][
  {0, 0}, (* starting point *)
  Function[{\[FormalX]}, (* function that takes 2D points as a list *)
     Evaluate @ Log @ PDF[dist, {Indexed[\[FormalX], 1], Indexed[\[FormalX], 2]}]],
    0.5 * DiagonalMatrix @ {1, 1}, (* starting covariance matrix for generating proposals *)
    10 (* the step after which the chain starts generating adaptive proposals *)

You can now draw samples from the chain:

Statistics`MCMC`MarkovChainIterate[chain, 1000]; (* Burn in for 1000 steps *)
samples = Statistics`MCMC`MarkovChainIterate[chain, 
  {1000, 10} (* generate 1000 samples, taking a sample every 10th step in the chain *)

Since this is an MCMC method, you need to make sure you burn in the chain long enough before you start generating samples and be generally aware of the limitations of MCMC and how to diagnose problems with correlations between samples. It's also good to keep in mind that it's generally better to work with log-densities if you can, since that makes numerical underflows significantly less likely to happen.


In Mathematica 11.0 (perhaps earlier), RandomVariate[] works with multivariate distributions. E.g.

   0}, {{0.1^2, 0}, {0, (0.3*10^-3)^2}}], 10]

0.000138816}, {3.11169, -0.000342556}, {3.12357, -0.0000204628},
{3.18824, -0.000355518}, {3.21631, -0.000137532}, {3.12506, 
3.81474*10^-6}, {3.16097, -0.0000302336}, {3.20877, 
0.0000405303}, {3.29134, 0.000316891}, {3.15165, 0.000308802}}
  • 2
    $\begingroup$ But it does not work with a multivariate ProbabilityDistribution. This does not solve the OP's problem. $\endgroup$
    – Szabolcs
    Jan 15, 2020 at 14:13

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