probability or statistics - Generate a random PDF

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?

Answer

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},
    MixtureDistribution[
      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]