International Journal of Stochastic Analysis

Volume 2013 (2013), Article ID 842981, 7 pages

http://dx.doi.org/10.1155/2013/842981

Research Article

A Stochastic Diffusion Process for the Dirichlet Distribution

J. Bakosi and J. R. Ristorcelli

Los Alamos National Laboratory, Los Alamos, NM 87545, USA

Received 19 December 2012; Accepted 1 March 2013

Academic Editor: Hong K. Xu

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.

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 [1–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 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 [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 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).

Table 1: Initial and final states of the Monte Carlo simulation starting from a triple delta. The coefficients, , , , , , , of the system of SDEs (18)–(20) determine the distribution to which the system converges. The Dirichlet parameters, implied by the SDE coefficients via (15)–(17), are in brackets. The corresponding statistics are determined by the well-known formulae of Dirichlet distributions [2].

Initial state:

triple delta,

see Figure 1
SDE coefficients and the statistics of

their implied Dirichlet distribution in

the statistically stationary state

(a)

(b)

(c)

(d)

Figure 1: Time evolution of the joint probability, , extracted from the numerical solution of (18)–(20). The initial condition is a triple-delta distribution, with unequal peaks at the three corners of the sample space. At the end of the simulation, , the solid lines are those of the distribution extracted from the numerical ensemble, and the dashed lines are those of a Dirichlet distribution to which the solution converges in the statistically stationary state, implied by the constant SDE coefficients, sampled at the same heights.

(a)

(b)

(c)

(d)

Figure 2: Time evolution of the joint probability, , extracted from the numerical solution of (18)–(20). The initial condition is a box with diffused sides. By , the distribution converges to the same Dirichlet distribution as in Figure 1.

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.

Figure 3: Time evolution of the means, extracted from the numerically integrated system (18)–(20), starting from the triple-delta initial condition. Dotted-solid lines: numerical solution, dashed lines: statistics of the Dirichlet distribution determined analytically using the constant coefficients of the SDE, see Table 1.

(a)

(b)

Figure 4: Time evolution of the second central moments, extracted from the numerically integrated system (18)–(20), starting from the triple-delta initial condition. The legend is the same as in Figure 3.

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.

References

N. L. Johnson, “An approximation to the multinomial distribution some properties and applications,” Biometrika, vol. 47, no. 1-2, pp. 93–102, 1960.

J. E. Mosimann, “On the compound multinomial distribution, the multivariate-distribution, and correlations among proportions,” Biometrika, vol. 49, no. 1-2, pp. 65–82, 1962.

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.

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.

F. Chayes, “Numerical correlation and petrographic variation,” The Journal of Geology, vol. 70, pp. 440–452, 1962.

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.

K. Lange, “Applications of the Dirichlet distribution to forensic match probabilities,” Genetica, vol. 96, no. 1-2, pp. 107–117, 1995.

C. Gourieroux and J. Jasiak, “Multivariate Jacobi process with application to smooth transitions,” Journal of Econometrics, vol. 131, no. 1-2, pp. 475–505, 2006.

S. S. Girimaji, “Assumed beta-pdf model for turbulent mixing: validation and extension to multiple scalar mixing,” Combustion Science and Technology, vol. 78, no. 4, pp. 177–196, 1991.

M. Steinrucken, Y. X. Rachel Wang, and Y. S. Song, “An explicit transition density expansion for a multi-allelic Wright-Fisher diffusion with general diploid selection,” Theoretical Population Biology, vol. 83, pp. 1–14, 2013.

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 variable-density turbulent mixing,” Journal of Turbulence, vol. 11, no. 37, pp. 1–31, 2010.

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.

S. B. Pope, “PDF methods for turbulent reactive flows,” Progress in Energy and Combustion Science, vol. 11, no. 2, pp. 119–192, 1985.

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.