Simulation versus stochastic perturbation method in symbolic computing with Gaussian random variables


  • Author: Marcin Kamiński
  • Date: 05 Aug 2014
  • Copyright: Image is copyright of Professor Kamiński


Modern uncertainty analysis in engineering includes quite various techniques like analytical calculus, traditional and modified Monte-Carlo simulation, fuzzy sets theory, polynomial chaos or Karhunen-Loeve expansions as well as, of course, a variety of perturbation methods. The most interesting issue is the availability of higher order probabilistic characteristics, reliable determination, the overall computational error and also, the time and computer power effort necessary to accomplish these goals. The second, extremely important issue is the availability of these methods for a broader spectrum of computer users with their slower and less powerful notebooks and desktops, including very basic computer software. Recommendation for the stochastic perturbation technique is justified by the fact that it can be implemented in any symbolic software, even with lower level programming languages standard procedures. It does not demand fast and powerful computers (even a very slow computer is fast enough) and the user may extend this method by himself. One may build his own perturbation-based applications in any symbolic algebra system like MAPLE, MAXIMA, MATLAB or MATHEMATICA as well as even engage Microsoft EXCEL for this purpose with the finite difference calculus.

Further, what is important for non-mathematicians, is that it does not involve larger portions on any theory except the knowledge of basic formulas from the probability theory and differentiation rules. This method has, of course, some apparent limitations like single random input source (algebraic difficulties with development of additional formulas with higher order cross-correlations), unknown precision for larger levels of the input uncertainty as well as unavailability for stochastic processes, but maybe some of them would be overcome in the nearest future.

This article shows some aspects together with numerical illustration on the fundamental problem in astronomy to give a comparison against the well-established Monte-Carlo method and to document a convergence of both methods – non-statistical and statistical – in case of the first four probabilistic characteristics. The reader interested in its further applications together with classical discrete numerical techniques like the Finite Element Method, the Boundary Element Method or the Finite Difference Method is advised to study the computational experiments contained in [2]. (Some MAPLE and MAXIMA resources concerning statistics, stochastic perturbation method as well as time series analysis with Gaussian coefficients are available in this book and may be obtained after email request directly from the author).

thumbnail image: Simulation versus stochastic perturbation method in symbolic computing with Gaussian random variables

Stochastic perturbation method

As it is known, the basic idea of the stochastic perturbation method is to expand all random variables or processes with the use of Taylor series with random coefficients about the expected values of the input uncertainty source. It is usually provided for a finite length of such an expansion like the tenth order and, for instance, it has been recently widely explored in structural engineering and safety analysis. It holds for the function f=f(b), where b represents the input uncertainty having mean value b0 and the probability distribution function (PDF) gb(x)


This expansion needs the perturbation parameter ε, taken most frequently equal to 1 and, what is definitely more important, a prior knowledge of the partial derivatives of the objective function with respect to the input random parameter. Let us note that these derivatives are calculated in a traditional, deterministic way, for the mean value of the random input. This can be achieved by numerical solution to the given boundary value problem or, alternatively, by the Least Squares Method approximation of this function coming from the series of the trial solution with varying b. Further, we derive probabilistic central moments of the function f(b) like the pth one, for instance, by inserting this expansion into its classical definition [1].


Further possible simplifications depends upon the specific type of the initial PDF and are the most transparent for the Gaussian one, where all odd order terms vanish and where higher order central moments can be replaced with specific powers of the input standard deviation. Derivation of the objective function expectation is quite straightforward and leads for the Gaussian input combined with the tenth order Taylor expansion to


One may remark that this expansion brings a single additional component in the expected value for each successive order of the expansion. Similar expansion is definitely more complex for higher central probabilistic moments of such function, even with Gaussian variables, so that it would be desirable at longer expansions to employ the computer algebra program for an automatic derivation process. An increase of the extra components for the variance is illustrated with the tenth order term (a difference in-between tenth and eight order expansions of the same variance)



where one denotes for simplicity


So that it is expected that each new order results in the new term that is with single element longer than the previous order expansion. It needs to be underlined that the iterative derivation of the particular moments brings also some negative and mixed order and moments terms that do not appear in case of the expectation. The Reader interested in some specific results of such a derivation (including additional theoretical studies related to its probabilistic convergence) may find tenth order formulas in [2], also with non-symmetric PDF of an input including both even and odd order terms. Let us remember that the overall computational error of probabilistic moment’s determination is correlated with the length of an expansion and simply cannot be too short.

Computational illustration

The perturbation method is widely explored in various scientific disciplines related to mechanics, specifically dynamics, where its deterministic version has been frequently employed. Therefore, some probabilistic problem concerning elementary astronomy [3] is numerically modelled with the use of both stochastic perturbation techniques as well as, for a comparison, classical Monte-Carlo simulation approach. It should be underlined that the statistical estimators implemented (together with skewness β and kurtosis κ) are not corrected in any way by the total number of random trials, i.e.


and follow standard equations summarized at least in [4]. Such an illustration is provided for the third Kepler law, where a distance d in-between the mass center of the star and the given planet results from the following equation:


where G denotes the gravity constant, Ms is the mass of the star, m stands for the mass of the planet and T means the period of the full rotation. The data corresponding to the Earth and the Sun have been adopted in this computational illustration treating T as the Gaussian random variable with the given expected value and the coefficient of variation adopted as the additional variable during final visualization of the results. The expected values, coefficients of variation, skewness and kurtosis of the distance d are contrasted with these two methods for the entire variability interval of , when the total number of Monte-Carlo simulation trials equals 3*105. The entire computational analysis, both in statistical and perturbation-based contexts, has been implemented and carried out in the symbolic algebra environment MAPLE, v. 14 [2,5]. Some trial towards programming of the analytical formulas for probabilistic characteristics following directly the definitions of probability theory has been also undertaken, but MAPLE returns the Laguerre polynomials as the representation for this equation roots. So that it is impossible to derive symbolically the same probabilistic characteristics in this particular case as the functions of the input coefficient of variation α(T) of the given Gaussian period.

Computational analysis enables to compare the input PDF versus the output one (Fig. 1) coming from the Monte-Carlo analysis by only as well as the expected values (Fig. 2), the coefficients of variation (Fig. 3), the skewness (Fig. 4) and also kurtosis (Fig. 5) of the computed orbit radius d. The last four graphs collected the continuous series determined with the use of the tenth order perturbation technique (tenth order polynomials of the input coefficient of variation) denoted in Figs. 2-5 by SPT and discrete sets of data coming from the estimation procedure (abbreviated by MCS). The first general impression from Fig. 1 is that both input and output PDFs are Gaussian, however a precise verification of the skewness (Fig. 4) indicates that the higher the input coefficient of variation, the more distant the output PDF from the Gaussian bell shaped curve. Successfully, stochastic perturbation technique returns continuous and smooth enough functions for all the basic characteristics with neither singularities nor local end oscillations. Further, one may conclude after Figs. 2-5 that both methods coincide very well, not only for the relatively small uncertainty margin and the first two probabilistic moments (as in its second order and second moment limited version), but for remarkable stochastic fluctuations and the first four probabilistic characteristics also.

Figure 1. Histograms of rotation period (left) versus the orbit radius (right)

Figure 2. Expected values of the Earth orbit radius

Figure 3. Coefficients of variation of the Earth orbit radius

Figure 4. Skewness of the Earth orbit radius

Figure 5. Kurtosis of the Earth orbit radius

It is interesting that the expected value of the orbit radius shows some fluctuations while increasing the input coefficient of variation and the estimated expected values (MCS solution) diverge a little from the principal trend given by the stochastic perturbation method (SPT solution). One may also notice that the increasing uncertainty in the rotation period decreases effectively the resulting expectation of the orbit and this dependence is nonlinear and convex. Further (see Fig. 3), the output coefficient of variation depends almost linearly upon the input one and as smaller from its input counterpart shows probabilistic damping in this equation (a negative probabilistic entropy variation); a linear interrelation of α(d) and α(T) corresponds very well with a similarity of both diagrams shown in Fig. 1.

Final illustration concerns the relative error analysis inherent in the Monte-Carlo simulation - it is computed for random trials and shown in Fig. 6 as a function of this N. It is discussed in the literature in addition to the expectation as well as the coefficient of variation (alternatively for variance or standard deviation) but rarely – for skewness and kurtosis.

Figure 6. Relative error analysis for the first four probabilistic characteristics

This relative error is defined as a difference of the current probabilistic characteristics and its counterpart determined for the largest number of random trials and scaled by this counterpart (it is treated as the exact result). It is apparent that the slowest statistical convergence is seen in case of the skewness which equals 0 for the Gaussian distribution unlike in case of its certain transforms. Nevertheless, satisfactory accuracy for determination of up to the fourth order probabilistic moments and characteristics in this problem with 7 digits precision using the traditional crude Monte-Carlo simulation needs about 10^7 random trials, which costs almost an hour of MAPLE operation on the laptop equipped with the processor Intel Core i7-3630QM, CPU 2.40 GHz, 7,87 GB RAM and operating system Windows 8.0.

Concluding remarks

The Stochastic perturbation method in its general order version has great potential and needs more theoretical and computational studies to develop all necessary formulas for the basic probabilistic characteristics for all possible input probability density functions, including correlated or uncorrelated multiple uncertainty sources. One could be interested in application of this apparatus to time series, where the Taylor expansion is provided with respect to both time and spatial variables and, maybe, in the area of stochastic processes. There is no doubt that extremely small time and computer power consumption associated with this approach and the fact that the additional computations may be programmed in any computer algebra system and used on any computational platform can make this method universal in the uncertainty engineering. Determination of higher order statistics in engineering is still very exceptional and the numerical illustration shows that traditional Monte-Carlo simulation needs many times more trials to determine precisely skewness and kurtosis than it is mandatory for the expectations and coefficients of variation. As one may finally notice, the stochastic perturbation method offers determination of the basic probabilistic characteristics as continuous functions of the input uncertainty level unlike the Monte-Carlo simulation, where the discrete numerical data are available. The modern stochastic perturbation method operates also on large multi-processor computers to determine probabilistic moments and their histories in nonlinear and coupled transient processes in the large scale structures.


[1] W. Feller, An introduction to probability theory and its applications. Wiley-VCH Verlag GmbH, New York, 1965.
[2] M. Kamiński, The stochastic perturbation method for computational mechanics. Wiley, Chichester, 2013.
[3] E. Belbruno, edr., New trends in astrodynamics and applications. Annals of the New York Academy of Sciences, vol. 1065, New York, 2005.
[4] S. Brandt, Data analysis. Statistical and computational methods for scientists and engineers. Springer Verlag, New York, 1999.
[5] J.M. Cornil, P. Testud, An introduction to Maple V. Springer Verlag, Berlin-Heidelberg, 2001.

Related Topics

Related Publications

Related Content

Site Footer


This website is provided by John Wiley & Sons Limited, The Atrium, Southern Gate, Chichester, West Sussex PO19 8SQ (Company No: 00641132, VAT No: 376766987)

Published features on are checked for statistical accuracy by a panel from the European Network for Business and Industrial Statistics (ENBIS)   to whom Wiley and express their gratitude. This panel are: Ron Kenett, David Steinberg, Shirley Coleman, Irena Ograjenšek, Fabrizio Ruggeri, Rainer Göb, Philippe Castagliola, Xavier Tort-Martorell, Bart De Ketelaere, Antonio Pievatolo, Martina Vandebroek, Lance Mitchell, Gilbert Saporta, Helmut Waldl and Stelios Psarakis.