High Accurate Simple Approximation of Normal Distribution Integral

Abstract

The integral of the standard normal distribution function is an integral without solution and represents the probability that an aleatory variable normally distributed has values between zero and 𝑥. The normal distribution integral is used in several areas of science. Thus, this work provides an approximate solution to the Gaussian distribution integral by using the homotopy perturbation method (HPM). After solving the Gaussian integral by HPM, the result served as base to solve other integrals like error function and the cumulative distribution function. The error function is compared against other reported approximations showing advantages like less relative error or less mathematical complexity. Besides, some integrals related to the normal (Gaussian) distribution integral were solved showing a relative error quite small. Also, the utility for the proposed approximations is verified applying them to a couple of heat flow examples. Last, a brief discussion is presented about the way an electronic circuit could be created to implement the approximate error function.

1. Introduction

Normal distribution is considered as one of the most important distribution functions in statistics because it is simple to handle analytically, that is, it is possible to solve a large number of problems explicitly; the normal distribution is the result of the central limit theorem. The central limit theorem states that in a series of repeated observations, the precision of the approximation improves as the number of observations increases [1]. Besides, the bell shape of the normal distribution helps to model, in a practical way, a variety of random variables. The normal (or gaussian) distribution integral has a wide use on several science branches like: heat flow, statistics, signal processing, image processing, quantum mechanics, optics, social sciences, financial mathematics, hydrology, and biology, among others. Normal distribution integral has no analytical solution. Nevertheless, there are several methods [2] which provide an approximation of the integral by numerical methods: Taylor series, asymptotic series, continual fractions, and some other more. Numerical integration is expensive in computation time and it is only good for numerical simulations; in several cases the power series require many terms to obtain an accurate approximation; this fact makes them a numerical tool exclusive for computational calculations, which are not adequate for quick calculations done by hand.

The Homotopy Perturbation Method [3–10] was proposed by He, it was introduced like a powerful tool to approach various kinds of nonlinear problems. As is well known, nonlinear phenomena appear in a wide variety of scientific fields, such as applied mathematics, physics, and engineering. Scientists in these disciplines are constantly faced with the task of finding solutions of linear and nonlinear ordinary differential equations, partial differential equations, and systems of nonlinear ordinary differential equations. In fact, there are several methods employed to find approximate solutions to nonlinear problems such as variational approaches [11–16], Tanh method [17], exp-function [18], Adomian’s decomposition method [19, 20], and parameter expansion [21]. Nevertheless, the HPM method is powerful, is relatively simple to use, and has been successfully tested in a wide range of applications [15, 22–29]. The homotopy perturbation method can be considered as combination of the classical perturbation technique and the homotopy (whose origin is in the topology), but has eliminated the limitations of the traditional perturbation methods [30]. The HPM method does not need a small parameter or linearization, in fact it only requires few iterations for getting highly accurate solutions. Besides, the HPM method has been used successfully for solving integral equations, for instance, the case of the Volterra integral equations [31]. The method requires an initial approximation, which should contain as much information as possible about the nature of the solution. That often can be achieved through an empirical knowledge of the solution. Therefore, we propose the use of the homotopy perturbation method to calculate an approximate analytical solution of the normal distribution integral.

In this work the HPM method is applied to a problem having initial conditions. Nevertheless, it is also possible to apply it to the case where the problem has boundary conditions; here, the differential equation solutions are subject to satisfy a condition for different values of the independent variable (two for the case of a second order equation) [32].

This paper is organized as follows. In Section 2, we solve the normal distribution integral (NDI) by HPM method. In Section 3, we use the approximated solution of NDI to establish an approximate version of error function, and the result is compared to other analytic approximations reported in the literature. In Section 4, a series of normal distribution related integrals are solved. In Section 6, two application problems for the error function are solved. In Section 7, we summarize our findings and suggest possible directions for future investigations. Finally, a brief conclusion is given in Section 8.

2. Solution of the Gaussian Distribution Integral

This section deals with the solution of the Gaussian distribution integral; it is represented as follows:𝑦(𝑥)=^𝑥₀exp−𝜋𝑡²𝑑𝑡,𝑥∈ℜ. (2.1)

This integral can be reformulated as a differential equation𝑦^′−exp−𝜋𝑥²=0,𝑥∈ℜ. (2.2)

Evaluating the integral (2.1) at 𝑥=0, the result will be zero, therefore, the initial condition for the differential equation (2.2) should be 𝑦(0)=0.

To apply the homotopy perturbation method, the next equation is created1−𝑝𝐺𝑦−𝐺𝑦_𝑖+𝑝𝑦^′−exp−𝜋𝑥²=0, (2.3) where 𝑝 is the homotopy parameter. The solution for (2.2) is similar, qualitatively, to a hyperbolic tangent because when 𝑥 tends to ±∞, the derivative 𝑦^′ tends to zero, hence by symmetry, 𝑦 tends to the same constant (absolute value) on both directions. Therefore, it is desirable that the first approach of the HPM method contains a hyperbolic term. In consequence, a differential equation is established (𝐺(𝑦)) that may be solved using hyperbolic tangent, which is given by𝐺𝑦=1+𝑐²𝑥²𝑑𝑦(𝑥)−𝑦(𝑥)𝑑𝑥𝑑−²𝑑𝑎+𝑎𝑐²𝑥²+𝑏𝑐. (2.4)

Also, the initial approximation for homotopy would be𝑦_𝑖(𝑥)=𝑑tanh(𝑎𝑥+𝑏arctan(𝑐𝑥)), (2.5) where 𝑎, 𝑏, 𝑐, and 𝑑 are adjustment parameters.

Now, by the HPM it is assumed that (2.3) has the following form:𝑦(𝑥)=𝑦₀(𝑥)+𝑝𝑦₁(𝑥)+𝑝²𝑦₂(𝑥)+…. (2.6)

Adjusting 𝑝=1, the approximate solution is obtained𝑦(𝑥)=𝑦₀(𝑥)+𝑦₁(𝑥)+𝑦₂(𝑥)+…. (2.7)

Substituting (2.6) into (2.3) and equating terms with powers equals to 𝑝, it can be solved for 𝑦₀(𝑥), 𝑦₁(𝑥), 𝑦₂(𝑥), and so on. In order to fulfill the initial condition from (2.2) (𝑦(0)=0), it follows that 𝑦₀(0)=0,𝑦₁(0)=0, 𝑦₂(0)=0, and so on. Thus, the result is𝑦₀(𝑥)=𝑑tanh(𝑎𝑥+𝑏arctan(𝑐𝑥)), (2.8) where 𝑎=−39/2, 𝑏=111/2, 𝑐=35/111, and 𝑑=1/2; these parameters are adjusted (using the NonlinearFit command from Maple Release 15) to obtain a good approximation; this adjustment allows to ignore successive terms. (Given a total of 𝑘-samples from the exact model, the NonlinearFit command finds values of the approximate model parameters such that the sum of the squared 𝑘-residuals is minimized.) Therefore, the proposed approximation for (2.1) is𝑦(𝑥)=^𝑥₀exp−𝜋𝑡²1𝑑𝑡≈2tanh39𝑥2−1112arctan35𝑥111,−∞≤𝑥≤∞. (2.9)

Also, it is known that solution for (2.1) may be expressed in terms of the error function𝑦(𝑥)=12erf√𝜋𝑥, (2.10) where the error function is defined by(𝑥)=2erf√𝜋^𝑥₀exp−𝑡²𝑑𝑡. (2.11)

Figure 1(a) shows the approximate solution (2.9) (continuous line) and the exact solution (2.1) (diagonal cross), and Figure 1(b) shows the relative error curve.

(a)

(b)

Figure 1 

(a) Gaussian distribution integral and approximation. (b) Relative error.

In Figure 1 it can be seen that the maximum relative error is negligible for many practical applications in fields like engineering and applied sciences. In fact, the precision of this approximation allows to derive (2.9) with respect to 𝑥 and graphically compare it to the exact function exp(−𝜋𝑥²), as shown in Figure 2(a), reaching an acceptable relative error (see Figure 2(b)) in the range of −1.3≤𝑥≤1.3. Therefore,𝑦^′(𝑥)≈316428+15925𝑥²exp(−39𝑥+111arctan(35𝑥/111))12321+1225𝑥²exp(−39𝑥+111arctan(35𝑥/111))+1²,−∞≤𝑥≤∞. (2.12)

(a)

(b)

Figure 2 

(a) (2.11) and exact function exp(−𝜋𝑥2). (b) Relative error.

From this approximation it is possible to generate normal distribution functions like the error function or the cumulative distribution function.

3. Error Function

From (2.9) and (2.10) it can be concluded that(𝑥)𝑥erf=2𝑦√𝜋≈tanh39𝑥2√𝜋−1112arctan35𝑥111√𝜋,−∞≤𝑥≤∞. (3.1)

Now, function (3.1) is readjusted using the NonlinearFit command, then the command “convert” (with option “rational”) (by using Maple 15) is applied in order to obtain an optimized version of (3.1):erf(𝑥)≈tanh4907𝑥−44617753arctan34𝑥191,−∞≤𝑥≤∞. (3.2)

Next, three different approximations for the error function are presented; all of them are compared to the approximations proposed in this work ((3.1) and (3.2)).

(1) In [33] it was presented the following approximationerf(𝑥)≈−41−exp𝜋+0.14𝑥²1+0.14𝑥²𝑥²,𝑥≥0. (3.3)

Due to the fact that (3.3) is only useful for values 𝑥⩾0, it is considered a disadvantage because it limits variation over 𝑥 axis.

(2) In [34], it was reported an approximation related to the function error, which is expressed as𝐹(𝑥)=^∞_𝑥−𝑡exp²2𝑑𝑡,𝑥∈ℜ. (3.4)

The approximation for (3.4) is𝑃(𝑥)=𝑃₀+𝑥⁻¹−𝑥exp²2−𝑃²₀𝑥²+exp−𝑥²/21+𝑏𝑥^21/21+𝑎𝑥^21/2, (3.5) where𝑃₀=𝜋2^1/2,𝑎=1+1−2𝜋²+6𝜋^1/22𝜋,𝑏=2𝜋a². (3.6)

In order to compare our work with  (3.5), consider𝐹(𝑥)=^∞_𝑥−𝑡exp²2𝑑𝑡=−𝜋2𝑥erf√2+𝜋2. (3.7)

Replacing 𝑃(𝑥) for 𝐹(𝑥), it is obtained that𝑃(𝑥)≈−𝜋2𝑥erf√2+𝜋2, (3.8)

for solving erf(𝑥/√2), 𝑥erf√2≈1−2𝜋𝑃(𝑥). (3.9)

From (3.9), the conclusion iserf(𝑥)≈1−2𝜋𝑃√2𝑥, (3.10) where 𝑃(√2𝑥) is calculated from (3.5).

(3) In [2] was reported an approximation of the normal distribution integral, which is transformed into an error functionerf(𝑥)≈tanh2𝑥1+0.089430𝑥²√𝜋. (3.11)

Figure 3(a) shows the exact error function contrasted to (3.1), (3.2), (3.3), (3.10), and (3.11), finding a high level of accuracy for all five approximations. Besides, Figure 3(b) shows the relative error for the five approximations, where (3.10) has the lowest relative error followed by (3.1) and (3.2), while (3.11) has the highest relative error.

(a)

(b)

Figure 3 

(a) Error function erf(𝑥) and approximations. (b) Relative error for each approximation.

4. Cumulative Distribution Function

The cumulative distribution function [2] is defined asΦ(𝑥)=^𝑥_−∞𝜙(𝑡)1𝑑𝑡=2𝑥1+erf√2, (4.1) where 𝜙(𝑥) is the standard normal probability density function defined as follows:𝜙(𝑥)=1√2𝜋exp−𝑥²2. (4.2)

From (3.1) and (4.1) it can be concluded thatΦ(𝑥)≈12tanh39𝑥2√−2𝜋1112arctan35𝑥111√+12𝜋2,−∞≤𝑥≤∞. (4.3)

Then, the process applied to (3.1) is repeated to convert coefficients of (4.3) into fractions. The result is an approximate version of (4.3) now in fractions, which is given byΦ(𝑥)≈12tanh179𝑥−231112arctan37𝑥+12942,−∞≤𝑥≤∞, (4.4)where, converting the result into exponential terms and performing simplification, the expression becomesΦ(𝑥)≈exp−358𝑥23+111arctan37𝑥294+1⁻¹,−∞≤𝑥≤∞. (4.5)

Figure 4(a) shows the cumulative distribution function versus (4.3) and (4.5), it can be seen a higher level of accuracy in both approximations, except for region 𝑥<−3, where values of Φ(𝑥) are close to zero.

(a)

(b)

(c)

Figure 4 

(a) Cumulative distribution function (4.1) and approximations (4.3) and (4.5). (b) Relative error for each approximation 𝑥∈[−3,0]. (c) Relative error for each approximation 𝑥∈[0,4].

5. Table for Normal Distribution Function-Related Integrals

In [2] are shown a series of integrals without exact solution, these involve the use of the cumulative distribution function (4.1), and the standard normal probability density function (4.2). The cumulative distribution function is replaced for its approximate version (4.5) to generate the integrals in Table 1.

Table 1 

Gaussian function-related integrals.

For illustrative purposes, we have chosen values for 𝑎, 𝑏, 𝑐, 𝐶, 𝑘, and 𝑛 to be able to calculate the relative error between the exact and approximate solutions. Thus, to obtain a new solution, it can be done assigning a new value to the desired constant, or constants, and perform the required calculations. In this case, the approximate version with coefficients expressed as fractions (4.5) for the cumulative distribution function (4.1) was employed in order to make simple manual calculations easier. Nevertheless, it is also possible to use approximation (4.3).

6. Application Cases

As an example to apply the error function, two cases are considered for the unidimensional heat flow equation.

6.1. Example 1

Consider the case for a thin semi-infinite bar (𝑥⩾0) whose surface is isolated (see Figure 5); it has a constant initial temperature 𝑇_𝑜. Suddenly, zero temperature is applied at the 𝑥=0 end and is kept at that value. It is attempted to determine the temperature distribution for the bar, 𝑇(𝑥,𝑡), at any point 𝑥 and time 𝑡. This problem will be solved using techniques from Fourier integral [35].

Figure 5 

Scheme for Example  1.

From heat flow theory, it is known that 𝑇(𝑥,𝑡) should satisfy (heat conduction equation)𝑎²𝜕²𝑇𝜕𝑥²=𝜕𝑇,𝜕𝑡 (6.1) where 𝑎²=𝑘, known as thermal diffusivity. It is subject to initial condition 𝑢(𝑥,0)=𝑇_𝑜 and boundary condition 𝑢(0,𝑡)=0. The method of separation of variables is applied. This technique assumes that it is possible to obtain solutions in a product fashion𝑇(𝑥,𝑡)=𝑋(𝑥)𝑌(𝑡). (6.2)

Replacing (6.2) in (6.1) results in 𝑇(𝑥,𝑡) becoming𝑇(𝑥,𝑡)=exp−𝑘𝜆²𝑡(𝐴cos(𝜆𝑥)+𝐵sin(𝜆𝑥)), (6.3) where 𝜆 is a separation constant, 𝐴 and 𝐵 are constants to be determined.

From boundary conditions follows that 𝐴=0, and since there is no restriction for 𝜆, it is possible to integrate over 𝜆, replacing 𝐵 in function 𝐵(𝜆) such that solution to the problem takes the shape𝑇(𝑥,𝑡)=^∞₀𝐵(𝜆)exp−𝑘𝜆²𝑡sin(𝜆𝑥)𝑑𝜆, (6.4) where 𝐵(𝜆) is determined from the following integral equation:𝑇_𝑜=^∞₀𝐵(𝜆)sin(𝜆𝑥)𝑑𝜆, (6.5)

now, it is possible to express the heat distribution as follows:𝑇𝑥,𝑦=2𝑇₀𝜋^∞₀^∞₀exp−𝑘𝜆²𝑡sin(𝜆𝑉)sin(𝜆𝑥)𝑑𝜆𝑑𝑉, (6.6)

using()()=1sin𝜆𝑉sin𝜆𝑥2[cos𝜆(𝑉−𝑥)−cos𝜆(𝑉+𝑥)]_∞0,exp−𝛼𝑥²1cos𝛽𝜆𝑑𝜆=2𝜋𝛼−𝛽exp².4𝛼 (6.7)

It is possible to express (6.6) in terms of the error function𝑇(𝑥,𝑡)=2𝑇_𝑜√𝜋^√𝑥/2₀^𝑘𝑡exp−𝑤²𝑑𝑤, (6.8) or𝑇(𝑥,𝑡)=𝑇_𝑜𝑥erf2√.𝑘𝑡 (6.9)

Expressing (6.9) in terms of the approximate solution, given by (3.1), we obtain𝑇(𝑥,𝑡)≈𝑇_𝑜tanh19.5𝑥2√𝜋𝑘𝑡−55.5arctan35𝑥222√.𝜋𝑘𝑡 (6.10)

Figure 6 shows how temperature varies for a certain distance (expressed in meters) and time (measured in seconds); temperature is expressed in Kelvin. Thermal diffusivity parameter (𝑘) employed for this example is 1.6563𝐸−4m²/s which belongs to silver, pure (99.9%).

(a) Time range between 0 s and 2000 s

(b) Time range between 0 s and 250 s

(c) 𝑥=0.1m

(d) 𝑥=0.2m

(e) 𝑥=0.3m

(f) 𝑥=0.4m

(g) 𝑥=0.5m

(h) 𝑥=0.6m

Figure 6 

(a) Graphical comparison between exact (6.9) and approximate (6.10) results. (c)–(h) Relative error for the considered distance 𝑥={0.1m,0.2m,0.3m,0.4m,0.5m,0.6m}.

6.2. Example 2

Another interesting example of heat flow is the nonstationary flux in an agriculture field due to the sun (Figure 7). Suppose that the initial distribution of temperature on the field is given by 𝑇(𝑥,0)=𝑇_𝑓, and the superficial temperature 𝑇_𝑠 is always constant [36].

Figure 7 

Scheme for Example  2.

Let the origin be on the surface of the field, in such a way that the positive end for 𝑥 axis points inward the field. For symmetry reasons, it is possible to consider 𝑇 as a function of 𝑥 and time 𝑡𝑇(𝑥,𝑡). The temperature distribution is governed, again, by (6.1). It is subject to conditions 𝑇(0,𝑡)=𝑇_𝑠 and 𝑇(𝑥,0)=𝑇_𝑓. Unlike Example  1, (6.1) is expressed as an ordinary single variable second-order differential equation using the substitution𝑉(𝑢)=𝑇(𝑥,𝑡), (6.11) where()=𝑥𝑢=𝑢𝑥,𝑡2𝑎√𝑡. (6.12)

This substitution immediately converts (6.1) into𝑑²𝑉(𝑢)𝑑𝑢²(𝑢)=−2𝑢𝑑𝑉.𝑑𝑢 (6.13)

Integrating (6.13) leads to the following:𝑉(𝑢)=𝐶₁^𝑢_𝑜exp−𝑝²𝑑𝑝+𝐶₂, (6.14) where 𝐶₁ and 𝐶₂ are integration constants.

Using (2.11) and (6.12), it is possible to express (6.14) in terms of the error function by𝑇(𝑥,𝑡)=𝐶₁+𝐶₂1erf2𝑥𝑎√𝑡. (6.15)

Determining 𝐶₁ and 𝐶₂ constants results from applying the initial and boundary conditions, such that the final result is given as𝑇(𝑥,𝑡)=𝑇_𝑓−𝑇_𝑠1erf2𝑥𝑎√𝑡+𝑇_𝑠. (6.16)

The approximate solution is obtained by substituting (3.1) in (6.16) as follows:𝑇(𝑥,𝑡)≈𝑇_𝑓−𝑇_𝑠tanh19.5𝑥2𝑎√𝜋𝑡−55.535𝑥222𝑎√𝜋𝑡+𝑇_𝑠. (6.17)

Figure 8 shows how temperature varies in depth; values are in meters. The temperature is in Kelvin, and time is measured in seconds. For this case, field temperature is 𝑇_𝑓=285^∘K, surface temperature is 𝑇_𝑠=300^∘K, and the thermal conductivity value is 𝑘=𝑎²=0.003m²/s.

(a) Time range between 0 s and 1800 s

(b) 𝑥=0.01m

(c) 𝑥=0.05m

(d) 𝑥=0.1m

(e) 𝑥=0.2m

(f) 𝑥=0.3m

Figure 8 

(a) Graphical comparison between exact (6.16) and approximate (6.17) results. (b)–(f) Relative error for the considered distance 𝑥={0.01m,0.05m,0.1m,0.2m,0.3m} expressed in meters.

7. Discussion

This paper presents the normal distribution integral as a differential equation. Then, instead of using a traditional linear function 𝐿 in HPM, a nonlinear differential equation 𝐺 is used which has analytic solution (qualitatively related to the solution of the normal distribution integral). This is done to initiate the HPM method at the “closest”‌ point to the solution. In fact, in the research area of homotopy continuation methods [37–40], it is well known that when the homotopy parameter is 𝑝=0, the solution of the homotopy function must be trivial or simple to solve. Hence, it is possible to use a nonlinear differential equation 𝐺 instead of 𝐿 as long as the differential equation 𝐺 has an analytic solution (at 𝑝=0). The arctan function nested into the tanh function helps to establish a better approximation for the normal distribution integral. Then, the HPM method was successfully applied on the treatment of the normal distribution integral. From this result other related integrals were calculated, finding that the maximum relative error for the Gaussian integral was less than 180𝐸−6 (Figure 1(b)). For the case of the approximate error function (3.1) and (3.2), the relative error was less than 200𝐸−6(Figure 3(b)), and for the cumulative error function the maximum error for region 𝑥>0 was less than 90𝐸−6 (Figure 4(c)). Therefore, those approximations have a high level of accuracy, comparable to other approximations found in the literature; nevertheless, the proposed approximations in this work have such mathematical simplicity that allows to be used on practical engineering applications where the relative error does not represent a severe constraint.

The approximate solution of the cumulative distribution function was used to solve some defined and undefined integrals without known analytic solution, showing a low order relative error. Besides, the approximate error function (3.1) was employed to express, analytically, the solution for two heat flow problems, obtaining results with low order relative error.

Approximations (2.9), (3.2), and (4.4) have a low order of mathematical complexity so they become susceptible to be implemented in hardware for analog circuits. Therefore, approximations for the Gaussian, error, and cumulative functions may be part of a circuit for analog signal processing. Finally, the necessary blocks to implement some of the normal distribution functions in a circuit using the current-current mode are as follows:(1)The hyperbolic tangent was implemented in analog circuits [41].(2)Arctangent function was implemented in [42].(3)To multiply a current by a factor, it is only necessary to use current mirrors [43].(4)The addition or subtraction of constant values is achieved by connecting to the terminal a positive or negative current, respectively [43].(5)The addition of constants is equivalent to add constant current sources by using current mirrors [43].

Finally, after analyzing qualitatively the proposed approximation of this work, it is possible to establish a better approximation (see Figure 9) to the error function with hyperbolic tangents nested, resulting in the following:erf(𝑥)≈tanh77𝑥+7511625tanh147𝑥−73767tanh51𝑥,278 (7.1) where values for constants were calculated by means of the same numerical adjustment procedure employed to calculate (2.9). This approximation for the error function may be implemented in a circuit easily by means of the repeated use of the analog block for the hyperboloidal tangent [41] and some current mirrors [43].

Figure 9 

Relative error for approximation (7.1) of error function erf(𝑥).

8. Conclusion

This work presented an approximate analytic solution for the Gaussian distribution integral (by using HPM method), the error function and the cumulative normal distribution providing low order relative errors. Besides, the relative error for the error function is comparable to other approximations found in the literature and has the advantage of being a simple expression. Also, it was possible to solve, satisfactorily, a series of normal distribution-related integrals, which may have potential applications in several areas of applied sciences. It was also demonstrated that the approximate error function can be employed on the practical solution of engineering problems like the ones related to heat flow. Besides, given the simplicity of the approximations, these are susceptible of being implemented in analog circuits focusing on the analog signal processing area.