A stochastic diffusion process for the Dirichlet distribution

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.


Objective
We develop a Fokker-Planck equation whose statistically stationary solution is the Dirichlet distribution [1,2,3]. 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 non-negative 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].

Preview of results
The Dirichlet distribution [1,2,3] for a set of scalars 0 ≤ Y α , α = 1, . . . , N − 1, N −1 α=1 Y α ≤ 1, is given by where ω α > 0 are parameters, Y N = 1 − N −1 β=1 Y β , and Γ(·) denotes the gamma function. We derive the stochastic diffusion process, governing the scalars, Y α , where dW α (t) is an isotropic vector-valued Wiener process [12], and b α > 0, κ α > 0, and 0 < S α < 1 are coefficients. We show that the statistically stationary solution of Eq. (2) is the Dirichlet distribution, Eq. (1), provided the SDE coefficients satisfy The restrictions imposed on the SDE coefficients, b α , κ α , and S α , ensure reflection towards the interior of the sample space, which is a generalized triangle or tetrahedron (more precisely, a simplex) in N −1 dimensions. The restrictions together with the specification of the N th scalar as Indeed, inspection of Eq. (2) shows that for example when Y 1 = 0, the diffusion is zero and the drift is strictly positive, while if Y 1 = 1, the diffusion is zero (Y N = 0) and the drift is strictly negative.

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, Y α , with drift, a α (Y), diffusion, b αβ (Y), and the isotropic vector-valued Wiener process, dW β (t), where summation is implied for repeated indices. Using standard methods given in [12] the equivalent Fokker-Planck equation governing the joint probability, F (Y, t), derived from Eq. (5), is with diffusion B αβ = b αγ b γβ . Since the drift and diffusion coefficients are time-homogeneous, a α (Y, t) = a α (Y) and B αβ (Y, t) = B αβ (Y), Eq. (5) is a statistically stationary process and the solution of Eq. (6) converges to a stationary distribution, [12] Sec. 6.2.2. Our task is to specify the functional forms of a α (Y) and b αβ (Y) so that the stationary solution of Eq. (6) is D(Y), defined by Eq. (1). A potential solution of Eq. (6) exists if is satisfied, [12] Sec. 6.2.2. Since the left hand side of Eq. (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 φ(Y). This puts a constraint on the possible choices of a α and B αβ and on the potential, as φ, αβ = φ, βα must also be satisfied. The potential solution is Now functional forms of a α (Y) and B αβ (Y) that satisfy Eq. (7) with F (Y) ≡ D(Y) are sought. The mathematical constraints on the specification of a α and B αβ are as follows: 1. B αβ must be symmetric positive semi-definite. This is to ensure that • the square-root of B αβ (e.g. the Cholesky-decomposition, b αβ ) exists, required by the correspondence of the SDE (5) and the Fokker-Planck equation (6), • Eq. (5) represents a diffusion, and • det(B αβ ) = 0, required by the existence of the inverse in Eq. (7).
2. For a potential solution to exist Eq. (7) must be satisfied.
With F (Y) ≡ D(Y) Eq. (8) shows that the scalar potential must be It is straightforward to verify that the specifications satisfy the above mathematical constraints, 1. and 2. Here b α > 0, κ α > 0, 0 < S α < 1, and Summation is not implied for Eqs. (9-11). Substituting Eqs. (9-11) into Eq. (7) yields a system with the same functions on both sides with different coefficients, yielding the correspondence between the N coefficients of the Dirichlet distribution, Eq. (1), and the Fokker-Planck equation (6) with Eqs. (10)(11) as Eq. (7) then becomes the system 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 D(Y, ω) for N = 3. The above development generalizes to N variables, yielding Eqs. (12)(13), and reduces to the beta distribution, a univariate specialization of D for N = 2, where Y 1 = Y and Y 2 = 1 − Y , see [13]. If Eqs. (12)(13) hold, the stationary solution of the Fokker-Planck equation (6) with drift (10) and diffusion (11) is the Dirichlet distribution, Eq. (1). Note that Eqs. (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 Eqs. (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, Eq. (1), is determined by the N coefficients, ω α . Eqs. (12)(13) show that in the stochastic system, different combinations of b α , S α , and κ α may yield the same ω α and that not all of b α , S α , and κ α may be chosen independently to keep the invariant Dirichlet.

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 Eqs. (5-6) [12]. Therefore, a way of computing the (discrete) numerical solution of Eq. (6) is to integrate Eq. (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, Eq. (1).
The time-evolution of an ensemble of particles, each with N = 3 variables (Y 1 , Y 2 , Y 3 ), is numerically computed by integrating the system of equations (5), with drift and diffusion (10)(11), for N = 3 as for each particle i. In Eqs. (20-21) dW 1 and dW 2 are independent Wiener processes, sampled from Gaussian streams of random numbers with mean dW α = 0 and covariance dW α dW β = δ αβ dt. 400,000 particle-triplets, (Y 1 , Y 2 , Y 3 ), are generated with two different initial distributions, displayed in the upper-left of Figures 1 and 2, a triple-delta and a box, respectively. Each member of both initial ensembles satisfy 3 α=1 Y α = 1. Eqs. (20-22) are advanced in time with the Euler-Maruyama scheme [16] with time step ∆t = 0.05. Table 1 shows the coefficients of the stochastic system (20-22), 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 are determined by the SDE coefficients, constant for these exercises, also given in Table 1, the same for both simulations, satisfying Eq. (19).
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 at the bottom-right of both figures: 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 evolution of the first two moments, µ α = Y α = , 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 non-positive covariances [2], in the statistically stationary limit, Figure 4(b). Note that during the evolution of the process, 0 < t 80, the solution is not necessarily Dirichlet, but the stochastic variables sum to one at all times. The point (Y 1 , Y 2 ), governed by Eqs. (20-21), can never leave the (N −1)-dimensional (here N = 3) convex polytope and by definition Y 3 = 1 − Y 1 − Y 2 . The rate at which the numerical solution converges to a Dirichlet distribution is determined by the vectors b α 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.

Relation to other diffusion processes
It is useful to relate the Dirichlet diffusion process, Eq. (2), to other multivariate stochastic diffusion processes with linear drift and quadratic diffusion.
A close relative of Eq. (2) is the multivariate Wright-Fisher (WF) process [11], used extensively in population and genetic biology, where δ αβ is Kronecker's delta, ω = N β=1 ω β with ω α defined in Eq. (1) and, Y N = 1− N −1 β=1 Y β . Similarly to Eq. (2), the statistically stationary solution of Eq. (23) is the Dirichlet distribution [17]. It is straightforward to verify that its drift and diffusion also satisfy Eq. (7) with F ≡ D, i.e. WF is a process whose invariant is Dirichlet and this solution is potential. A notable difference between Eqs. (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 it is full.

Summary
The method of potential solutions of Fokker-Planck equations has been used to derive a transport equation for the joint distribution of N fluctuating variables. The equivalent stochastic process, governing the set of random variables, 0 ≤ Y α , α = 1, . . . , where Y N = 1− N −1 β=1 Y β , and b α , κ α and S α are parameters, while dW α (t) is an isotropic Wiener process with independent increments. Restricting the coefficients to b α > 0, κ α > 0 and 0 < S α < 1, and defining Y N as above ensure N α=1 Y α = 1 and that individual realizations of (Y 1 , Y 2 , . . . , Y N ) are confined to the (N−1)-dimensional convex polytope of the sample space. Eq. (26) can therefore be used to numerically evolve the joint distribution of N fluctuating variables required to satisfy a conservation principle. Eq. (26) is a coupled system of nonlinear stochastic differential equations whose statistically stationary solution is the Dirichlet distribution, Eq. (1), provided the coefficients 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.    Table 1.  Tab. 1: Initial and final states of the Monte-Carlo simulation starting from a triple-delta. The coefficients, b 1 , b 2 , S 1 , S 2 , κ 1 , κ 2 , of the system of SDEs (20-22) determine the distribution to which the system converges. The Dirichlet parameters, implied by the SDE coefficients via Eqs. (17)(18)(19), are in brackets. The corresponding statistics are determined by the wellknown formulae of Dirichlet distributions [2].