Referring to Fig. 6.6, a rectangular region Ш = {(x, y)a < x < b, 0 < y < c} is superimposed to enclose the area Ф = {(x, y)a < x < b, 0 < y = g(x) < c} represented by Eq. (6.49). By the hit-and-miss method, the rectangular region Ш containing the area under g(x), that is, Ф, is hung on the wall, and one is to throw n darts on it. Assume that the darts fly in a random fashion and that all n darts hit within the rectangular region. The area under g(x), then, can be estimated as the proportion of n darts hitting the target multiplied by the known area of rectangular region Ш, that is,

G = A(6.50)

where G is the estimate of the true area G under g(x), A = c(b – a) is the area of the rectangular region, and nh is the number of darts hitting the target out of a total of n trials.

The hit-and-miss method can be implemented numerically on a computer. The two coordinates (Xi, Yi) on the rectangular region Ш, which represents the location where the ith dart lands, are treated as two independent random variables that can be generated from two uniform distributions. That is, Xi is generated from U(a, b) and Yi from U(0, c). When Yi < g(Xi), the dart hits its target; otherwise, the dart misses the target. A simple hit-and-miss algorithm is given as follows:

1. Generate 2n uniform random variates from U(0, 1). Form them arbitrarily into n pairs, that is, (u1, u[), (u2, u’2),…, (un, u’n).

2. Compute xi = a + (b – a)ui and g(xi), for i = 1,2,…, n.

3. Count the number of cases nh that g(xt) > cu.

4. Estimate the integral G by Eq. (6.50).

Note that G is an estimator of the integral G; it is therefore also a random variable. It can be shown that G is unbiased, namely,

E(G) = A x = Ap = A (A = G (6.51)

where nh/n, the proportion of n darts hitting the target, is an unbiased estimator of the true probability of hits, and p simply is the ratio of the area under g(x) to the area of the rectangular region. Furthermore, the standard error associated with the estimator G is

As can be seen, the precision associated with G, represented by its inverse of standard deviation, using the hit-and-miss Monte Carlo integration method increases with n1/2.

A practical question is how many trials have to be carried out so that the estimated G satisfies a specified accuracy requirement. In other words, one would like to determine a minimum number of trials n such that the following relationship holds:

P(|G – G<e) > a (6.53)

in which є is the specified maximum error between G and G, and a is the mini­mum probability that G would be within є around the exact solution. Applying the Chebyshev inequality, the minimum number of trials required to achieve Eq. (6.53) can be determined as (Rubinstein, 1981)

(1 – p)p[c(b – a)]2 (1 – p)pA2

n >———— / 2——— = n * 2 (6.54)

(1 – а)є2 (1 – а)є2

Note that the required number of trials n increases as the specified error level є decreases and as the confidence level a increases. In addition, for the specified є and a, Eq. (6.54) indicates that the required number of trials n can be reduced by letting p approach 1. This implies that selecting an enclosed region Ш as close to Ф as possible would reduce the required number of trials. However, consideration must be given to the ease of generating random variates for U’ in the algorithm.

When the number of trials n is sufficiently large, the random variable T,

G — G

T = ^^ (6.55)

sG

approximately, has the standard normal distribution, that is, T ~ N(0, 1), with sG being the sample estimator of aG, that is,

Hence the (1 – 2a)-percent (a < 0.5) confidence interval for G then can be ob­tained as

G ± sgza

with za = Ф 1(1 – a).

Example 6.6 Suppose that the time to failure of a pump in a water distribution system follows an exponential distribution with the parameter в = 0.0008/h (i. e., 7 failures per year). The PDF of the time to failure of the pump can be expressed as

ft(t) = 0.0008e-0 0008t for t > 0

Determine the failure probability of the pump within its first 200 hours of operation by the hit-and-miss algorithm with n = 2000. Also compute the standard deviation associated with the estimated failure probability and derive the 95 percent confidence interval containing the exact failure probability.

Solution The probability that the pump would fail within 200 hours can be computed as

r-200 ,■ 200

ft (t) dt = 0.0008e—00008t dt

00 which is the area under the PDF between 0 and 200 hours (Fig. 6.7). Using the hit – and-miss Monte Carlo method, a rectangular area with a height of 0.0010 over the interval [0, 200] is imposed to contain the area representing the pump failure proba­bility.

The area of the rectangle can be easily determined as A = 0.001(200) = 0.2. The hit-and-miss algorithm then can be outlined in the following steps:

1. Initialize i = 0 and щ = 0.

2. Let i = i + 1, and generate a pair of standard uniform random variates (ui, ui) from U(0, 1).

3. Let ti = 200ui, and compute ft(ti) = 0.0008e—00008ti, y = 0.001ui.

4. If ft(ti) > y, nh = nh + 1. If i = 2000, go to step 5; otherwise, go to step 1.

5. Estimate the pump failure probability as pf = A(nh/n) = 0.2(nh/n).

Using the preceding algorithm, 2000 simulations were made, and the estimated pump failure probability is pif = 0.2(nh/n) = 0.2(1500/2000) = 0.15. Comparing with the exact failure probability pf = 1 — exp(-0.16) = 0.147856, the estimated

failure probability by the hit-and-miss method, with n = 2000 and the rectangular area chosen, has a 1.45 percent error relative to the exact solution.

The associated standard error can be computed according to Eq. (6.56) as

Assuming normality for the estimated pump failure probability, the 95 percent confi­dence interval containing the exact failure probability pf is

pf ± Z0.975SPf = (0.1462,0.1538)

where zq.975 = 1.96.

In reliability analysis, computations of system and/or component reliability and other related quantities, such as mean time to failure, essentially involve inte­gration operations. A simple example is the time-to-failure analysis in which the reliability of a system within a time interval (0, t) is obtained from

where ft (t) is the failure density function. A more complex example of the reli­ability computation is by load-resistance interference in that the reliability is

Ps = P[R(Xr) > L(Xl)] = P [W(Xr, Xl) > 0] = P [W(X) > 0]

= f fx (x) dx (6.48)

JW (x)>0

where R(XR) and L(XL) are, respectively, resistance and load functions, which are dependent on some basic stochastic variables XR = (X1, X2,…, Xm) and XL = (Xm+1,Xm+2,…,XK), and W(X) is the performance function. As can be seen, computation of reliability by Eq. (6.48) involves K-dimensional integrations.

For cases of integration in one or two dimensions, such as Eq. (6.47), where the integrands are well behaved (e. g., no discontinuity), conventional numer­ical integration methods, such as the trapezoidal approximation or Simpson’s rule (see Appendix 4A), are efficient and accurate. For example, using Simp­son’s rule, the error in a one-dimensional integration is O(n-4), with n being the number of discretizations, and the error in a two-dimensional integration is O(n-2). Gould and Tobochnik (1988) show that, in general, if the error for the one-dimensional integration is O(n-a), the error with a K-dimensional in­tegration would be O(n~a/K). As can be seen, the accuracy of conventional nu­merical integration schemes decreases rapidly as the dimension of integration increases. For multiple integrals, such as Eq. (6.48), the Monte Carlo method becomes a more suitable numerical technique for integration.

To illustrate the basic idea of the Monte Carlo integration, consider a simple one-dimensional integration