Abstract

The method of potential solutions of Fokker-Planck equations is used to develop a transport equation for the joint probability of N coupled stochastic variables with the Dirichlet distribution as its asymptotic solution. To ensure a bounded sample space, a coupled nonlinear diffusion process is required: the Wiener processes in the equivalent system of stochastic differential equations are multiplicative with coefficients dependent on all the stochastic variables. Individual samples of a discrete ensemble, obtained from the stochastic process, satisfy a unit-sum constraint at all times. The process may be used to represent realizations of a fluctuating ensemble of N variables subject to a conservation principle. Similar to the multivariate Wright-Fisher process, whose invariant is also Dirichlet, the univariate case yields a process whose invariant is the beta distribution. As a test of the results, Monte Carlo simulations are used to evolve numerical ensembles toward the invariant Dirichlet distribution.

1. Objective

We develop a Fokker-Planck equation whose statistically stationary solution is the Dirichlet distribution [13]. The system of stochastic differential equations (SDE), equivalent to the Fokker-Planck equation, yields a Markov process that allows a Monte Carlo method to numerically evolve an ensemble of fluctuating variables that satisfy a unit-sum requirement. A Monte Carlo solution is used to verify that the invariant distribution is Dirichlet.

The Dirichlet distribution is a statistical representation of nonnegative variables subject to a unit-sum requirement. The properties of such variables have been of interest in a variety of fields, including evolutionary theory [4], Bayesian statistics [5], geology [6, 7], forensics [8], econometrics [9], turbulent combustion [10], and population biology [11].

2. Preview of Results

The Dirichlet distribution [13] for a set of scalars , , , is given by where are parameters, , and denotes the gamma function. We derive the stochastic diffusion process, governing the scalars, , where is an isotropic vector-valued Wiener process [12], and , , and are coefficients. We show that the statistically stationary solution of (2) is the Dirichlet distribution, (1), provided that the SDE coefficients satisfy The restrictions imposed on the SDE coefficients, , , and , ensure reflection towards the interior of the sample space, which is a generalized triangle or tetrahedron (more precisely, a simplex) in dimensions. The restrictions together with the specification of the Nth scalar as ensure Indeed, inspection of (2) shows that, for example, when , the diffusion is zero and the drift is strictly positive, while if , the diffusion is zero () and the drift is strictly negative.

3. Development of the Diffusion Process

The diffusion process (2) is developed by the method of potential solutions.

We start from the Itô diffusion process [12] for the stochastic vector, , with drift, , diffusion, , and the isotropic vector-valued Wiener process, , where summation is implied for repeated indices. Using standard methods given in [12], the equivalent Fokker-Planck equation governing the joint probability, , derived from (5), is with diffusion . Since the drift and diffusion coefficients are time homogeneous, and , (5) is a statistically stationary process and the solution of (6) converges to a stationary distribution, [12, Sec. ]. Our task is to specify the functional forms of and so that the stationary solution of (6) is , defined by (1).

A potential solution of (6) exists if is satisfied, [12, Sec. ]. Since the left-hand side of (7) is a gradient, the expression on the right must also be a gradient and can therefore be obtained from a scalar potential denoted by . This puts a constraint on the possible choices of and and on the potential, as must also be satisfied. The potential solution is Now functional forms of and that satisfy (7) with are sought. The mathematical constraints on the specification of and are as follows.(1) must be symmetric positive semidefinite. This is to ensure the following.(i)The square-root of (e.g., the Cholesky-decomposition, ) exists, required by the correspondence of the SDE (5) and the Fokker-Planck equation (6).(ii)Equation (5) represents a diffusion.(iii), required by the existence of the inverse in (7).(2)For a potential solution to exist (7) must be satisfied. With (8) shows that the scalar potential must be It is straightforward to verify that the specifications satisfy the above mathematical constraints, and . Here , , , and . Summation is not implied for (9)–(11).

Substituting (9)–(11) into (7) yields a system with the same functions on both sides with different coefficients, yielding the correspondence between the coefficients of the Dirichlet distribution, (1), and the Fokker-Planck equation (6) with (10)-(11) as For example, for one has and from (9) the scalar potential is The system in (7) then becomes which shows that by specifying the parameters, , of the Dirichlet distribution as the stationary solution of the Fokker-Planck equation (6) with drift (10) and diffusion (11) is for . The above development generalizes to variables, yielding (12) and reduces to the beta distribution, a univariate specialization of for , where and , see [13].

If (12) hold, the stationary solution of the Fokker-Planck equation (6) with drift (10) and diffusion (11) is the Dirichlet distribution, (1). Note that (10)-(11) are one possible way of specifying a drift and a diffusion to arrive at a Dirichlet distribution; other functional forms may be possible. The specifications in (10)-(11) are a generalization of the results for a univariate diffusion process, discussed in [13, 14], whose invariant distribution is beta.

The shape of the Dirichlet distribution, (1), is determined by the coefficients, . Equation (12) shows that in the stochastic system, different combinations of , , and may yield the same and that not all of , , and may be chosen independently to keep the invariant Dirichlet.

4. Corroborating That the Invariant Distribution Is Dirichlet

For any multivariate Fokker-Planck equation there is an equivalent system of Itô diffusion processes, such as the pair of (5)-(6) [12]. Therefore, a way of computing the (discrete) numerical solution of (6) is to integrate (5) in a Monte Carlo fashion for an ensemble [15]. Using a Monte Carlo simulation we show that the statistically stationary solution of the Fokker-Planck equation (6) with drift and diffusion (10)-(11) is a Dirichlet distribution, (1).

The time evolution of an ensemble of particles, each with variables (), is numerically computed by integrating the system in (5), with drift and diffusion (10)-(11), for as for each particle . In (18)-(19) and are independent Wiener processes, sampled from Gaussian streams of random numbers with mean and covariance . 400,000 particle triplets, , are generated with two different initial distributions, displayed in Figures 1(a) and 2(a), a triple-delta and a box, respectively. Each member of both initial ensembles satisfy . Equations (18)–(20) are advanced in time with the Euler-Maruyama scheme [16] with time step . Table 1 shows the coefficients of the stochastic system (18)–(20), the corresponding parameters of the final Dirichlet distribution, and the first two moments at the initial times for the triple-delta initial condition case. The final state of the ensembles is determined by the SDE coefficients, constant for these exercises, also given in Table 1, the same for both simulations, satisfying (17).

The time evolutions of the joint probabilities are extracted from both calculations and displayed at different times in Figures 1 and 2. At the end of the simulations two distributions are plotted in Figures 1(d) and 2(d): the one extracted from the numerical ensemble and the Dirichlet distribution determined analytically using the SDE coefficients—in excellent agreement in both figures. The statistically stationary solution of the developed stochastic system is the Dirichlet distribution.

For a more quantitative evaluation, the time evolutions of the first two moments, are also extracted from the numerical simulation with the triple-delta-peak initial condition as ensemble averages and displayed in Figures 3 and 4. The figures show that the statistics converge to the precise state given by the Dirichlet distribution that is prescribed by the SDE coefficients, see Table 1.

The solution approaches a Dirichlet distribution, with nonpositive covariances [2], in the statistically stationary limit, Figure 4(b). Note that during the evolution of the process, , the solution is not necessarily Dirichlet, but the stochastic variables sum to one at all times. The point (, ), governed by (18)-(19), can never leave the -dimensional (here ) convex polytope and by definition . The rate at which the numerical solution converges to a Dirichlet distribution is determined by the vectors and .

The above numerical results confirm that starting from arbitrary realizable ensembles, the solution of the stochastic system converges to a Dirichlet distribution in the statistically stationary state, specified by the SDE coefficients.

5. Relation to Other Diffusion Processes

It is useful to relate the Dirichlet diffusion process, (2), to other multivariate stochastic diffusion processes with linear drift and quadratic diffusion.

A close relative of (2) is the multivariate Wright-Fisher (WF) process [11], used extensively in population and genetic biology, where is Kronecker’s delta, with defined in (1), and . Similar to (2), the statistically stationary solution of (23) is the Dirichlet distribution [17]. It is straightforward to verify that its drift and diffusion also satisfy (7) with ; that is, WF is a process whose invariant is Dirichlet and this solution is potential. A notable difference between (2) and (23), other than the coefficients, is that the diffusion matrix of the Dirichlet diffusion process is diagonal, while that of the WF process is full.

Another process similar to (2) and (23) is the multivariate Jacobi process, used in econometrics, of Gourieroux and Jasiak [9] with , , , and .

In the univariate case, the Dirichlet, WF, and Jacobi diffusions reduce to see also [13], whose invariant is the beta distribution, which belongs to the family of Pearson diffusions, discussed in detail by Forman and Sørensen [14].

6. Summary

The method of potential solutions of Fokker-Planck equations has been used to derive a transport equation for the joint distribution of fluctuating variables. The equivalent stochastic process, governing the set of random variables, , , , reads as where and , , and are parameters, while is an isotropic Wiener process with independent increments. Restricting the coefficients to , , and and defining as above ensure and that individual realizations of ( are confined to the ()-dimensional convex polytope of the sample space. Equation (26) can therefore be used to numerically evolve the joint distribution of fluctuating variables required to satisfy a conservation principle. Equation (26) a coupled system of nonlinear stochastic differential equations whose statistically stationary solution is the Dirichlet distribution, (1), provided that the coefficients satisfy

In stochastic modeling, one typically begins with a physical problem, perhaps discrete, then derives the stochastic differential equations whose solution yields a distribution. In this paper we reversed the process: we assumed a desired stationary distribution and derived the stochastic differential equations that converge to the assumed distribution. A potential solution form of the Fokker-Planck equation was posited, from which we obtained the stochastic differential equations for the diffusion process whose statistically stationary solution is the Dirichlet distribution. We have also made connections to other stochastic processes, such as the Wright-Fisher diffusions of population biology and the Jacobi diffusions in econometrics, whose invariant distributions possess similar properties but whose stochastic differential equations are different.

Acknowledgments

It is a pleasure to acknowledge a series of informative discussions with J. Waltz. This work was performed under the auspices of the U.S. Department of Energy under the Advanced Simulation and Computing Program.