Multi-Gaussian Random Variables for Modeling Optical Phenomena

A generalization of the classic Gaussian random variable to the family of multi-Gaussian (MG) random variables characterized by shape parameter M > 0, in addition to the mean and the standard deviation, is introduced. The probability density function (PDF) of the MG family members is an alternating series of Gaussian functions with suitably chosen heights and widths. In particular, for integer values of M, the series has a finite number of terms and leads to flattened profiles, while reducing to the classic Gaussian PDF for M = 1. For non-integer, positive values of M, a convergent infinite series of Gaussian functions is obtained that can be truncated in practical problems. For all M > 1, the MG PDF has flattened profiles, while for 0 < M < 1, the MG PDF has cusped profiles. Moreover, the multivariate extension of the MG random variable is obtained and the log-multi-Gaussian random variable is introduced. In order to illustrate the usefulness of these new random variables for optics, the application of MG random variables to the characterization of novel speckle fields is discussed, both theoretically and via numerical simulations.


Introduction
The most famous probability density function (PDF) of a continuous random variable-Gaussianstemming from the early works of de Moivre [1] and Gauss [2] can be generalized in a number of ways for inclusion of desired shape details such as flattening, skewing, splitting, etc. One well-known generalization was proposed by Subbotin [3] (see also Ref. [4]), who extended the PDF curve to flatter or sharper versions by varying the power law of the exponential function, hence the family is sometimes termed exponential-power or super-Gaussian. Subbotin's PDF was later rescaled by Lunetta [5] and the resulting family has been well explored (cf. Reference [6]). While originally used for analysis of astrophysical data, super-Gaussian random variables are currently employed for characterizing a wide range of statistical phenomena: big data analysis [7], finance [8,9], genetics [10], and scientific impact assessment [11], to name a few. However, this seemingly transparent generalization often relies on the use of special functions, such as hypergeometric functions, as is the case for evaluation of its characteristic function [12] (see also Ref. [13]).
In this paper, we introduce a novel family of continuous random variables that serves a similar purpose as the super-Gaussian family, i.e., it reshapes the Gaussian distribution to flat-top or cusped versions, depending on the value of the shape parameter. However, the main advantage of our family, vice the super-Gaussian, stems from the fact that the statistical properties of its members can be expressed as series of those for Gaussian random variables, with very simple expressions defining their heights and widths. Moreover, in the case when the shape parameter is an integer, flat-top distributions can be formed by a finite number of terms in the series. As we show, the ability to represent a PDF of the new random variable as a linear combination of Gaussian contributions leads to unprecedented tractability in the derivation of a number of its characteristics. For any value of the shape parameter, we will term our new family of random variables multi-Gaussian (MG), not to be confused with the well-known multivariate Gaussian. The MG functions of various dimensions have been previously used in optics for modeling of beam intensity profiles [14], aperture shapes [15], scattering potentials [16,17], and various correlation and coherence functions [18][19][20][21].
Starting with the finite series case (i.e., integer values of the shape parameter), we first derive the MG PDF and discuss its characteristics to gain insight into the behavior of MG random variables. We then provide calculations of the MG cumulative distribution function (CDF), moment generating function (MGF), characteristic function (CF), and the cumulant generating function (CGF). We derive the general expressions for the moments and the cumulants of any order and find explicit expressions for the first four members of each sequence. Next, we introduce the log-multi-Gaussian (LMG) random variable on assuming that its logarithm is MG-distributed. Such a distribution can be regarded as an extension of the classic log-normal distribution [22] to flat-topped profiles. The LMG distribution has an analog derived from the super-Gaussian family [23,24]; however, as we show, the calculations of major statistical characteristics of LMG random variables can be obtained almost effortlessly, as linear combinations of well-known results for log-normal variables. We also discuss the natural extension of the univariate MG random variable to the multivariate domain and, in particular, the bivariate domain. Such variables reduce to classic multivariate/bivariate Gaussian random variables if the shape parameter M = 1. We also show how to generalize all the aforementioned results to an MG distribution with shape parameter M taking on any positive value, not necessarily an integer, while discussing numerical examples for the special case when M is the reciprocal of an integer, leading to the formation of various cusped distributions. The features of LMG and bivariate MG random variables with any positive M are also briefly outlined and the corresponding numerical results are provided.
Lastly, as an example, we apply the aforementioned analysis on MG random variables to an optical application. Namely, we derive the first-order statistics (single-spatial-point statistics) of polarized, or scalar MG speckle fields. Analogous to traditional speckle fields, where the optical field real and imaginary parts are jointly Gaussian [25], MG speckle fields have real and imaginary parts that are jointly multi-Gaussian. From the bivariate MG PDF, we derive the joint PDF of the MG speckle amplitude and phase, from which, we further obtain the amplitude, phase, and intensity marginal PDFs. We describe how to simulate MG speckle fields; our simulation method can easily be applied to synthesize MG speckle fields using a spatial light modulator. Then we validate our MG speckle analysis by comparing Monte Carlo moments to the corresponding theory.

Probability density function and cumulative distribution function
Let us begin by recalling that the Gaussian PDF of a continuous real random variable X with mean µ (location parameter) and standard deviation σ (scale parameter), i.e., X ∼ N (µ, σ), has the form Consider now a function where M is a positive real number. For M = 1, f (x) reduces to the Gaussian function in Eq. (1), while for M>1 and 0<M<1 it describes flat-topped and cusped distributions, respectively. Hence, we may refer to M as a shape parameter.
We will first discuss the case when M is an integer, M = 1, 2, 3, . . .. On using the binomial theorem, namely, , we arrive at the finite series with equal mean µ for all terms and m th -term standard deviation Let us now introduce a PDF as and therefore, we finally obtain the multi-Gaussian PDF: We may also say that X ∼ N (µ, σ, M) as a generalization of a normal random variable to any value of the shape parameter M. The cumulative distribution function (CDF) of the MG random variable in Eq. (8) can be readily calculated from its definition and by changing the order of summation and integration: where Erf stands for the error function. Figure 1(A) and (B) show the PDF and CDF of a MG random variable for M = 1, 2, 10, 40 calculated from Eqs. (8) and (9), respectively. Figure 2 shows the PDF for M = 10 and different values of µ and σ.

Moment generating function, characteristic function, and moments
The moment generating function (MGF) of the MG random variable in Eq. (8) can also be directly evaluated from its definition:  The characteristic function (CF) of the MG distribution can be readily obtained either from its definition or its relation to the MGF: The statistical moment of order k can be evaluated from the MGF in Eq. (10) via the expression or directly via the PDF function in Eq. (8). Following the latter path, we find that where we recognize that each integral term is the k th moment of the Gaussian distribution [see Eq.
(1)] with mean µ and standard deviation σ m : The U(a, b, z) in the above expression is the confluent hypergeometric function [26]. In particular, on substituting the first four moments of the Gaussian distribution from Eq. (14) into Eq. (13), we find at once that µ (MG) The parameters ξ n (M) are defined via ratios where Also, ξ n (1) = 1 for any n = 1, 2, 3, . . ., hence, the sequence in Eq. (15) reduces to that in Eq. (14).

Cumulant generating function and cumulants
The cumulant generating function (CGF) of an MG random variable has the form as implied by Eq. (10). On expanding the exponential function in Eq. (18) in a Taylor series, interchanging the order of the two summations, and recognizing the coefficients in the sum as ξ n (M), we get Further, expanding the logarithmic function in Eq. (19) in a Taylor series, we arrive at the double power series The cumulants can be found as coefficients κ (MG) n of powers of h in Eq. (20): In particular, the first four cumulants are: with all odd-order cumulants greater than three equal to zero.

Log-multi-Gaussian (LMG) distribution
Let random variable X be MG-distributed, and Y is ln(Y) = X. If we let Y = exp(X) = g(X) and Substitution of the MG PDF from Eq. (8) into Eq. (23) leads to which we call the log-multi-Gaussian (LMG) PDF. For M = 1, the LMG PDF reduces to the classic log-normal (LN) PDF [22]. The CDF of a LMG random variable can be readily found from The MGF of the LMG distribution diverges but the moments can be found from the definition: The expression in the braces is the k th moment of the LN distribution with mean µ and standard deviation σ m : thus, In particular, the first four moments are Figure 3 shows the PDF and CDF of a LMG random variable with µ = 0 and σ = 1, for several values of the index M. As M increases, the PDF profiles become sharper with maxima occurring at smaller values of y. Figure 4 illustrates the PDF of a LMG random variable with M = 10 but different values of σ and µ.

Multivariate multi-Gaussian distribution
Let X = [X 1 , X 2 , . . . , X N ] T ∈ R be an N-dimensional vector of real random variables with T standing for the transpose. Also, let the N-dimensional vector of their mean values be Then, the PDF of the multivariate MG random vector can be defined by a straightforward generalization of the multivariate Gaussian PDF as where det stands for determinant of a matrix and superscript −1 denotes matrix inverse. In particular, for N = 2 (i.e., the bivariate MG case), Eq. (30) simplifies to To arrive at Eq. (31), we used

Multi-Gaussian moment theorem
Using Eq. (30), we can derive the MG version of the Gaussian moment theorem [25]. Assuming µ = 0 for simplicity, the cross-correlation of N MG random variables is The N-dimensional integral evaluates to where P is the set of distinct pair groupings of the N variables and ρ ij is the cross-correlation coefficient of x i and x j [25]. Note that if N is odd, ⟨x 1 x 2 · · · x N ⟩ = 0.

Generalization to any M>0
Let us now return to Eq. (2) and assume that M is any positive number, not necessarily an integer. The fractional binomial theorem is now applicable, and Eq. (3) becomes where is the Pochhammer symbol (falling factorial). We now set u = 1 and ]︁ and express f (x) in Eq. (2) via the generalized binomial series: The MG PDF can be obtained as the normalized version of f (x) in Eq. (40) [see Eq. (6)]: with The infinite series in Eq. (41) converges as m → ∞, and hence, in practical applications can be suitably truncated. All calculations relating to the CDF, MGF, CF, KGF, moments, cumulants, etc. can be carried out in the same manner as above (integer M) except with binomial coefficients of the form in Eqs. (38) and (39). For instance, the CDF for non-integer M>0 takes the form instead of that given in Eq. (17). The expression for the LMG PDF now takes the form where C 0 (M) is given in Eq. (42), and the bivariate MG PDF becomes where z is the same as in Eq. (32).

Theory for speckle fields
As an example optical application of MG random variables, we derive the first-order statistics of amplitude, phase, and intensity of polarized MG speckle fields. We begin the analysis with the joint PDFs of the real R and imaginary I parts of an MG speckle field. We assume that the means ⟨R⟩ = ⟨I⟩ = 0 hereafter.

Rectangular MG speckle fields
If R and I are statistically independent (analogous to fields obeying circular Gaussian statistics), then the joint PDF is the product of the marginal R and I PDFs and equal to where M and N are the shape parameters of the MG distributions for R and I, respectively, σ Rm = σ R / √ m, σ In = σ I / √ n, and C 0 (M) is a normalization factor given generally in Eq. (44). As discussed above, M and N can be any real numbers greater than zero. When either is an integer, the corresponding series is finite with all terms greater than M or N equal to zero. If M = N = 1 and σ R = σ I , Eq. (47) simplifies to the joint PDF of two independent Gaussian random variables. This PDF, having a circular shape, describes fully developed speckle fields [also called circular complex Gaussian (CCG) fields] generated by scattering spatially coherent light from rough surfaces or diffusers [25]. In this way, Eq. (47) can be considered a generalization of CCG speckle fields, although its shape is rectangular for M, N ≠ 1. By correlating R and I, we can develop a p R,I that is circular or, more generally, elliptical in shape. We present that theory in the Section 4.2.
We note that there is no single joint MG PDF that incorporates both rectangular-and ellipticalshaped p R,I . The one exception to this statement is when M = N = 1, but this is just the joint Gaussian PDF which has already been discussed.
The joint PDF of the amplitude and phase can be found by performing a random variable transformation on Eq. (47), namely, where R = A cos ϕ and I = A sin ϕ [25]. Evaluating Eq. (48) produces The marginal PDFs of ϕ and A are respectively. Substituting Eq. (49) into Eq. (50) and evaluating the integral using substitution yields The integral in Eq. (51) is more difficult to evaluate; the details are shown in Appendix A. Equation (51) simplifies to where I 0 (x) is a zeroth-order modified Bessel function of the first kind. In contrast to CCG fields, here p A,φ (A, ϕ) ≠ p A (A) p φ (ϕ). Thus, statistically independent MG R and I do not produce speckle fields with independent amplitudes and phases.
Lastly, the PDF of the speckle intensities can be derived via a random variable transformation of p A (A), i.e., We investigate these PDFs further in Section 5. The speckle contrast, i.e., C = σ I /⟨I⟩ [25], can be numerically computed from Eq. (55); however, we can obtain an analytical result using the MG moments given in Eq. (15). We begin by expressing C as The mean and variance of intensity are Referring back to Eq. (15), the moments of R and I are Substituting Eqs. (57) and (58) into Eq. (56) and simplifying yields If σ R = σ I and M = N (i.e., a "square" MG speckle field), Eq. (59) simplifies to This further simplifies to C = 1 when M = 1, as expected for CCG fields. Figure 9 shows a plot of Eq. (60). This figure can be used to produce square MG speckle fields with desired contrasts. As previously stated, square MG speckle fields are equal to CCG speckle fields when M = 1. All M<1 yield speckle fields with C>1. For M>1, we obtain C<1.

Elliptical MG speckle fields
Having derived the speckle phase, amplitude, and intensity PDFs for independent MG R and I, we now consider MG R and I that are correlated. The joint PDF of R and I takes the form where ρ is the Gaussian distribution correlation coefficient and C 0 (M) is given in Eq. (44). In general, the shape of p R,I is elliptical. If ρ = 0 and σ R = σ I , p R,I becomes circular; nevertheless, R and I are not statistically independent unless M = 1. The joint PDF of the amplitude and phase can be found using Eq. (48). Substituting Eq. (61) into Eq. (48) and simplifying yields Like above, we are interested in the phase and amplitude marginal PDFs defined in Eqs. (50) and (51), respectively. The phase PDF integral can be evaluated using substitution, namely, and interestingly, does not depend on the MG shape parameter M. The amplitude PDF integral is again much more difficult to evaluate. The details are shown in Appendix B: With p A , the PDF of intensity follows immediately from Eq. (54): Inspection of Eqs. (63) and (64) reveals that, in general, p A,φ (A, ϕ) ≠ p A (A) p φ (ϕ) and therefore, A and ϕ are not statistically independent. However, if ρ = 0 and σ R = σ I , i.e., Eq. (61) is circular, A and ϕ are independent.
We can obtain an expression for the speckle contrast C using Eq. (15) and the MG moment theorem discussed in Section 3.3. The mean intensity and intensity variance for elliptical MG speckle fields take the form The only difference between Eq. (66) and Eq. (57) is the term in parentheses. For rectangular MG speckle fields this term is zero as R and I are statistically independent. This is not the case here, however. The single moments in Eq. (66) are defined in Eq. (58) with N = M. We obtain the joint moment by using the MG moment theorem and Eq. (37), namely, Substituting Eqs. (58), (66), and (67) into Eq. (56) and simplifying yields the desired result: If σ R = σ I , Eq. (68) simplifies to .
(69) Figure 10 shows a plot of Eq. (69). Like the associated square MG speckle contrast plot (Fig. 9), Fig. 10 can be used to design speckle fields with specific contrasts. Note that Eq. (69) is the square MG speckle contrast [Eq. (60)] multiplied by a term containing ρ. Because of this ρ term, elliptical MG speckle fields have a larger range of possible contrasts than square MG speckle fields. In the next section, we generate MG speckle fields and compare the Monte Carlo PDFs to the theoretical expressions derived above.

Simulation
Here, we present simulation results of polarized MG speckle fields. First, we discuss the simulation details.

Setup
For the simulations, we used grids that were 256 points per side with a side-length equal to 1 m. We generated 500 each rectangular and elliptical MG field instances, from which, we computed p R,I , p R , p I , p A , p φ , and p I . The simulated rectangular MG parameters were M = 20, N = 1/20, σ R = √︁ 1/2, and σ I = 1; the elliptical MG parameters were M = 1/20, σ R = √︁ 1/2, σ I = 1, and ρ = 0.6.
We used the following process to generate rectangular and elliptical MG speckle fields: 1. Using the spectral representation method [27][28][29], we "colored" 256 × 256 matrices of zero-mean, unit-variance, independent Gaussian random numbers (representing R and I) to simulate partial spatial coherence and ultimately control the speckle size. We used a Gaussian-shaped correlation function for this purpose, viz., where δ = 5 cm.
2. To generate rectangular and elliptical MG fields, we used a technique known by several different names in the literature: translation [30,31], NORTA (NORmal To Anything) [32], and the inverse CDF transform method [28,29]. Translation (NORTA, etc.) starts with a correlated Gaussian-distributed sequence of random numbers and maps it, via a non-linear transformation, to a correlated non-Gaussian sequence.
(a) For rectangular MG fields, the input Gaussian sequences R G and I G were translated to MG R and I using the following: where F G is the Gaussian CDF and F −1 is the inverse MG CDF [the forward CDF is given in Eq. (9)]. In Eq. (71), F G is parameterized by the mean and standard deviation of the input sequence (either R G or I G ). The output of F G is a standard uniform random sequence, which subsequently serves as the input to F −1 . The inverse MG CDF (computed numerically) is parameterized by the desired MG shape parameter and Gaussian standard deviation. The output of Eq. (71) is two independent sequences of MG random numbers with shape parameters M and N and Gaussian standard deviations σ R and σ I , respectively.
(b) Since elliptical MG R and I are correlated, the translation process was more complicated. We first converted R G and I G to A G and ϕ G . Then, we applied the following sequence of non-linear transforms: where F U and F R are the uniform and Rayleigh CDFs, respectively. F U and F R are parameterized by the interval and scale parameter of the respective distributions, and the outputs of both are standard uniform sequences. These serve as inputs to F −1 φ and F −1 A |φ , which are the inverse MG phase CDF [computed from Eq. (63)] and the inverse conditional (A given ϕ) CDF, respectively. The latter is computed numerically from the corresponding forward conditional CDF: where p A,φ and p φ are given in Eqs. (62) and (63). In the first step of Eq. (72), a ϕ sequence is generated. F −1 φ is parameterized by the desired MG σ R , σ I , and ρ. Since A depends on ϕ in elliptical MG fields (A and ϕ are correlated), the ϕ sequence must be input into the second step of Eq. (72). F −1 A |φ is parameterized by the desired MG shape parameter M, σ R , σ I , and ρ. The output of Eq. (72) is an MG ϕ sequence and a correlated A sequence.
3. Lastly, we formed rectangular or elliptical MG speckle field realizations by U = R + jI or U = A exp (jϕ), respectively.

Results
Figures 11-13 show the results for the rectangular MG speckle field simulation. Figure 11 compares rectangular MG and CCG speckle field realizations. Figure 12 reports the R and I The agreement between the theoretical and Monte Carlo PDFs in these figures is excellent. It should be noted that although MG speckle fields look qualitatively similar to CCG speckle fields (Figs. 11 and 14), their statistics are quite different, especially in phase (Figs. 12, 13,  15, and 16). There has been recent theoretical and experimental work on generating speckle fields with customizable non-Rayleigh amplitude statistics for use in microscopy and imaging applications [33][34][35]. Although we presented MG speckle fields as an example application of MG random variables, the theory and synthesis method developed in this section might find use in engineered-speckle-pattern applications.  In the legend, "G" stands for Gaussian, "MG" for multi-Gaussian, "Thy" for theory, and "Sim" for simulation. Fig. 13. A, ϕ, and I PDFs-(A) p A theory versus simulation, (B) p φ theory versus simulation, and (C) p I theory versus simulation. In the legend, "G" stands for Gaussian, "MG" for multi-Gaussian, "Thy" for theory, and "Sim" for simulation.  In the legend, "G" stands for Gaussian, "MG" for multi-Gaussian, "Thy" for theory, and "Sim" for simulation. Fig. 16. A, ϕ, and I PDFs-(A) p A theory versus simulation, (B) p φ theory versus simulation, and (C) p I theory versus simulation. In the legend, "G" stands for Gaussian, "MG" for multi-Gaussian, "Thy" for theory, and "Sim" for simulation.

Summary
We introduced a new family of continuous, real random variables whose PDF represents flattened and cusped deviations from the Gaussian PDF depending on a single shape parameter taking positive values. Gaussian random variables are a particular case of our family when the shape parameter equals unity. While our new PDF can be expressed in closed form, it is possible and much more analytically convenient to express it as a linear combination of Gaussian PDFs with alternating signs, equal means, and monotonically decreasing standard deviations. Hence, the name of our new random variable, multi-Gaussian. It was shown that for integer values of the shape parameter, the series of Gaussian functions has a finite number of terms, while for non-integer positive values, the series has infinite terms, but is convergent and can be suitably truncated. Due to this feature, the calculations of the statistical properties of multi-Gaussian random variables reduce to algebraic operations over the well-known statistics of Gaussian random variables. In addition, we developed the log-multi-Gaussian random variable and extended univariate multi-Gaussian random variables to their multivariate or joint counterparts. Other straightforward extensions can readily be made, for instance, relating to complex multi-Gaussian random variables.
Our motivation for introducing the new type of the PDF function was driven by a search for a tractable analytical model to be used for prescribing the single-point statistics of a scalar optical speckle field with a deformed (flattened, stretched, or skewed) shape in the central region of the distribution, and a Gaussian-like edge profile. Such fields are quite common in optics and might appear in cases when filtering (absorption or gain) of the field is not uniform but rather depends on its local amplitude or phase values at incidence. In addition, digitally controlled devices such as spatial light modulators are now capable of practically unlimited manipulation of the amplitude and phase statistics of laser beams. This experimental field seems however little explored because of the lack of analytical models for speckle PDFs. In particular, to our knowledge, a simple PDF model describing the aforementioned situation does not exist in the literature. We have therefore used the introduced MG PDF (with integer values of the shape parameter) for modeling of the real and imaginary parts of a speckle field and found that two PDF families were possible, the rectangular and elliptical versions. These families can be regarded as two extensions of the classic circular-Gaussian PDF, both reducing to it when the shape parameter equals unity. We then investigated in detail the first-order statistics of these scalar multi-Gaussian speckle fields for both families, reporting the expressions for the amplitude, phase and intensity PDFs. In addition, for both families, we obtained simple analytic expressions for the high-order speckle statistics, e.g., the speckle contrast or equivalently, the scintillation index. We further simulated the generation of such random optical fields and compared the Monte-Carlo-computed amplitude, phase, and intensity PDFs to the corresponding theory. The agreement was found to be excellent, thereby validating our theory and synthesis approach.
We envision the multi-Gaussian family of random variables to be useful in the same applications as Subbotin's family (exponential-power or super-Gaussian family), but in situations where simple, analytical results are of importance. As we have illustrated, this statement is valid in general and, in particular, in application to the theory of optical speckle.