How can I generate a random probability density function $p(x)$ in the domain $[0,1]$? Is there a single Mathematica function to do so?


2 Answers 2


There are many ways to generate random distributions on the unit interval. I do not know of a function in Mathematica that will do that by itself, but it is straightforward to write a simple function that will.

The approach presented here is based on random Bernstein polynomials, which in this context become random mixtures of beta distributions. Here is a reference:

Sonia Petrone (1999) Bayesian Density Estimation Using Bernstein Polynomials, The Canadian Journal of Statistics, Vol. 27, No. 1 (Mar., 1999), pp. 105-126

Both the number of mixture components and the mixture weights are random. There is a certain amount of flexibility in choosing the underlying random distributions. In what follows, I use the geometric distribution to determine the number of mixture components and the Dirichlet distribution to determine the mixture weights.

As I have implemented this, there are two parameters to be chosen by the user, one for each of the underlying distributions. The first parameter is the mean geometric distribution (which must be greater than or equal to one).

The second parameter is the concentration parameter for the Dirichlet distribution, which determines how concentrated the distribution is around its mean. (The concentration parameter must be positive.) I have fixed the mean at $(1/k, \ldots, 1/k)$ where $k$ is the number of mixture components. As a consequence, the average random distribution is uniform. (There are a number of ways to produce non-uniform average random distributions. An important way involves applying a suitable transformation to the beta distributions.)

Programming notes:

I do not use the built-in DirichletDistribution because it only generates $k-1$ random weights, instead of $k$. In principle the missing weight could be computed via $w_k = 1 - \sum_{j=1}^{k-1} w_j$, but in practice this fails because $w_k$ ends up being negative too often to be useful. (It's a very small negative number, but negative nonetheless.) Instead I generate the gamma variates that Dirichlet variates are based on. (There is no need to normalize these gamma variates, since MixutureDistribution does that automatically.)

For efficiency, I use Chop to kill off those mixture components with trivially small mixture weights. (This is not necessary.)

randomDistribution[xi_, alpha_] :=
  With[{k = RandomVariate[GeometricDistribution[1/xi]] + 1},
      RandomVariate[GammaDistribution[alpha/k, 1], k], 
      Table[BetaDistribution[j, k - j + 1], {j, k}]
    ] // Chop

Here is some usage. I have set the mean of the geometric distribution to 100 and the concentration parameter to 1. (The use of Evaluate in the plot can speed things up significantly.)

dist = randomDistribution[100, 1];
Plot[PDF[dist, x] // Evaluate, {x, 0, 1}, PlotRange -> {0, All}, Filling -> 0]

For a given function $p(x)$ to satisfy your requirement it needs to meet the constraints:

$$\int_0^1 p(x)dx = 1; \quad p(0 \le x \le 1) \ge 0; \quad p(x < 0 \cup x > 1)=0$$

How to make it random is an open question. One approach is a function that passes through several random points in the unit square 0<=x<=1; 0<=y<=1, as well as necessarily through some point at x=0 and x=1 for good measure. That's easily done by

list = {{0, RandomReal[]}}~Join~RandomReal[{0, 1}, {5, 2}]~
   Join~{{1, RandomReal[]}}; (*get random points *)
func = Interpolation[list, InterpolationOrder -> 1]; (*construct interpolating function *)
norm = Integrate[func[x], {x, 0, 1}]; (* find the coefficient to normalize it *)
Plot[func[x]/norm, {x, 0, 1}] (* plot the result *)

enter image description here

The randomness of this function is... questionable. If we take a bunch of random points with their x and y coordinates independently distributed according to this, we'll get the following plot:

  ProbabilityDistribution[func[x]/norm, {x, 0, 1}], {100000, 2}];

enter image description here

Of course, we may want more "randomness, as the behavior of the random function in between the random points it goes through is very much deterministic.

So let's run the above code with RandomReal[{0,1},{100,2}]... except it won't work. With a function like that RandomVariate takes 3 seconds to generate just two(!) points.

So I'll simply compute the CDF to make a point:

list = {{0, RandomReal[]}}~Join~RandomReal[{0, 1}, {100000, 2}]~
   Join~{{1, RandomReal[]}};
list = SortBy[list, First];
{{0, 0}}~Join~({#3 - #1, (#3 - #1) (#4 + #2)/2} & @@@ 
     Flatten /@ Partition[list, 2, 1]) // Accumulate;
cdf = Interpolation[%, InterpolationOrder -> 1];
Plot[{cdf[x]/cdf[1], x - 0.01}, {x, 0, 1}]

enter image description here

As you can see, in the limit of "very many points" the CDF is hardly different from the CDF of a uniform distribution (shown here with a slight offset to distinguish between the two functions). To progress further, one will need a more rigorous definition of a "random function".


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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