Monte Carlo method for fractional-order differentiation extended to higher orders

In this work the Monte Carlo method, introduced recently by the authors for orders of differentiation between zero and one, is further extended to differentiation of orders higher than one. Two approaches have been developed on this way. The first approach is based on interpreting the coefficients of the Grünwald–Letnikov fractional differences as so called signed probabilities, which in the case of orders higher than one can be negative or positive. We demonstrate how this situation can be processed and used for computations. The second approach uses the Monte Carlo method for orders between zero and one and the semi-group property of fractional-order differences. Both methods have been implemented in MATLAB and illustrated by several examples of fractional-order differentiation of several functions that typically appear in applications. Computational results of both methods were in mutual agreement and conform with the exact fractional-order derivatives of the functions used in the examples.


Introduction
This paper is a continuation of our paper [8], where the Monte Carlo method for fractional differentiation was introduced and used for the approximation and evaluation of the Grünwald-Letnikov fractional derivative of order 0 < α < 1. The goal of the present work is the extension of the Monte Carlo method for fractional-order differentiation to higher orders.
Let us recall that the Grünwald-Letnikov fractional derivative is a non-local operator on L 1 (R) defined as [10] For f (t) = 0 for t < 0, the fractional-order difference (1.2) can be used for numerical evaluation of the Grünwald-Letnikov fractional order derivative, the Riemann-Liouville fractional derivative, and the Caputo fractional derivative when it is equivalent to the Riemann-Liouville one, [10].

A bit of history: signed probabilities and fractions of a coin
Since for the coefficients w k = γ (α, k), k = 1, 2, . . ., we have [8,9], ∞ k=1 (−γ (α, k)) = 1, (2.1) we can look at the coefficients w k as probabilities. However, not always all coefficients w k = γ (α, k) are positive, and in such situations we can consider them as signed probabilities, that can be positive or negative, while preserving the property (2.1). The notion of negative probability goes back to 1932 to Wigner's remark [16] that some considered expressions "…cannot be really interpreted as the simultaneous probability for coordinates and momenta, as is clear from the fact, that it may take negative values. But of course this must not hinder the use of it in calculations as an auxiliary function which obeys many relations we would expect from such a probability." In 1942, Dirac [4] who used negative energy and negative probabilities in quantum mechanics, emphasized that negative energy states occur with negative probability and positive energy states occur with positive probability and that it is possible to develop a theory would allow its application "to a hypothetical world in which the initial probability of certain states is negative, but transition probabilities calculated for this hypothetical world are found to be always positive", and concluded that "Negative energies and probabilities should not be considered as nonsense. They are well-defined concepts mathematically, like a negative sum of money, since the equations which express the important properties of energies and probabilities can still be used when they are negative. Thus negative energies and probabilities should be considered simply as things which do not appear in experimental results." It appears that Dirac's paper directly motivated Bartlett in 1945 to introduce signed probabilities [1], and one of the important consequences of allowing negative probabilities was that "…a negative probability implies automatically a complementary probability greater than unity…" This idea was also elaborated by Feynman [6] in 1987. Similarly to Bartlett, he mentioned that using negative probabilities leads to simultaneously allowing probabilities greater than one, but "probability greater than unity presents no problem different from that of negative probabilities", because those are used in intermediate calculations and the final result is always expressed in positive numbers between 0 and 1.
This development of the idea of allowing probabilities to be negative and positive naturally leads to the algebraic theory of probability, see [14]. As quoted above, negative probabilities can be used in the intermediate computations, while the final result must be physical (which means -positive, measurable). This is similar to using, say, negative resistors or capacitors in electrical circuits: there are no passive elements having negative resistance, but it is possible to create circuits exhibiting negative resistances, negative capacitances, and negative inductances. This was pointed out by Bode in his classical book [2, Chapter IX]. Such circuits are called negative-impedance converters ( [5,12]) and are based on using operational amplifiers. For example, the negative resistance of −10 kΩ means that if such an element is connected in series with a classical 20 kΩ resistor, then the resistance of the resulting connection is 10 kΩ.
A beautiful example of the interpretation of negative probabilities was developed by Székely [15], who introduced half-coins as random variables that take the values n with signed probabilities -positive for odd values of n and negative for even values. A fair coin is a random variable X that take the values 0 and 1 with probability 1/2. Its probability generating function is E z X = (1 + z)/2. The addition of independent random variables corresponds to multiplication of their probability generating functions. Therefore, the probability generating function of the sum of two fair coins is equal to ((1 + z)/2) 2 , and it is natural to define a half-coin via its probability generating function where for μ = 1/2 the coefficients p n = 1/2 n have alternating signs. When "flipping" two half-coins, the sum of the outcomes is 0 or 1 with probability 0.5, like in the case of flipping a normal coin.
Taking μ = 1/3 yields third-coins and "flipping" three of them makes a normal coin; μ = 1/4 defines a quarter-coin, etc. Non-integer values of μ produce μ-coins, and to get a normal coin it is necessary to "flip" 1/μ of such μ-coins, which depending on μ can be a finite or even infinite number of such "fractions of coins".
Extending the Monte Carlo method introduced in our paper [8] for the case of 0 < α < 1 to orders higher than one leads to working with signed probability distributions, because in such case the coefficients w k = γ (α, k) in the fractional-order difference (1.2) have different signs, and some of them have values outside of the interval (0, 1), as mentioned by Bartlett and Feynman, while still satisfying the property (2.1).

Monte Carlo fractional differentiation using signed probabilities
Since we deal with signed probabilities, it is necessary to follow the ideas of Dirac, Bartlett, and Feynman, which means that for using the Monte Carlo method we have to transform the signed probabilities to the positive probabilities. As a result, the expression used for Monte Carlo simulations will consist of two parts -one of them might contains terms with different signs and is independent of the random variables used in the Monte Carlo draws (trials), and the other contains the terms of the same (positive) sign and is dependent of those random variables. Below we first illustrate this for the cases 1 < α < 2 and 2 < α < 3, and then present a framework for the general case.

The case of order 1 <˛< 2
We have the following: Let Y 1 , Y 2 , …, Y n , …be independent copies of the random variable Y , then by the strong law of large numbers with probability one for any fixed t and h, and hence, as N → ∞, then by the central limit theorem and Slutsky lemma, as N → ∞, we have where → D denotes convergence in distributions and N (0, 1) is the standard normal law, and v N is a sample variance of random variables This allows us to build asymptotic confidence intervals (see details in [8]). The above results can be used as the basis of the Monte Carlo method for numerical approximation and evaluation of the Grünwald-Letnikov fractional derivatives.
Indeed, we can replace the sample Y 1 , Y 2 , …, Y N by its Monte Carlo simulations. For the simulation of the random variable Y with the distribution (3.1), we define Then If U is a random variable uniformly distributed on [0, 1], then

Using the binomial series
we easily obtain and denoting Thus, putting p 2 = q 1 + q 2 and p k = q k = −w k for k ≥ 3, we have the following: Let Y 1 , Y 2 , …, Y n , …be independent copies of the random variable Y , then for for fixed t and h such that and by the strong law of large numbers we obtain with probability one as N → ∞. Thus under assumption (3.5) we have the following convergence with probability one: and this can be used for evaluation of the Grünwlad-Letnikov fractional derivative by the Monte Carlo method. The simulation of Y 1 , Y 2 , …, Y n , …is similar to the previous cases of α ∈ (0, 1) and α ∈ (1, 2), with the only difference that and if F k−1 < U < F k , then Y = k + 1 (k = 1, 2, …), where U is a random variable uniformly distributed on [0, 1].

The general case for any˛> 0
It turns out that the methodology presented above for α ∈ (1, 2) and α ∈ (2, 3) can be extended to the case of arbitrary α > 0 in terms of generalized (signed) probabilities. Let us consider a generalized random variableȲ ∈ 1, 2, . . . such that where we interpret q k as generalized probabilityP(Ȳ = k), which are allowed to be negative, but for which there exists an integer r ≥ 1 such that q 1 + q 2 ∈ (0, 1), Then we can write the following formal identity for the Borel function g(k), k = 1, 2, . . .: which can be finally written as Thus, denoting π j+1 = q 2 j−1 + q 2 j , j = 1, 2, . . . , r , π r +1 = q r + j , j = r + 1, r + 2, . . . , we define the generalized expectation (3.12) where the ordinary discrete random variable Y is defined as andḡ(k), k ∈K , be the above function. If Y 1 , Y 2 , …, Y n , …are independent copies of the discrete random variable Y such that Eḡ(Y ) < ∞, then by the strong law of large numbers with probability one, and hence for the generalized random variableȲ with probability one but π 2 = q 1 + q 2 ∈ (0, 1), as well as for the rest q j = −w j ∈ (0, 1), j = 3, 4, . . ., and g(k) = f (t − kh), k = 0, 1, 2, . . ., for a given f and fixed r and h.

Monte Carlo method for fractional differentiation using the semi-group property of fractional differences
Another approach can be based on using the semigroup property [3,7] of the fractionalorder finite difference operator defined in (1.2). In fact, in this way we can simply follow our paper [8], where the function f (t) should be replaced by an integer-order fractional difference. Using the semi-group property of fractional-order differences, we can write Hence [8], Let Z ∈ {1, 2, . . .} be a discrete random variable with Then from (4.1) we have If Z 1 , Z 2 , …, Z m , …are independent copies of the random variable Z , then, as N → ∞, by the strong law of large numbers with probability one for α ∈ (n, n + 1), assuming that the last expectation exists, and hence with probability one as N → ∞ we have convergence to A α h f (t) defined by (1.2): Moreover, if for fixed f , t, and h the following inequality holds then by the central limit theorem and Slutsky lemma, as N → ∞, we have where → D denotes convergence in distributions and N (0, 1) is the standard normal law, and v N is a sample variance of random variables This allows us to build asymptotic confidence intervals (see details in [8]). The above results can be used as the basis of the Monte Carlo method for numerical approximation and evaluation of the Grünwald-Letnikov fractional derivatives. Indeed, we can replace the sample Z 1 , Z 2 , …, Z N by its Monte Carlo simulations.
For the simulation of the random variable Z with the distribution (4.3), we define Then If U is a random variable uniformly distributed on [0, 1], then and to generate Z ∈ 1, 2, . . .,

Examples
Both proposed methods have been implemented in MATLAB [13]. This allows their mutual comparison, as well as useful visualizations and numerical experiments with the functions that are frequently used in the applications of the fractional calculus and fractional-order differential equations. The Mittag-Leffler function [10] E α,β (z) = ∞ n=0 z n Γ (αn + β) , z ∈ C, α,β > 0, that appears in some of the provided examples, is computed using [11]. In all examples the considered interval is sufficiently large, namely t ∈ [0, 10]. The exact fractional derivatives are plotted using solid lines, the results of the proposed Monte Carlo method are shown by bold points, the results of K individual trials (draws) are shown by vertically oriented small points (in all examples, K = 200), and the confidence intervals are shown by short horizontal lines above and below the bold points.

Example 1. The power function
The particular case of ν = 0 is the Heaviside unit-step function, and its derivatives of orders α = 1.7 and α = 2.6 are shown in Fig. 1.
The derivatives of the power function for ν = 1.3 and orders α = 1.7 and α = 2.6 are shown in Fig. 2.

Example 2. The exponential function
Since the difference of the exponential function and the first terms of its power series can be expressed in terms of the Mittag-Leffler function, we can easily obtain the explicit expression for the corresponding fractional derivative [10]: The derivative of order α = 1.7 of the function y(t) = e λt − 1 − λt for λ = −1.4 and λ = −0.4 are shown in Fig. 3.

Example 3. The sine function
Similarly to the previous example, the fractional-order derivatives of the sine function can be obtained using its representation in terms of the Mittag-Leffler function: The results of evaluation of the derivatives of order 1.7 of the sine function using both methods are shown in Fig. 4 (method P is based on the signed probabilities approach, method S is based on using the semi-group property of fractional differences). The computed values of the fractional-order derivative are in mutual agreement and conform with the exact fractional-order derivatives. The confidence intervals are also of similar size, while the variance is much smaller in the case of the method based on the semi-group property of fractional-order differences; however, the computational cost is much higher. This observation holds for all other examples.

Example 4. The Mittag-Leffler function
The product of the power function and the Mittag-Leffler function appears frequently in solutions of fractional differential equations with Riemann-Liouville fractional derivatives, so it is also a suitable example: The results of computations for such function are shown in Fig. 5 for α = 1.7, μ = 1.5, β = 3.5, and λ = −1.4, resp. λ = −0.4.

Conclusions
Our extension of the Monte Carlo approach to fractional differentiation of orders higher than one led to working with signed probabilities that are not necessarily positive. We have demonstrated how this can be processed and used for computations using the Monte Carlo method.
The results of computations by the method based on signed probabilities are compared with the results obtained by the method that uses the semi-group property of fractional-order differences. Both methods produce practically the same values of fractional derivatives, which are in agreement with the values computed using the explicit formulas. The method based on the semi-group property is, by its construction, significantly less efficient as it requires more computations; at the same time it is characterized by smaller variance of the outputs of trials at the points of evaluations. The method based on signed probabilities is, on the contrary, faster and is characterized by larger values of variance of the trials at the points of evaluations. In both cases the confidence intervals are sufficiently small for practical purposes. The presented method can be further enhanced using standard approaches for improving the classical Monte Carlo method, such as reduction of variance, importance sampling, stratified sampling, control variates or antithetic sampling.

Conflict of interest The authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.