The calculation of the probability density and distribution function of a strictly stable law in the vicinity of zero

The problem of calculating the probability density and distribution function of a strictly stable law is considered at $x\to0$. The expansions of these values into power series were obtained to solve this problem. It was shown that in the case $\alpha<1$ the obtained series were asymptotic at $x\to0$, in the case $\alpha>1$ they were convergent and in the case $\alpha=1$ in the domain $|x|<1$ these series converged to an asymmetric Cauchy distribution. It has been shown that at $x\to0$ the obtained expansions can be successfully used to calculate the probability density and distribution function of strictly stable laws.


Introduction
The major inconvenience of using stable laws is the absence of expressions for probability density and distribution function in terms of elementary functions. There are only five cases known when the density is expressed in terms of elementary functions: the Lévy distribution (α = 1/2, θ = 1) symmetric Lévy distribution (α = 1/2, θ = −1), Cauchy distribution (α = 1, θ = 0), The Gaussian distribution (α = 2, θ = 0) and asymmetric Cauchy distribution (α = 1, −1 θ 1) (see formulas (5) and (8)). Here α is characteristic exponent of a stable law, θ -is a parameter of asymmetry. The latter distribution first came out in the book by V.M. Zolotarev [1] (see formula (2.3.5a)) and later was examined in the works [2,3]. Different representations for stable laws are required to calculate the probability density or distribution function in other cases.
The paper [4] shows that if values of the characteristic exponent α and the asymmetry parameter β are limited by values of rational numbers (α = P/Q, β = U/V , where P, Q, V are positive integers), then in this case it is possible to express the probability density of a strictly stable law in terms of special functions. The papers [4][5][6][7][8][9][10] are devoted to obtaining such representations. The limitation of this approach lies in the fact that it is possible to obtain an expression for the probability density only for rational values of the parameters α and β, and only for strictly stable laws. The application of the Fast Fourier Transform algorithm is another method of calculating density. This approach has been examined in the papers [11,12]. However, this method gives an opportunity to calculate the probability density on a grid of equidistant points. In the paper linear interpolation must be used to calculate the density at intermediate points or at irregularly spaced points.
The use of integral representations is the main method for calculating the probability density and the distribution function of stable laws. This approach is based on the inversion formula (2). There are two possible ways of inverting the characteristic function. The first way is to directly calculate the integral in (2). As a result, the probability density is expressed in terms of the integral of the oscillating function [13,14]. However, since the integrand is an oscillating function, this leads to difficulties in numerical integration in the cases α < 0.75, β = 0 and 0 < |α − 1| < 0.001 and in the case of large values x [13]. Modernization of the standard quadrature method of numerical integration makes it possible to reduce the lower boundary of the parameter α from the value 0.75 to the value 0.5 [14]. It is proposed to use the representation of the density in the form of a power series to calculate the density for large values of x.
The second way of obtaining integral representations is the application of the stationary phase method when calculating the integral in (2) (see [1,3,15,16]). The advantage of this method of inverting the 1 characteristic function is that the resulting integral representation is expressed in terms of a definite integral of a monotonic function. Such integral representations were obtained for stable laws with different parameterizations of the characteristic function: for parameterization "B" in the works [1,15], for parameterization "M" in the paper [16], for parameterization "C" in the paper [3]. (Here, the notation of various parameterizations of the characteristic function is given in accordance with the designations introduced in the book by V.M. Zolotarev [1].) These integral representations are more convenient from a practical point of view and allow calculating the density in a wide range of parameter values α, β and coordinates x. The integral representation obtained in the work [16] served as a foundation for developing several software products [17][18][19][20][21].
From a theoretical point of view, these integral representations are valid for all values of x. However, in practice, it is not possible to calculate the probability density and distribution function for all values of x. The reason for this lies in the behavior of the integrand. The integrand has the form of a very sharp peak with small and large values of x. As a result, numerical integration algorithms cannot correctly calculate the integral in this range of x. To settle this issue in the papers [16,18,19] it is proposed to use various numerical methods to increase the accuracy of calculations. However, all proposed approaches increase the accuracy of the calculation, but do not completely eliminate the problem. To calculate the probability density and distribution function in this range of values of x it is expedient to use other representations for stable laws which do not have any specific features in the indicated areas. The approach used in the papers [12,14] seems to be the most suitable which consists in applying expansions in a power series for probability density and distribution function with x → 0 and x → ∞.
Such expansions are well known and are obtained, as a rule, for parametrization "B". Depending on the value of the parameter α the obtained power series is either convergent or asymptotic. The expansion of the probability density of a stable law into a convergent series in the case x → ∞ and 0 < α < 1, was firstly mentioned in the paper [22]. Later, in the paper [23] a generalization of this density expansion was given for x → ∞ in the case 1 < α < 2. In this range of values of the parameter α this series turns out to be asymptotic. In the same paper, the expansion of the density in a series in the vicinity of the point x → 0 was obtained for the case 0 < α < 2. The resulting power series is asymptotic in the case 0 < α < 1, and convergent in the case 1 < α < 2. Expansions for α > 1 in the cases x → 0 and x → ∞ were also obtained in the work [5] as a result of expansion into a power series of the probability density, expressed in terms of the Fox function. The same expansions were given in the books [24] (see Chapter 17, §7) and [1] (see §2.4 and §2.5). Expansions of the density of a stable law in a power series for the characteristic function in parameterization "M" were obtained in the paper [14]. An interesting result was obtained in the paper [25]. In this paper, expansions in power series were obtained for the probability density of a symmetric stable law at x → 0 and x → ∞ for the cases 0 < α < 1 and 1 < α < 2. A distinctive property of this expansion is that these power series for all α are convergent.
The purpose of this work is to obtain power series expansions of the probability density and distribution function of a strictly stable law with the characteristic function where α ∈ (0, 2], |θ| min(1, 2/α − 1), λ > 0. This parameterization of the characteristic function, according to the book [1], is called parameterization "C". Obtaining such expansions turns out to be necessary in connection with the problem of calculating the probability density and distribution function of stable and fractionally stable laws. In fact, in the article [3] integral representations were obtained for the probability density and distribution function of a strictly stable law with the characteristic function (1). Since these integral representations were obtained using the stationary phase method, then with small and large values of the coordinate x the integrand has the form of a very sharp peak. This causes difficulties for numerical integration algorithms and leads to incorrect integration results. Therefore, to calculate the probability density and distribution function in these coordinate regions, it is expedient to use representations in the form of a power series for the corresponding quantities. This work is devoted to obtaining such expansions.
The solution to this problem will turn out to be useful not only when calculating the density of strictly stable laws but also in the task of calculation the density and distribution function of a fractional-stable law [2,26,27]. These distributions are expressed in terms of the Mellin convolution of two strictly stable laws. Correct calculation of the probability density will make it possible to use an algorithm for statistical estimation of the parameters of these laws based on the maximum likelihood method. Such an algorithm for estimating parameters will give an opportunity to correctly describe various experimental data. It is known that the distribution of gene expression is described by laws with a power-law decrease in density [28][29][30]. Since the stable and fractionally stable densities decrease according to the power law x −α−1 at x → ∞, then these classes of distributions were used to describe the distribution of gene expression. In the works [31,32] fractional stable distributions were used to describe the expression of genes obtained using microarray technology. In the work [33] these distributions were used to describe the results obtained using the Next Generation Sequence technology. To describe these experimental data, it is necessary to have algorithms for statistical estimation of parameters, the most effective of which is the maximum likelihood method. To construct such an algorithm, it is necessary to be able to correctly calculate the density of a strictly stable law for any values of x.

Preliminary remarks
The major purpose is to obtain the expansion of the density and distribution function of a strictly stable law in a power series in the vicinity of the point x = 0. The paper deals with strictly stable laws with the characteristic function (1). Without loss of generality, we will assume that the scale parameter λ = 1. Strictly stable laws with the parameter λ = 1 are commonly called standard strictly stable laws. Designation abbreviations are accepted for standard strictly stable laws. The characteristic function will be designated byg(t, α, θ, 1) ≡ĝ(t, α, θ), the probability density distribution will be designated by g(x, α, θ, 1) ≡ g(x, α, θ), the distribution function will be designated by G(x, α, θ, 1) ≡ G(x, α, θ).
To perform the inverse Fourier transform and obtain the probability density distribution, the following lemma is useful, which defines the inversion formula Lemma 1. The probability density function g(x, α, θ) for any admissible set of parameters (α, θ) and any x can be obtained using the inversion formulas The proof of this lemma can be found in the paper [3]. To obtain the probability density, there is no fundamental difference which of the formulas to use on the right side (2). The result will differ only in the sign of the parameter θ. Without loss of generality, in this paper we will use the first formula (2). Such a choice results from the fact that in the works [1,3,34] this formula was used to invert the characteristic function. This will give us an opportunity to compare the results obtained below with the results of the mentioned papers without any additional transformations.
In the article [3] the inverse Fourier transform of the characteristic function (1) was performed and expressions for the probability density and distribution function of a strictly stable law were obtained. In the case α = 1 and x = 0 for any admissible θ the following integral representation is true for the probability density where θ * = θ sign(x) and If α = 1, then for any admissible −1 θ 1 the probability density has the form g(x, 1, θ) = cos(πθ/2) π (x 2 − 2x sin(πθ/2) + 1) .
If we use Euler's formula cos β + i sin β = e iβ , then this integral can be represented in the form 3. Representation of the probability density in the form of a power series We obtain the expansion of the probability density g(x, α, θ) in a series at x → 0. The following theorem is valid Theorem 1. In the case x → 0 for any admissible set of parameters (α, θ) except for the values α = 1, θ = ±1 for the probability density g(x, α, θ) the following representation in the form of a series is valid where Proof We will perform the inverse Fourier transform of the characteristic function (1). To do this we make use of the first relation in (2). We have Since the considered case x → 0, then we expand exp{itx} in a series in the vicinity of the point x = 0. As a result, we get where the N -th partial sum g 0 N (x, α, θ) and the remainder R 0 N (x, α, θ) of a series have the form Here R N (itx) = (itx) N N ! e itxζ , (0 < ζ < 1) is the remainder in the Lagrange form. We consider the N -th partial sum g 0 N (x, α, θ). To calculate the integral in (15), we will change in some places the order of summation and integration and we will substitute the integration variable t α = τ . As a result, we obtain Next, we examine the range of valid values of the argument π 2 αθ. The range of admissible values of the parameter θ is determined by the inequality |θ| min(1, 2/α − 1). Hence, if 0 < α 1, then −1 θ 1, if 1 < α 2, then −(2/α − 1) θ 2/α − 1. Thus, π 2 αθ ∈ π 2 α, π 2 α ∈ − π 2 , π 2 , if 0 < α 1.
Combining (18) and (19), we obtain As we can see, extreme values of this range are reached in the case α = 1 and θ = ±1.
Now we consider the remainder R 0 N (x, α, θ). From the expression (16) we get It is not possible to calculate this integral since the exact value of the quantity ζ is not known. It is only known that 0 < ζ < 1. However, one can obtain an estimate for this integral. We have To obtain the third inequality, it was taken into consideration that | exp{itxζ}| 1. Next, the integration variable was substituted t α = τ . To calculate the resulting integral, the formula (10) was used. The obtained expression completely proves the theorem.
We need to make one small remark. When proving the theorem, it was pointed out that it was necessary to exclude the case π 2 αθ = ± π 2 from consideration, which corresponds to the values of parameters α = 1, θ ± 1. As part of the proof of the theorem, this was done so that the range of admissible values of the argument π 2 αθ of the integral in (17), should coincide with the range of admissible values of the argument β, included in the integral (10). However, the exception of the case β = ± π 2 from the integral (10) is related with the fact that in these points the integral (10) will diverge (for details see [36]). Therefore, it should be assumed that α = 1 and θ = ±1, then the integral in (17) will diverge. This in its turn leads to a degenerate probability density at that point. As a result, we arrive at the well-known fact that the probability density with the characteristic function (1) is degenerate in the points α = 1, θ ± 1.
As noted in the Introduction, depending on the value of the parameter α, the expansion of the probability density of the stable law in a power series turns out to be either convergent or divergent. Absolutely the same situation occurs in the considered case. The answer to the question under what values of α the expansion (11) is convergent, and for which it is divergent we formulate as a corollary Corollary 1. In the case α < 1 the series (12) is divergent at N → ∞. In this case for the density g(x, α, θ) for any admissible θ the asymptotic expansion is valid.
In the case α = 1 the series (12) converge for any x, satisfying the condition |x| < 1. In this case for the density g(x, 1, θ) for any admissible θ = ±1 it is possible to represent in the form of an infinite series In the case α > 1 the series (12) converge for any x. In this case for the density g(x, α, θ) for any admissible θ it possible to represent in the form of an infinite series Proof We examine the convergence of the series (12). As we can see, this series is sign-alternating. Consequently We apply the Cauchy criterion in the limiting form to the resulting series.
Here the Stirling's formula was used From the result obtained we can see that for the values α < 1 the series (12) diverges for any x, with the value α = 1 the series (12) converges for any values of x, satisfying the condition |x| < 1, and in the case α > 1 the series (12) converges for any x.
We consider the case α < 1. In this case the series (12) diverges at N → ∞. However, from the expression (13) it follows that for some fixed N Thus, with every N we have As a result, we have obtained the definition of an asymptotic series. Consequently, We consider the case α > 1. From the expression (11) it follows that We will set some arbitrary x and consider the limit of the right-hand side of this inequality under the condition N → ∞ . We have Thus, the right-hand side (24) represents an element of an infinitesimal sequence. This means that the sequence g 0 N (x, α, θ) at N → ∞ converges to the density g(x, α, θ). Therefore, in the case α > 1 for any fixed x for the density g(x, α, θ) the representation in the form of an infinite series is valid Now we consider the case α = 1. From the expression (11) it directly follows We fix some arbitrary x and find the redistribution under the condition N → ∞. As a result, we obtain Thus, the right side of the previous expression at |x| < 1 is an element of an infinitesimal sequence. Therefore, the sequence g 0 N (x, 1, θ) converges to the density g(x, 1, θ) at N → ∞ and |x| < 1. Now substituting the value α = 1, in the series (12) we obtain (22).
The proved corollary shows that in the case α = 1 in the interval −1 < x < 1 the series (22) converges to the density g(x, 1, θ). It is important to show that in this case the series (22) converges to the density (5). We formulate this result in the form Remark 1. In the case α = 1 for any −1 < θ < 1 in the region −1 < x < 1 the series (22) converges to the density (5).
Proof To prove this remark we will consider the density (5) and show that the expansion of this density in a Taylor series in the vicinity of the point x = 0 has the form (22). For the convenience of further presentation, we use the reduction formulas cos π 2 θ = sin π 2 (1 − θ) , sin π 2 θ = cos π 2 (1 − θ) and represent the density (5) in the form . Now we expand the density g(x, 1, θ) in a Taylor series in the vicinity of the point x = 0. Since this density is an infinitely differentiable function, we have We will draw attention that the function g(x, 1, θ) is a complex function. We will introduce the designa- In view of the introduced designations, the density (5) takes the form To calculate the n-th derivative we use the Bruno formula where Y n (f g 1 , f g 2 , . . . , f g n ) are the Bell polynomials (see [37]) Here f m ≡ f m , m = k 1 + k 2 + · · · + k n , the sum is taken over all solutions to the equation and Taking into account (26), we get For coefficients g j we have This shows that in the expression (28) the sum contains the summands that satisfy the equation Indeed, in the expression (28) the summation is done over all solutions to the equation (29). In case, if the solution k j = 0, j = 3, 4, . . . , n, then the corresponding term in the sum will be equal to zero, since g j = 0, j = 3, 4, . . . , n. If k j = 0, j = 3, 4, . . . , n, then the multiplier (g j /j!) k j = 1, since 0 0 = 1. Consequently, in the expression (28) there are the summands that satisfy the solution to the equation (32). This significantly simplifies the summation. It follows from the equation (32) that k 1 = n − 2k 2 . Taking into consideration that k 1 0 and k 2 0, we obtain k 2 = 0, 1, 2, . . . , n 2 , where [A] means the integer part of the number A. It gives an opportunity to introduce directly the summation index in the sum (28). In view of the foregoing, the formula (28) takes the form where the relation k 1 = n − 2k 2 is used and the summation index k ≡ k 2 is introduced. Now substituting this relation in (27) and using (30) and (31), we obtain Now we calculate the value of this derivative in the point x = 0. It is easy to see that where it is taken into account that (−1) 2n−3k = (−1) k . Next, we use the general formula for sin(nϕ) (see, for example, [38]) Using this formula in (34), we obtain Using this expression now in (25), we get Thus, the expansion of the density (5) into an infinite Taylor series in the vicinity of the point x = 0 agrees exactly with the series (22). It completely proves the corollary.
Theorem 1 gives an opportunity to determine the range of values x within which the absolute error of the density calculation using the series (12) and for some fixed N will not exceed the predetermined value ε. This turns out to be very convenient when calculating the probability density. Indeed, from the relations (11) and (13) we obtain If, for a given value of N we set the absolute value of the error |g(x, α, θ) − g 0 N (x, α, θ)| = ε, then it becomes possible to introduce the threshold coordinate This value shows that for coordinates |x| |x N ε | the absolute value of the density calculation error using the series (12) will not exceed ε, i.e.  Thus, to calculate the probability density, we can use the N -th partial sum (12) in the range of coordinates In this case, the magnitude of the absolute error at fixed N will not exceed the chosen value ε. Figures 1a and 2a show the results of calculating the probability density g(x, α, θ) using the series (12). In these figures, the solid curve corresponds to the exact density values g(x, α, θ), calculated with the help of (3), the dashed-dotted curve corresponds to the density calculation results using the series (12) for the specified values of N . Figures 1b and 2b show the results of the calculation of the absolute error |g(x, α, θ) − g 0 N (x, α, θ)|. In these figures, the solid curves correspond to the exact value of the absolute error |g(x, α, θ) − g 0 N (x, α, θ)|, where g(x, α, θ) -the exact density value calculated when using (3), g 0 N (x, α, θ) -the series (12), the dashed curve corresponds to the estimate of the remainder term (13). The calculation results are given for the specified values of N in the figures. In all these figures, the circles show the position of the threshold coordinate x N ε for the selected level of accuracy ε and each number of summands N . It is clear from Figure 1b and 2b that in the region |x| x N ε the absolute magnitude of the error does not exceed the specified level of accuracy ε for all N . This means that at |x| x N ε the expansion (12) can be used to calculate the density. Corollary 1 shows that in the case α < 1 the series (12) is divergent at N → ∞, and in the case α > 1 this series converges. The cause of this behavior lies in the ratio Γ((n + 1)/α)/Γ(n + 1), which is present in this series. At α < 1 this ratio turns out to be more than unity and, as n increases, this ratio only rises. Therefore, to achieve the specified calculation accuracy ε one has to decrease the value of x. It is clearly seen from the behavior of the threshold coordinates x N ε . Figures 1a and 1b show that at first the addition of summands in the expansion (12) leads to an increase in the range of x for which the inequality |g(x, α, θ − g 0 N (x, α, θ))| ε is satisfied. The fact that x 3 ε < x 10 ε < x 30 ε testifies to it. However, further addition of summands leads to an increase in the ratio Γ((n + 1)/α)/Γ(n + 1) and, thus, to an increase in the absolute calculation error. Therefore, to achieve the specified level of accuracy, it is necessary to decrease the value of the coordinate x. This causes the threshold coordinate x N ε to start decreasing and we see that x 100 ε > x 300 ε . In the case α > 1 the situation changes. In this case the ratio Γ((n + 1)/α)/Γ(n + 1) < 1 and as n increases this ratio only decreases. Consequently, the increase in the number of summands in (12) increases the accuracy of the density calculation. This leads to the fact that the range of values of x, for which the inequality |g(x, α, θ) − g 0 N (x, α, θ)| ε is met increases with the addition of the number of summands N in the sum (12). This is clearly seen from the location of the threshold coordinates x N ε , shown in Figures 2a and 2b. We can see from the figures that x 3 ε < x 10 ε < x 30 ε < x 100 ε < x 300 ε . Thus, in the case α > 1 the series (12) is convergent for all x at N → ∞.
The results of calculations for the case α = 1 are given in Fig. 3. This figure shows, ε . Thus, an increase in the number of summands N in the sum (12) leads to an increase in the interval x, within which the inequality |g(x, α, θ) − g 0 N (x, α, θ)| ε is satisfied. Here, the formula (5) was used to calculate the density. It is also seen from the figure that for all presented N at x → 1 the partial sums of the series diverge. This is in full agreement with the statement of corollary 1 which states that in the case α = 1 the series (12) converges at |x| < 1.

Representation of the distribution function in the form of a power series
Now we will try to obtain the representation of the distribution function in the case x → 0 in the form of a power series. We will formulate the result obtained as a theorem Theorem 2. In the case x → 0 for any admissible set of parameters (α, θ), except for the values α = 1, θ = ±1 for the distribution function G(x, α, θ) a representation in the form of a power series is valid. x n+1 (n + 1)! Γ n + 1 α sin π 2 (n + 1)(1 − θ) , Proof From the definition of the distribution function, it follows Here G(0, α, θ) is the value of the distribution function in the point x = 0 and is determined by the formula (9). Using the expansion (14) for the density g(x, α, θ) we obtain where Here g 0 N (x, α, θ) and R 0 N (x, α, θ) are determined by the expressions (12) and (16) respectively.
To calculate the partial sum G 0 N (x, α, θ) we will make use of the results of theorem 1. It was obtained in this theorem that the partial sum g 0 N (x, α, θ) has the form (12). Substituting the expression (12) in (41) and changing the order of integration and summation, we get Now we obtain the expression for the remainder R 0 N (x, α, θ). Substituting the expression (16) in (42) and changing the order of integration, we obtain where 0 < ζ < 1. This integral cannot be calculated directly, since the exact value of ζ is not known. It is only known that ζ ∈ (0, 1). However, one can obtain an estimate of this integral.
To obtain an estimate for the integral, we use the inequality Here, to calculate the outer integral, the integration variable t α = τ was first substituted, and then the formula (10) was used. It should be noted that the case α = 1, θ = ±1 must be excluded from consideration. Indeed, for such parameter values, the argument π 2 αθ = ± π 2 and integral (10) will diverge. Now substituting the expressions (43) and (44) in (43) we get the statement of the theorem.
The proved theorem shows that in the vicinity of the point x = 0 the expansion (37) is valid for the distribution function of a strictly stable law with the characteristic function (1). However, as in the case of the probability density, the obtained power series diverges for all x at α < 1, and in the case α > 1 is convergent for all x. In the case α = 1 this series converges at |x| < 1 and diverges at |x| 1. In this regard, for the values α < 1 the representation (37) is asymptotic, and for the values α > 1 the expansion G(x, α, θ) can be represented in the form of an infinite power series. We formulate this result as a corollary. Corollary 2. In the case α < 1 the series (38) diverges for all x at N → ∞. In this case the asymptotic expansion is valid for the distribution function G(x, α, θ) for any admissible θ x n+1 (n + 1)! Γ n + 1 α sin π 2 (n + 1)(1 − θ) , x → 0.
In the case α = 1 the series (37) converges at |x| < 1. In this case the distribution function G(x, 1, θ) for any θ = ±1 can be represented as an infinite series In the case α > 1 the series (38) at N → ∞ converges for any x. In this case the representation in the form of an infinite power series is true for the distribution function G(x, α, θ) for any admissible θ Proof We examine the convergence of the series (38). It is clear that this series is sign-alternating. Consequently We apply the Cauchy criterion in the limiting form to the obtained series. Using Stirling's formula (23) and taking into consideration that n + 2 ≈ n + 1 at n → ∞, we get This shows that in the case α < 1 the series (38) diverges for all x, in the case α > 1 the series converges for all x, and in the case α = 1 the series (38) converges if |x| < 1. Now we consider the case α < 1. In this case at N → ∞ the series (38) diverges. However, it follows from the expression (39) that for some fixed N Consequently, for each N we have Thus, we have obtained the definition of an asymptotic series. Consequently, Now we consider the case α > 1. In this case the series (38) is convergent. It follows from the expressions (37) and (39) that We will find the limit at N → ∞ of the right-hand side of this inequality. Using Stirling's formula (23) and taking into account that N + 2 ≈ N + 1 at N → ∞, we obtain Thus, in the two cases α > 1 and α = 1, |x| < 1 the right side of the inequality (46) is an element of an infinitesimal sequence. In its turn, this means that in the above two cases, for any fixed x, the sequences (1 − θ)/2 + G 0 N (x, α, θ) converge to the distribution function G(x, α, θ). Therefore, in the considered case α > 1 for any fixed x the distribution function can be represented as an infinite series.
Now we consider the case α = 1. As shown above, in this case, when the condition |x| < 1 is met, the right side (46) is an element of an infinitesimal series. Therefore, for any fixed |x| < 1 the representation in the form of an infinite series is true for the distribution function Thus, the corollary has been proved completely.
As in the case of the probability density, the proved property shows that in the case α = 1 and |x| < 1 the series (45) converges to the distribution function G (x, 1, θ). It is possible to show that this series converges to the distribution function (8). We will formulate this result as a remark Remark 2. In the case α = 1 for any −1 < θ < 1 in the region −1 < x < 1 the series (45) converges to the distribution function (8).
Proof To prove this remark, we proceed in the same way as in the proof of remark 1. Let us show that the expansion of the distribution function (8) into a Taylor series in the vicinity of the point x = 0 has the form (45). We will use the reduction formulas cos π 2 θ = sin π 2 (1 − θ) and sin π 2 θ = cos π 2 (1 − θ) and will write the distribution function (8) in the form .
Note that the function arctan(x) is infinitely differentiable, therefore, expanding it into an infinite series, we obtain For the derivative of the order n we have Thus, the problem has been reduced to calculating the derivative n − 1 of the probability density g(x, 1, θ). However, this problem has been solved by us when proving remark 1. Using the formula (33) we get Substituting this expression in (48) and calculating the value of the obtained derivative in the point x = 0 and then using (35), we obtain where it was taken into account that (−1) 2n−2−3k = (−1) k . Substituting now the obtained expression for the n-th derivative in (47) and taking into consideration (9), we get Here, in the last equality, the summation index n = k + 1 was changed. Thus, the expansion of the distribution function (8) If now, for a specified fixed N we set the absolute magnitude of the error G(x, α, θ) − 1 2 (1 − θ) − G 0 N (x, α, θ) = ε, then it is possible to introduce the threshold coordinate This value shows that for all x satisfying the condition |x| x N ε , the absolute magnitude of the error in calculating the distribution function using the expansion (37) will not exceed the value ε: Figures 4a and 5a show the calculation results of the distribution function using the integral representation (6) (solid curves) and using the expansion (37) (dash-dotted curves) for parameter values α = 0.6, θ = 0.5 and α = 1.2, θ = 0.5 respectively. These figures contain the results of calculating the distribution function using the expansion (37) for the values N = 3, 10, 30, 100, 300. Figures 4b and 5b show the results of calculating the absolute error. In these figures the dashed curve corresponds to the estimate of the Here G(x, α, θ) is the exact value of the distribution function calculated using the representation (6), G 0 N (x, α, θ) is determined by (38). From Figures 4b and 5b it is clear that the condition (50) is met for all given values of N . In these figures the location of the threshold coordinate x N ε is marked with circles and the dotted line corresponds to the specified level of accuracy ε. We can see from the presented figures that for all values of x, satisfying the condition |x| x N ε , both the estimate of the remainder (39) (dashed lines), and the exact value of the absolute error (solid curves) are below the specified level of accuracy ε. This confirms the validity of the condition (50) and shows that the formula (49) can be used to estimate the boundary value of the coordinate in the expansion (37) at which the specified level of accuracy is achieved.
It should be noted that in the case α < 1 and α > 1 the threshold coordinate x N ε behaves differently as the number of summands N in the expansion (37) increases. In the case α < 1 (Fig. 4) an increase in N first increases the threshold coordinate x N ε (x 3 ε < x 10 ε ), but with the further increase in N the threshold coordinate x N ε decreases (x 30 ε > x 100 ε > x 300 ε ). The threshold coordinate x N ε behaves quite differently in the case α > 1. In this case with an increase in N the value of the threshold coordinate increases: Fig. 5). Such behavior of the threshold coordinate x N ε is due to the fact that in the case (α < 1) the series (38) is divergent, and in the case α > 1 this series converges (see corollary 2).
In the case α = 1 the threshold coordinate behaves in the same way as the case α > 1. With an increase in the number of summands of N in the expansion (37) the value of the threshold coordinate x N ε increases. We can see it from Fig. 6, which contains x 3 ε < x 10 ε < x 30 ε < x 100 ε < x 300 ε . However, unlike the previous case, lim N →∞ x N ε = 1. Indeed, in the case α = 1 the formula (49) takes the form Such behavior of the threshold coordinate x N ε is a consequence of proved corollary 2. Indeed, in the case α = 1 the series (38) and, therefore, the representation (37) converges in the region |x| < 1.
The results of calculating the absolute error in the case α = 1 are given in Fig. 6b. In this figure, the value of the threshold coordinate x N ε for different values N is shown with a circle. We can see from the presented results that for the values |x| x N ε both the estimate of the remainder (39) (dashed lines), and the exact value of the absolute error (solid curves) turn out to be less than specified accuracy level ε (dotted line). This demonstrates that the use of the formula (49) to estimate the values of the boundary coordinate leads to the validity of the condition (50).

5.
Calculation of the probability density and distribution function for small x We return to the problem of calculating the probability density of a strictly stable law. As mentioned in the Introduction, the main approach to the calculation of the probability density is to use the integral representation. For a strictly stable law with the characteristic function (1) such an integral representation is determined by the formula (3). This formula is valid for any x = 0 and any admissible values of parameters α and θ except for α = 1. However, in practice, it is not possible to calculate the integral in (3) numerically for all values of x. The reason for this lies in the behavior of the integrand. Figure. 7 shows the graph of the integrand in the formula (3) depending on the integration variable ϕ for different values of x. The graph of the function is plotted on a semi-logarithmic scale. We can see from this figure that as the value of x decreases, the integrand turns into a function with a very narrow and sharp peak. With a further decrease in x this peak becomes even narrower and higher. The same behavior of the integrand is also observed for large values of x. This leads to the fact that for very small and for very large values of x numerical integration algorithms cannot calculate the integral of this function. Figures 8 and 9 show the results of calculating the probability density using the integral representation (3) (solid curves). Fig. 8 shows the case α < 1, Fig. 9 shows the case α > 1. The Gauss-Kronrod algorithm was used to calculate the integral in the formula (3). It is clear from the presented calculations, for small values of x the numerical integration method used is cannot calculate the integral in (3). The critical value of the coordinate x cr at which the numerical integration algorithm used begins to produce an incorrect result for α = 0.3 is x cr ≈ 9 · 10 −10 , for α = 0.6 x cr ≈ 7 · 10 −7 , for α = 0.9 x cr ≈ 3 · 10 −4 (see Fig. 8). In the case α > 1 (see Fig. 9) for the value α = 1.1 x cr ≈ 3 · 10 −6 , for the value α = 1.4 x cr ≈ 10 −10 , and for the value α = 1.7 x cr ≈ 10 −11 . Consequently, at x < x cr it is necessary to use other methods to calculate the probability density. The same problem exists for integral representations of the density of stable laws in other parameterizations of the characteristic function (see [14,16,18,19]). To solve this problem in these works, the authors used various numerical methods, which make it possible to increase the accuracy of the calculation. However, these methods increase the accuracy of calculations, but do not solve the problem completely.
To calculate the density for the values |x| < x cr one should use other representations that do not have any singularities in this area. The most suitable option for this purpose is the power series representation obtained in Theorem 1 for the probability density. The estimate of the remainder (13) obtained in the same theorem made it possible to obtain the formula for the threshold coordinate x N ε (36) at which the given value of the absolute error ε is achieved for fixed N . This means that in the region −x N ε x x N ε the absolute error of the density calculation using the series (11) will not exceed the specified value ε. In Figures 8 and 9 the dash-dotted curves show the results of calculating the probability density using the series (12) for the specified values of α. The position of the threshold coordinate x N ε is shown with circles. The values x N ε are calculated for the absolute error ε = 10 −5 and N = 10. These figures show in the region x cr x x N ε the results of calculating the probability density using the integral representation (3) and using the series (12) coincide. For the values |x| x cr the numerical integration algorithm no longer allows obtaining the correct density value. At the same time, the calculation of the probability density using the series (12) does not cause any difficulties. It follows that for the values |x| x N ε it is expedient to use the series (12) to calculate the probability density. Thus, using theorem 1 and, in particular, the series(12) completely solves the problem of calculating the probability density at x → 0.
Similar problems arise when calculating the distribution function using the integral representation Solid curves -the formula (3), dash-dotted curves -the representation (11) for N = 10, circles -the position of the threshold coordinate x N ε (36). While calculating ε = 10 −5 was used . Calculations are given for the specified values of α and θ = 0.9 Fig. 9. Probability density for the case α > 1. The notation is the same as in Fig. 8 (6). The integrand in this integral representation also has some singularities at x → 0. In the general case, the integrand in (6) (see also (7)) behaves in the following way. In the point of the lower limit ϕ = −πθ/2 the integrand is equal to 1, in the point of the upper limit ϕ = π/2 the value of the integrand is equal to 0. As the variable ϕ increases from the value −πθ/2 up to the value π/2 the integrand decreases monotonically from 1 to 0. However, for very small values x the integrand in (7) decreases very sharply from 1 to 0 in a very narrow range ϕ. As a result, some numerical integration algorithms cannot recognize such a sharp decrease in the function and give an incorrect integration result. To exclude the possibility of incorrect results completely for small values x, it is expedient to use the expansion of the distribution function into a series obtained in theorem 2. The estimate of the remainder obtained in this theorem made it possible to obtain the formula (49) for the threshold coordinate. The value x N ε enables us to determine the range of x at which the inequality (50) is satisfied. In other words, in the range of values |x| x N ε the absolute error of calculating the distribution function using the expansion (37) will not exceed the value ε, where ε is given by in advance. Therefore, when calculating the distribution function in the range of values |x| x N ε t is expedient to use the expansion (37), and for |x| > x N ε the integral representation (6).

Conclusion
The major approach to the calculation of the probability density and the distribution function of stable laws is the use of integral representations. Theoretically, these representations are valid for all values of the coordinate x. However, it is not possible to calculate the density numerically for all x. Problems arise in the domain of very small and very large values of x. Therefore, it is expedient to use other methods for numerical calculations.
The paper considers the problem of calculating the probability density and distribution function in the case of x → 0 for a strictly stable law with a characteristic function (1). To solve this problem, expansions of the probability density and distribution function in a power series and estimates of the residual terms for each of the expansions were obtained. Estimates of the threshold coordinates x N ε were obtained for the expansion of the probability density and the distribution function, which are defined by the expressions (36) and (49), respectively. The threshold coordinate allows one to determine the domain of coordinates |x| x N ε within which the absolute computational error will not exceed the specified accuracy level ε. The performed calculations showed that the value of the critical coordinate x cr , at which the numerical integration algorithm used starts giving an incorrect result, is significantly less than the threshold coordinate x N ε (see Fig. 8 and 9). This fact shows that in the domain |x| x N ε it is possible to use theorems 1 and 2 to calculate the probability density and distribution function. The analysis of the obtained series made it possible to confirm both the known properties of these series and to establish new properties, as well as to improve the known estimates of the remainder terms. It was shown that in the case α < 1 the power series were divergent for any x at N → ∞. In this case these series are asymptotic at x → 0. In the case α > 1 the obtained series are convergent for all admissible x. In this case, representations in the form of infinite series are valid for the probability density and distribution function (see corollaries 1 and 2). These results are known and were previously obtained for the characteristic function in parameterization B in the works [23], [24] (see Chapter 17, §7), [1] (see §2.4 and §2.5), [34] (see §4.2, §4.3). It can be shown that the expansions from corollaries 1 and 2 in the cases α < 1 and α > 1 completely correspond to the expansions in the mentioned works. Examining the case α = 1 helped us establish that the expansions of the probability density and distribution function converged to the probability density (5) and the distribution function (8) in the domain |x| < 1 for any −1 < θ < 1 (see Remarks 1 and 2). The paper improves the estimates of the remainder terms in the expansions of the probability density and the distribution function defined by the formulas (13) and (39). The estimate of the remainder term obtained earlier (see [1], formula (2.5.2)) refers to the expansion of the probability density in parameterization B and in the case of α < 1 has the form We will take the relation θ = βK(α)/α into account, where K(α) = α − 1 + sign(1 − α), which relates the asymmetry parameter β in parameterization B to the asymmetry parameter θ in parameterization C . In the case α < 1 this relation gives β = θ. Now comparing (51) and (13) we see that The sign of equality is achieved here only in the case β = 0.
Finishing this paper the following should be noted. In the previous article [3] it was noted that when calculating the integral in the representation (3) numerical integration algorithms have difficulties in the domain of small values of the coordinate x, in the domain of large values of the coordinate x and in the domain of values of the characteristic parameter α ≈ 1. The first problem out of these three ones has been solved in this paper. The assumption made in the work [3] that the cause of the problem is associated with the behavior of the integrand in (3), was correct. Indeed, this integrand at smaller values of x starts acting as a singular function which makes it impossible to use for numerical algorithms to calculate the integral of it. Therefore, to calculate the density in the indicated domain x it is necessary to use series expansions of the density. The reason for the calculation difficulties in the second case is also the behavior of the integrand in the representation (3). As shown in section 5 with large x it behaves as a singular function. Therefore, to calculate the density in this domain of the coordinate, it is also expedient to use expansions of the density in a series. To solve the third problem, one can use the method proposed in the paper [39]. In this paper, to calculate the density, the authors propose to use the series expansion of a strictly stable law in view of the parameter α. However, the solution of each of the remaining two problems requires further research, which is beyond the scope of this paper. It was mentioned in the introduction that the obtained expansions of the probability density of a strictly stable law are useful in problems related to the calculation of the probability density of a fractionally stable law. Indeed, the probability density of the fractionally stable law q(x, α, β, θ) is determined by the Mellin transform of two strictly stable laws q(x, α, β, θ) = ∞ 0 g(xy β/α , α, θ, λ)g(y, β, 1)y β/α dy, where g(x, α, θ) and g(y, β, 1) -densities of strictly stable and one-sided strictly stable laws with the characteristic function (1) [2,26,27]. Here, the characteristic exponents α and β vary within 0 < α 2 and 0 < β 1, the asymmetry parameter θ takes the values within the interval |θ| min(1, 2/α − 1) and λ > 0 is the scaling parameter.
From the formula (52) it is clear that it is necessary to be able to calculate densities of the strictly stable laws g(y, α, θ, λ) and g(y, β, 1) for the calculation of density q(x, α, β, θ, λ). The integral representation (3) is used to calculate these densities. In this regard, at this stage, certain problems may arise with the calculation of the improper integral in (52). Indeed, to calculate the integral in (52) the numerical integration algorithm calculates the integrand at some integration nodes y i . If it turns out that the next integration node y i is less than the value of the critical coordinate y cr , i.е. |y i | < y cr , then the numerical integration algorithm will be unable to calculate the densities g(y i , α, θ) and g(y i , β, 1). This will lead to an incorrectly calculated density value q(x, α, β, θ). To eliminate this problem and to calculate the densities g(y i , α, θ) and g(y i , β, 1) in the case |y i | < y cr it is expedient to use the expansion (12). Thus, the use of the expansion (12) will give an opportunity to exclude the integration error associated with the singular behavior of the integrand in the representation (3) at small values of y.
Fractionally stable distributions turn out to be a convenient tool to describe the probability density distribution of gene expression obtained by means of NGS technology. Fig. 10 shows the approximation results of the probability density distribution of gene expression obtained with NGS technology using the density of a fractional stable law. In these figures, the dots are the histogram of the probability density of gene expression, the solid curve is the density (52). The density parameters q(x, α, β, θ, λ) are shown in the figures. As we can see from these figures the density q(x, α, β, θ, λ) approximates the experimental data quite well in a very wide range of values.
It should be noted that the parameters α, β, θ, λ, given in Fig. 10 were estimated from the experimental data using the minimum distance method, which is based on the distance χ 2 [40]. However, this method of parameter estimation is not effective. To build an effective estimate of the parameters of fractionally stable distributions, it is necessary to build an estimate based on the maximum likelihood method. Until now, the creation of such an estimate has met with some difficulties.They are related to the fact that to calculate the likelihood function, it is necessary to calculate the density of the fractionally stable law at the points determined by the original data sample. Taking into account that the density q(x, α, β, θ, λ) is calculated according to the formula (52), then at small values of the coordinate x or or the integration variable y the algorithm numerical integration gave the wrong result. This, in turn, led to an error in estimating the distribution parameters. Thus, the use of the expansion (12) when calculating the integral in (52) will give an opportunity to calculate the density correctly, which in turn will allow implementing the methods for estimating the parameters of fractionally stable and strictly stable distributions based on the maximum likelihood method.