Research Article  Open Access
J. Bakosi, J. R. Ristorcelli, "A Stochastic Diffusion Process for the Dirichlet Distribution", International Journal of Stochastic Analysis, vol. 2013, Article ID 842981, 7 pages, 2013. https://doi.org/10.1155/2013/842981
A Stochastic Diffusion Process for the Dirichlet Distribution
Abstract
The method of potential solutions of FokkerPlanck 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 unitsum 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 WrightFisher 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 FokkerPlanck equation whose statistically stationary solution is the Dirichlet distribution [1–3]. The system of stochastic differential equations (SDE), equivalent to the FokkerPlanck equation, yields a Markov process that allows a Monte Carlo method to numerically evolve an ensemble of fluctuating variables that satisfy a unitsum 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 unitsum 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 [1–3] 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 vectorvalued 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 vectorvalued Wiener process, , where summation is implied for repeated indices. Using standard methods given in [12], the equivalent FokkerPlanck 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 lefthand 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 squareroot of (e.g., the Choleskydecomposition, ) exists, required by the correspondence of the SDE (5) and the FokkerPlanck 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 FokkerPlanck 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 FokkerPlanck 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 FokkerPlanck 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 FokkerPlanck 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 FokkerPlanck 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 tripledelta and a box, respectively. Each member of both initial ensembles satisfy . Equations (18)–(20) are advanced in time with the EulerMaruyama 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 tripledelta 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).

(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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 tripledeltapeak 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.
(a)
(b)
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 WrightFisher (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 FokkerPlanck 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 FokkerPlanck 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 WrightFisher 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.
References
 N. L. Johnson, “An approximation to the multinomial distribution some properties and applications,” Biometrika, vol. 47, no. 12, pp. 93–102, 1960. View at: Google Scholar
 J. E. Mosimann, “On the compound multinomial distribution, the multivariatedistribution, and correlations among proportions,” Biometrika, vol. 49, no. 12, pp. 65–82, 1962. View at: Google Scholar
 S. Kotz, N. L. Johnson, and N. Balakrishnan, Continuous Multivariate Distributions: Models and Applications, Wiley Series in Probability and Statistics: Applied Probability and Statistics, Wiley, 2000.
 K. Pearson, “Mathematical contributions to the theory of evolution. On a form of spurious correlation which may arise when indices are used in the measurement of organs,” Royal Society of London Proceedings Series I, vol. 60, pp. 489–498, 1896. View at: Google Scholar
 C. D. M. Paulino and C. A. de Bragança Pereira, “Bayesian methods for categorical data under informative general censoring,” Biometrika, vol. 82, no. 2, pp. 439–446, 1995. View at: Publisher Site  Google Scholar
 F. Chayes, “Numerical correlation and petrographic variation,” The Journal of Geology, vol. 70, pp. 440–452, 1962. View at: Google Scholar
 P. S. Martin and J. E. Mosimann, “Geochronology of pluvial lake cochise, Southern Arizona, [part] 3, pollen statistics and pleistocene metastability,” American Journal of Science, vol. 263, no. 4, pp. 313–358, 1965. View at: Publisher Site  Google Scholar
 K. Lange, “Applications of the Dirichlet distribution to forensic match probabilities,” Genetica, vol. 96, no. 12, pp. 107–117, 1995. View at: Publisher Site  Google Scholar
 C. Gourieroux and J. Jasiak, “Multivariate Jacobi process with application to smooth transitions,” Journal of Econometrics, vol. 131, no. 12, pp. 475–505, 2006. View at: Publisher Site  Google Scholar
 S. S. Girimaji, “Assumed betapdf model for turbulent mixing: validation and extension to multiple scalar mixing,” Combustion Science and Technology, vol. 78, no. 4, pp. 177–196, 1991. View at: Google Scholar
 M. Steinrucken, Y. X. Rachel Wang, and Y. S. Song, “An explicit transition density expansion for a multiallelic WrightFisher diffusion with general diploid selection,” Theoretical Population Biology, vol. 83, pp. 1–14, 2013. View at: Google Scholar
 C. W. Gardiner, Stochastic Methods, A Handbook For the Natural and Social Sciences, Springer, Berlin, Germany, 4th edition, 2009.
 J. Bakosi and J. R. Ristorcelli, “Exploring the beta distribution in variabledensity turbulent mixing,” Journal of Turbulence, vol. 11, no. 37, pp. 1–31, 2010. View at: Google Scholar
 J. L. Forman and M. Sørensen, “The Pearson diffusions: a class of statistically tractable diffusion processes,” Scandinavian Journal of Statistics, vol. 35, no. 3, pp. 438–465, 2008. View at: Publisher Site  Google Scholar
 S. B. Pope, “PDF methods for turbulent reactive flows,” Progress in Energy and Combustion Science, vol. 11, no. 2, pp. 119–192, 1985. View at: Google Scholar
 P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Third corrected printing, Springer, Berlin, Germany, 1999.
 S. Karlin and H. M. Taylor, A Second Course in Stochastic Processes, vol. 2, Academic Press, 1981.
Copyright
Copyright © 2013 J. Bakosi and J. R. Ristorcelli. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.