Determining the First Probability Density Function of Linear Random Initial Value Problems by the Random Variable Transformation (RVT) Technique: A Comprehensive Study

and Applied Analysis 3 Table 1: List of the thirteen different cases considered to perform the full study depending on whether the IVP is homogeneous (H) or nonhomogeneous (NH) and the form that uncertainty is considered into the problem: just through one RV:Z 0 , B, orA (Cases I-II-III.1, I-II-III.2, and III.3), two RVs: (Z 0 , B), (Z 0 , A), or (B,A) (Cases I-II.3 and III.4-5-6), and the three involved RVs: (Z 0 , B, A) (Case III.7). Type IVP Case


Introduction and Motivation
Over the last few decades, random differential equations (RDEs) have been demonstrated to be powerful tools to model numerous problems appearing in many different areas such as physics, engineering, economics, epidemiology, and hydrology.The consideration of randomness into their formulation through initial/boundary conditions, source terms, and/or coefficients adapts better than their deterministic counterpart to model the uncertainty associate to the experimental measurement required to set the above inputs as well as the inherent complexity involved in many real modelling problems.This approach leads to face new and exciting questions different from the corresponding ones appearing in the deterministic scenario.Indeed, instead of obtaining just the solution stochastic process (SP) of RDEs, the theory is also concerned with its probabilistic properties, mainly the computations of the expectation and variance functions.The computation of the first probability density function  of the solution SP, say ( ), is much more desirable since, from it, one can compute the previous statistical functions as simple particular cases and, in addition, it provides a comprehensive probabilistic description of the solution SP for each time instant .However, the computation of the 1-PDF constitutes a major challenge that can been achieved in only a few cases.
The aim of this paper is to determine the 1-PDF, 1 ( , ), of the solution SP ( ) to the linear random initial value problem (IVP): where the data 0 , and are assumed to be continuous random variables (continuous RVs) defined on a common respectively.Hereinafter, in order to avoid cumbersome notation, we will hide the sample dependence when writing domains of continuous RVs.In this way, for instance, the domain 0 will be written as 0 = { 0 : 0,1 ≤ 0 ≤ 0,2 } rather than the first expression in (2).The same can be said for , , and the domain of any other RV throughout this paper.We allow the left (right) endpoint of each interval of the domain takes the value −∞ (+∞); that is, we also consider unbounded continuous RVs.Throughout the paper, we will denote by 0 ( 0 ), ( ), and ( ) the PDFs of the continuous RVs 0 , , and , respectively.The case where 0 , and are pairwise dependent continuous RVs will also be treated.In such case, 0 , ( 0 , ), 0 , ( 0 , ), and , ( , ) will denote the joint PDFs of the random vectors: ( 0 , ), ( 0 , ), and ( , ), respectively.Finally, we will also deal with the case where 0 , , and are dependent continuous RVs, then 0 , , ( 0 , , ) will represent their joint PDF.Notice that the domains of these two-and three-dimensional PDFs often can be written directly as products of the sets 0 , , and given by (2).In order to compute the 1-PDF 1 ( , ), random variable transformation (RVT) method will be applied.RVT is a probability technique that allows us to calculate the PDF ( ) of a RV resulting after the algebraic transformation of another RV, say , whose PDF, ( ), is known.In its simplest scalar formulation, the method reads as follows: if is a continuous RV lying on the domain or support = { : 1 ≤ ≤ 2 }, whose PDF is ( ) > 0 and = ( ) being : ⊆ R → R a monotone mapping on , then where ( ) = is the inverse function of on , which is assumed to have a continuous derivative on and |d ( )/d | denotes the modulus of the derivative of ( ).In the particular case that increases (decreases) on , the domain of = ( ) is determined by study a linear oscillator assuming randomness just in the two initial conditions related to the position and velocity.Most of the subsequent contributions have focused on the study of particular equations assuming specific probabilistic distributions for the involved uncertainty which facilitates the analysis.Here, we point out some recent contributions that illustrate quite well the current trends of RVT method in dealing with RDEs.In [2], authors solve the radiative transfer equation in a semi-infinite continuous stochastic medium with Rayleigh scattering.RVT method is applied to obtain the 1-PDF of the solution when the optical depth space variable is assumed to be a RV belonging to the following particular distributions: exponential and Gaussian.Higher order statistical moments of the solution stochastic process are also computed.An analogous study on the stochastic transport equation of neutral particles with anisotropic scattering can be found in [3].In [4,5], authors apply RVT technique to develop a stochastic finite element method for solving some stochastic problems with random excitation.
The application of RVT technique to the exact determination of the 1-PDF of the solution SP of RDEs requires the previous computation of the exact solution of the RDE under study.However, in the outstanding contribution [6], author takes advantage of RVT method together with classical numerical techniques to illustrate, through a wide range of examples, the potentiality of this method to approximate the 1-PDF for the solution SP of some RDEs.
As it has been announced previously, in this paper, we will compute the 1-PDF of the IVP (1) whose exact solution is available.For the sake of clarity, in the presentation, we will divide the study in the three main IVPs (I)-(III) listed in Table 1.In Case (I), we consider the homogeneous (H) problem whereas Cases (II) and (III) deal with the nonhomogeneous (NH) cases.Within each case, we distinguish, in a systematic manner, the different possibilities regarding the randomness of each of the involved parameters 0 , , and .These casuistries include the situations where the parameters are statistically dependent.In this context, and as it has been pointed out previously, if for example 0 and are statistically dependent, then 0 , ( 0 , ) will denote the joint PDF of the random vector ( 0 , ), and the same can be said for the rest of the possibilities.The IVPs (I) and (II) can be seen as particular cases of the IVP (III) when P[{ ∈ Ω : ( ) = 0}] = 1 and P[{ ∈ Ω : ( ) = 0}] = 1, respectively.When uncertainty can only be attributted to 0 and , and the parameter can be set in a deterministic way, it is more realistic and convenient to assume that the joint PDF 0 , ( 0 , ) is known rather than 0 , , ( 0 , , ).Notice that the construction of the joint PDF of only two continuous RVs from measured data can become a difficult problem which accuracy can deteriorate severely if one includes a new and inappropriate RV into its formulation [7].The most accuracy of the PDF of the random input parameter, the best approximation of the 1-PDF of the solution SP of the IVP (1).Therefore, the consideration of all the thirteen separate cases listed in Table 1 turns out more recommendable from a practical point of view.
For the sake of clarity in the development of all the cases listed in Table 1, throughout this paper, the input parameters, Table 1: List of the thirteen different cases considered to perform the full study depending on whether the IVP is homogeneous (H) or nonhomogeneous (NH) and the form that uncertainty is considered into the problem: just through one RV: 0 , , or (Cases I-II-III.1,I-II-III.2,and III.3), two RVs: ( 0 , ), ( 0 , ), or ( , ) (Cases I-II.3 and III.4-5-6), and the three involved RVs: ( 0 , , ) (Case III.7).

Type
IVP Case which are assumed to be continuous RVs, will be denoted by upper cases, while deterministic magnitudes will be written by lower cases.More precisely, for instance, in Case III.4,we will denote by 0 and the random inputs, while the multiplicative coefficient in the RDE will be denoted by rather than .
The paper is organized as follows.In Section 2 we summarize the main results concerned with RVT that will be applied throughout the paper.We particularly state different versions of this useful technique including its application to different transformations that will facilitate the presentation of the results.Sections 3, 4, and 5 provide a detailed study where the 1-PDF of the solution SP, ( ), corresponding to the IVPs (I), (II), and (III) listed in Table 1, respectively, is computed.For each one of the thirteen casuistries, the study shows, at least, an illustrative example.The choice of the PDFs considered in the examples has been made to show the ability of the method to deal with both standard and non-standard dependent probability distributions.In Section 6 we include some considerations related to the application and better understanding of RVT method that we found particularly useful.Conclusions are drawn in the closing section.

Preliminaries
Below, we state several versions of the RVT technique as well as some related results emerging from its applications that will play a relevant role in our subsequent developments.Most of these results can be found in [8,9] or they are a direct consequence of them.

Theorem 1 (RVT technique: scalar version).
Let be a continuous RV with PDF ( ) and domain = { : ( ) > 0}.Let = ( ) be a new RV generated by the map : R → R which is assumed to be continuously differentiable on and such that ὔ ( ) ̸ = 0 except at a finite number of points.
Let one suppose that, for each ∈ R, there exist ( ) ≥ 1 points: 1 ( ), 2 ( ), . . ., ( ) ( ) ∈ such that Then Although Theorem 1 unifies the treatment of the different cases that one can present to compute the PDF, ( ), in practice, it is easier to determine ( ) by dividing the domain of RV into subintervals where the mapping is monotone and then applying formula (3) on each subinterval.The process to compute ( ) on the whole domain of is completed by adding the corresponding expressions calculated previously for each subinterval.
In the simplest but significant case where the map is linear, Theorem 1 reads as follows.
The following result is a direct application of Theorem 1 in the case that ( ) = exp( ) + .

Proposition 3 (RVT technique: exponential transformation).
Let be a continuous RV with domain = { : 1 ≤ ≤ 2 } and PDF ( ).Then the PDF ( ) of the exponential transformation = exp( ) + , with ̸ = 0, is given by If = 0 or = 0, then = + with probability 1 and The computation of the joint PDF of two or more continuous RVs using the RVT method can also be performed by using the following generalization of formula (3).

Theorem 4 (RVT technique: multidimensional version)
. Let X = ( 1 , . . ., ) be a random vector of dimension with joint PDF X (x).Let r : R → R be a one-to-one deterministic map which is assumed to be continuous with respect to each one of its arguments and with continuous partial derivatives.Then, the joint PDF Y (y) of the random vector Y = r(X) is given by where s(y) is the inverse transformation of r(x): x = r −1 (y) = s(y) and is the Jacobian of the transformation; that is, which is assumed to be different from zero.
As we will see later, the analysis of Cases I-3, II-3, and III-3-6 requires the computation of the PDF of the sum and product of two continuous RVs which turns out by the application of Theorem 4 in its two-dimensional version.Thus, for the sake of clarity in the exposition, we specialize Theorem 4 in this significant case.

Theorem 5 (RVT technique: two-dimensional version). Let
be a one-to-one deterministic map from R 2 to R 2 ; that is, there exists its inverse transformation: on the range of the map (12).Let one assume that both maps (12) and (13) are continuous.Let further assume that the following partial derivatives exist and are continuous and the Jacobian 2 of the inverse map satisfies on the range of the transformation (12).Then, the joint PDF Next, we apply Theorem 5 in the particular case that transformation 1 only depends on variable 1 and 2 only depends on variable 2 .As it will be seen later, this result will be crucial in further applications.Proposition 6.Let X = ( 1 , 2 ) be a two-dimensional RV with joint PDF 1 , 2 ( 1 , 2 ).Let be a one-to-one deterministic map from R 2 to R 2 ; that is, there exists its inverse transformation: on the range of the map (17).Let one assume that both maps (17) and (18) are continuous and the two following derivatives that exist are continuous and satisfy on the range of the transformation (17).Then, the joint PDF On the other hand, applying Theorem 5 to we obtain the PDF of the sum and product of two continuous RVs, respectively.We state both results in the two following propositions.
Let ( 1 , 2 ) be a continuous random vector with joint PDF or, equivalent, by If 1 and 2 are independent continuous RVs, since , where ( ) denotes the PDF of , = 1, 2, (21) and (22) inform one that the PDF of the sum of two independent continuous RVs is just the convolution of their respective PDFs: Proposition 8 (RVT technique: product of two continuous RVs).Let ( 1 , 2 ) be a continuous random vector with joint PDF 1 , 2 ( 1 , 2 ) with respective domains Equivalently, if If 1 and 2 are independent continuous RVs with PDF's 1 ( 1 ) and 2 ( 2 ), respectively, then (24) and (25) become respectively.
As usual we have not specified the domain of variation of 1 in ( 24) and (25) since it is cumbersome.However, later we will detail it in some illustrative cases where it appears (see for instance Case III.6).
We close this section by extending Proposition 7 for the case of three terms since it will be required to deal with Case III.7.This result comes directly from the application of Theorem 4.

Case Study: Initial Value Problem (I)
This section is devoted to obtain the 1-PDF 1 ( , ) of the solution SP ( ) to the IVP (I) in each of the three cases listed in Table 1.Notice that ( ) has the following expression: 3.1.Case I.1: 0 Is a Random Variable.As we pointed out previously for the sake of clarity in the presentation, we rewrite (30) by distinguishing the deterministic character of parameter (which is written with a lower case letter): Next, we first fix : ≥ 0 and denote = ( ).Then we apply Proposition 2 to Then, taking into account that the domain of RV 0 is given by (2), one gets where We illustrate the previous development in the following example where 0 is assumed to be a standard continuous RV, although further distributions, not necessarily standard, could be considered.

Case I.2:
Is a Random Variable.First, we rewrite (30) by highlighting the deterministic character of the initial condition 0 : Next, we fix : > 0 and denote = ( ).In the first analysis, we will assume that 0 ̸ = 0. Then we apply Proposition 3 to Then, taking into account that / 0 = e ( − 0 ) > 0 and the domain of RV is given by ( 2), one gets where According to (2), if 0 = 0, then For = 0 , from (36) ( ) = ( 0 ) = 0 , which is a deterministic initial condition.Then its 1-PDF can be written through the Dirac delta function as follows: (41) In the example below, we illustrate the previous development in the case where 0 > 0.
3.3.Case I.3: ( 0 , ) Is a Random Vector.Let us denote by 0 , ( 0 , ) the joint PDF of the random vector ( 0 , ).Now, we rewrite (30) in the following equivalent form: In order to apply RVT, we fix : > 0 and denote 1 = 1 ( ), first we will determine the joint PDF of 1 and 2 by applying Proposition 6 to Then taking into account one gets Once the joint PDF of ( 1 , 2 ) has been determined, the computation of the PDF of = 1 2 follows directly by applying Proposition 8. Indeed, as 2 = e ( − 0 ) ̸ = 0, we will apply formula (25 or equivalently by using (46): where Notice that if = 0 , then ( ) = ( 0 ) = 0 and 1 ( , is just the PDF of RV 0 , which is easily obtained from the datum 0 , ( 0 , ) as the following marginal distribution: Example 12. Let ( 0 , ) be a two-dimensional RV whose PDF is defined by Notice that we are implicitly assuming independence between 0 and since the joint PDF 0 , ( 0 , ) factorizes as the product of the marginal PDFs 0 ( 0 ) = 2 0 , 0 < 0 < 1, and ( ) = 2 , 0 < < 1.Then, taking into account (48)-( 50) and (52), the 1-PDF of ( ) is given by ( 2 ) For the sake of clarity in the graphical representation of 1 ( , ), Figure 3 shows two equivalent plots of 1 ( , ) in the case that 0 = 0. Notice that the plot on the left side is the resultant surface where some 1-PDF 1 ( , ) have been highlighted for different fixed times, whereas these 1-PDFs have been represented in these fixed times on the right side.As the two integrals appearing into the expression (53) can be computed explicitly, an equivalent expression to 1 ( , ) with 0 = 0 is given by (54) 5} (corresponding to the solid lines in both plots) in the case that ( 0 , ) is a two-dimensional RV whose PDF is given by ( 52) and 0 = 0.

Case Study: Initial Value Problem (II)
This section is addressed to determine the 1-PDF 1 ( , ) of the solution SP ( ) to the IVP (II) in each of the three cases listed in Table 1.Notice that ( ) has the following expression: 4.1.Case II.1: 0 Is a Random Variable.As we did in Case I.1, for the sake of clarity in the presentation, we rewrite (56) by emphasizing the deterministic character of parameter (which is written with a lower case letter): Next, we fix : ≥ 0 and denote = ( ).Then we apply Proposition 2 to Taking into account that the domain of RV 0 is given by (2), one gets where Example 13.Let us assume that 0 has a gamma distribution, 0 ∼ Ga( ; ), , > 0. Therefore, according to (59)-( 60), the 1-PDF of ( ) is given by: where Γ( ) denotes the deterministic gamma special function.For each ≥ 0 , the domain of has been determined taking into account in (60) that, in this case, 0,1 = 0 and 0,2 = +∞.It can be checked that 1 ( , ) is a PDF for each ≥ 0 .Figure 4 shows 1 ( , ) at different values of in the particular case that 0 ∼ Ga(2; 1), = 1, and 0 = 0.

Case II.2:
Is a Random Variable.First, we rewrite (56) by highlighting the deterministic character of the initial condition 0 : Next, we fix : > 0 and denote = ( ).Then we apply Proposition 2 to This yields where For = 0 , as it also occurred in Case I.2, ( ) = ( 0 ) = 0 and therefore Example 14.Let us assume that has a 2 -distribution with ] degrees of freedom, ∼ 2 (]), and ] > 0. Therefore, according to (64)-(65), the 1-PDF of ( ) is given by For each : > 0 , the domain of has been determined taking into account in (65) that in this case 1 = 0 and 2 = ∞.

Case Study: Initial Value Problem (III)
This section deals with the computation of the 1-PDF 1 ( , ) of the solution SP ( ) to the IVP (III) in each of the seven cases listed in Table 1.Notice that ( ) has the following expression: Throughout the sequent analysis depending on whether is considered to be a RV (Cases III.3, III.5, III.6, and III.7) or a deterministic constant (Cases III.1, III.2, and III.4), see Table 1; we will assume that P[{ ∈ Ω : ( ) ̸ = 0}] = 1 or ̸ = 0, respectively.
Before closing this case, we provide an example where the usefulness of the 1-PDF 1 ( , ) to determine the statistical moments of the solution ( ) and to compute the probability of sets of interest is shown.To illustrate these applications, we will choose the context of Example 16, although it could be applied to any example throughout this paper.
Example 17.Let us consider the context of Example 16 where the 1-PDF of ( ), 1 ( , ), has been computed (see expression (86)).Then, the statistical moment of order of ( ) with respect to the origin can be computed directly in terms of 1 ( , ) as follows: As a consequence, the following expressions for the mean and the variance of ( ) are obtained: Figure 8 shows the expectation and variance of ( ) for the same values of , 0 , , and considered in Example 16.

Case III.3:
Is a Random Variable.So far we have observed that the application of RVT method to determine the PDF of a RV, say , generated by a transformation = ( ), relies strongly on the feasibility of computing the inverse of the map = ( ).Fortunately, this has been done exactly in the previous cases, but, in general, it is not possible in the current case where just is assumed to be a RV.In fact, once : > 0 has been fixed, we must isolate in the equation: which is not possible to perform in an exact manner.To circumvent this drawback, we will apply the Lagrange-Bürmann formula which gives the Taylor series expansion of the inverse of an analytic function.
Theorem 19 (Lagrange-Bürmann formula, see [10]).Suppose that is defined as a function of the variable by an equation of the form: = ( ) where is analytic about the point 0 and ὔ ( 0 ) ̸ = 0.Then, it is possible to invert (or to solve) the equation for : = ( ) on a neighbourhood N( ( 0 ); ), > 0 of ( 0 ): Although this result permits obtaining, from a theoretical stand point, the inverse function of ( ), in practice, often this can only be achieved in an approximate manner since the infinite series (97) must be truncated to be kept computationally feasible.Moreover, the representation (97) of the inverse is only valid in a certain neighbourhood N( ( 0 ); ), > 0 of ( 0 ), whose diameter > 0 must be determined carefully in each case study.
Let us apply Theorem 19 to determine the 1-PDF 1 ( , ) of ( ).Fixing : > 0 and denoting = ( ), we consider the map (96) which is analytic about any numerical value 0 ̸ = = ( ), ∈ Ω on the assumed domain to the RV .As the map (96), in general, is not monotone,we    have to apply Theorem 1.In a first step, we divide the domain of the map (or equivalently, the domain of the RV ) into subintervals: A 1 , A 2 , . . ., A where is monotone.Then, we will fix a such subinterval A =] ini, , end, [, 1 ≤ ≤ where the contribution to the total 1-PDF 1 ( , ) is going to be calculated.For this, we select a point 0, ∈ A where condition ὔ ( 0, ) ̸ = 0 is met.By applying the Lagrange-Bürmann formula, we construct the inverse of the map ( ) = ( ) on A that, in the sequel, will be denoted by ( ): In practice, this infinite series could only converge on a subset of ] ( ini, ), ( end, )[ (if increases on A ) or ] ( end, ), ( ini, )[ (if decreases on A ).In such case, the function ( ) can be completed on the whole interval by taking another (or other, if necessary) appropriate point(s), say 0 , ∈ A , 1 ≤ ≤ and then repeating the above process.In this manner, a piecewise inverse function ( ) on A will be constructed.
In accordance with RVT method, besides constructing the inverse ( ), one requires to compute its derivative.We again take advantage of Lagrange-Bürmann formula to complete this computation.In fact, notice that once ( ) has been constructed, from (98), one gets d ( ) Notice that if ( ) has been defined by means of a piecewise function, then d ( )/d will also be defined in the same way.The process will be culminated by repeating again the previous argument on all the subintervals A 1 , A 2 , . . ., A .Following the previous development, the 1-PDF 1 ( , ) of ( ) can be computed as follows: where ( ) and d ( )/d are defined by ( 98) and (99), respectively.
Often, the infinite series (98) has to be truncated at the term to control computational burden.In this way, we obtain an approximation of the inverse: Thus, an approximation of its derivative is d , ( ) The infinite series (98) and (99) have the same convergence radius; however, often in practice, we need to handle their corresponding truncations (101) and (102), respectively.In this case, the selection of the appropriate 0 , depends also on the quality of the approximations provided by both truncated series.
Repeating the foregoing process on each interval A , 1 ≤ ≤ , one gets the corresponding approximation of 1 ( , ) given by In the following example, we illustrate the previous development.
Once the joint PDF of ( 1 , 2 ) has been determined, the computation of the PDF of = 1 + 2 follows directly by applying Proposition 7. Indeed, as 2 = − / ̸ = 0 (since into Case III, ̸ = 0; otherwise we would be in Case I), we will apply formula (22 where 2 denotes the domain of variation of 2 which, according to (120)-(124), depends on both the sign of the products 1 2 and .Using (119), it is equivalent to The range of variation of , denoted by , can be straightforwardly computed taking into account = 1 + 2 and the domains of 1 and 2 which have been determined previously in (120)-(124).As its practical determination in specific examples is simple from the previous exposition, in order to avoid an unwieldy notation, we do not rewrite the final general expression.
Finally, we consider the case = 0 that implies ( ) = ( 0 ) = 0 .Therefore, Example 23.With the aim of showing the generality of the obtained results to deal with different classes of continuous RVs, in this example, we will consider that the joint PDF of the input continuous RVs and is constructed by a copula function.Let and be two uniform continuous RVs on the interval ]0, 1[; that is, , ∼ Un(]0, 1[).Using the Farlie-Gumbel-Morgenstern copula [11], we construct the random vector ( , ) with joint PDF hence, and are statistically dependent.As a characteristic of copula functions, notice that the two marginal distributions of , ( , ) are just the individual distributions of and .In the following, we will consider the case previously developed for 0 = 0 and we will take 0 = 0. First, notice that according to (141) and taking into account that 0 < , < 1, the domains of variation of 1 and 2 are 1 ∈]0, ∞[ and 2 ∈]0, e − 1[, respectively; thus, = 1 2 lies on ∈]0, ∞[.However, we must refine the domain of integration 2 to 2 of (148), in such a way that the two arguments of the PDF , (⋅, ⋅) lie inside the interval ]0, 1[, where the input continuous RVs and are defined.Thus, fixing > 0 and > 0, we must determine 2 such as  the probabilistic behaviour of the solution SP ( ) at each time instant , it must be pointed out that the computation of higher PDFs is also possible and desirable.In fact, for instance, the 2-PDF, say, 2 ( 1 , 1 ; 2 , 2 ) provides a whole probabilistic information of ( ) in two time instants 1 and 2 .Therefore, from it, we can calculate relevant probabilistic properties such as the correlation function (164) As it can be guested, from the previous development, the computation of 2 ( 1 , 1 ; 2 , 2 ) follows in broad outline that 1 ( , ).In Example 27, we compute the 2-PDF within the framework of foregoing Case II.3 including the correlation function.
Example 27.Let us fix 1 , 2 such as 2 > 1 ≥ 0 and denote 1 = ( 1 ) and 2 = ( 2 ).We want to determine the joint PDF of continuous RVs 1 and 2 .For it, we will apply Theorem 5 to where the expression of the solution SP given by (56) has been considered.Now, taking into account that then, the involved Jacobian is given by

Figure 14 :
Figure14: Plots of 1 ( , ) given by (159) in Example 24 at different values of with 0 = 0 in the case that ( 0 , , ) is a three-dimensional Gaussian vector with mean and variance-covariance matrix given by (161).