Biasing the transition probabilities in direct Monte Carlo

Reliability Engineering and System Safety 47 (1995) 59-63 ELSEVIER

0951-8320(94)00035-2

© 1995 Elsevier Science Limited Printed in Northern Ireland. All rights reserved 0951-8320/95/$9.50

Technical Note Biasing the transition probabilities in direct Monte Carlo E. Zio Centro Studi Nucleari Enrico Fermi--CESNEF, Dipartimento di Ingegneria Nucleare del Politecnico di Milano, Via Ponzio 34/3, 20133 Milano, Italy (Received 8 June 1994)

The mathematical scheme of the biased direct Monte Carlo method is presented in detail and compared with the corresponding biased analogues Monte Carlo method. With reference to the inhomogeneous Poisson processes, an example in the reliability and safety engineering field is given. This refers to sampling from a three-mode Weibull distribution by forcing each of the three modes to be exponential with decay constants 10 times larger than the natural ones. Since the natural distribution under analysis shows a rather uniform behaviour over most of the range of interest, biasing to a uniform distribution has also been investigated. Both examples are shown to lead to a simplification in the sampling procedure and an improvement in the computational efficiency.

The Monte Carlo m e t h o d o l o g y provides an efficient simulation tool which can be exploited in different fields of engineering. In this note we refer to the important field of reliability and safety engineering, which has revealed some similarities with the particle transport problem. H o w e v e r , it is worth noticing that what follows is general of the methodology. The Monte Carlo simulation m e t h o d o l o g y may be applied to reliability problems following the analogue or the direct approach. The first approach consists of sampling a transition time for the system and then selecting which c o m p o n e n t has undergone the transition. The direct approach consists of sampling the transition times for all the c o m p o n e n t s of the system and then ordering them in increasing order to produce the sequence of transitions occurring in the simulation, within the limits of the system's mission time. 1-3 If the system under study is made of a large n u m b e r of components, the direct method seems to be an inefficient procedure, since the great majority of the c o m p o n e n t s are unlikely to fail before the end of the mission time, and thus m a n y samplings would be carried out for no use. Biasing of the Monte Carlo

simulation is of p a r a m o u n t importance in reliability assessments for in most cases they deal with situations in which the c o m p o n e n t s ' transitions are rare events. In this context it appears that while the biased analog Monte Carlo is a well known and widely adopted method, 4s not m a n y applications can be found regarding the biased direct technique. This seems due to the fact that the latter approach is more intricate, less intuitive and not yet fully described, to the author's knowledge. In this note we analyze the biased direct Monte Carlo and show that this method leads to the same expressions of the biased analog case, as it should do for consistency. Let us consider a system made of n independent components, each of which m a y be in one of m a n y possible states and is characterized by a certain transition probability distribution. Let us suppose, for simplicity's sake, that all the states other than the functioning one are absorbing states. At time t = 0 the c o m p o n e n t s are assumed to be in their functioning state. Consider a generic c o m p o n e n t j and let Sj(t) denote the probability of finding the c o m p o n e n t still functioning at time t. Then Fj(t)= 1 - S j ( t ) is the 59

E. Zio

60

probability of component j failing to an absorbing state before time t, and let ~(t) be the corresponding probability density function (pdf). In the case of biasing of the transitions, the pdfs of the components are modified so as to favour those transitions which are more relevant for the quantities to be evaluated. Let Sj(t), [rift) and ~(t) be the symbols denoting the corresponding biased probabilities. Simulation of a system life-history in the biased direct Monte Carlo method begins with the sampling of the transition times tj of all the components from the biased ~(t), j = l . . . . . n. The pair (j, tj) defines the event consisting of the transition of the jth component at time ti. The sampled events are then reordered according to their time occurrences, thus leading to the sampled random walk of transitions

pair (i2, z'2)

given that the pair occurred, which is

p(i2,

r2]il,

rl)

--

(il, "(1)

has already

p(i2, r2;il, r,) p(i~, 1-1)

S.(r9 (5)

= f,.~(r2) fi Sil,(rl) I, 2

which in the biased form reads,

p(i2, Z'2 [i,, r,) = f,~(v9

k 3

L,(r9 (6)

h L,(r,)

/, :2

(i,, r,), (6, r:) . . . . . (i,,, r,,) for which r~ < r2 < • " " < r,,. The unbiased pdf for the first transition is given by

p(il, r l ) = f , ( z ' , )

I£I S,,(v,)

k=2

(1)

whereas the corresponding biased pdf from which the pair is actually sampled is

Again, in order to keep the estimator unbiased, a weight adjustment must take place in correspondence with the second transition,

p(i2, r2 l i,, r,) w2 = w(i2, r2) = w, fi(i2, r2 [ i,, r,)"

The natural extension to the general case for the weight adjustment in correspondence to the kth transition denoted by the pair (ik, r~) is

Wk=W(ik, z'k)

fi(i,, rl)=f/,(z',) lEI 3,.~(f,).

(2)

k=2

In order to retain an unbiased estimator of the quantities of interest, the initial statistical weight w, of the history must be updated as follows,

W1 = W ( i l ,

p(&, 1-,) r l ) = w o f i ( i l ' rl)

f,(1-,) h &(r,) Jill, 1J k = 2

i~

1-2;

(3)

I)

il, 1-1)=f,(1-1)f2(r2) [I Si,(1-2)

p(ik, rk l i,, r,:.., i~ ,, r~ ,) Wk ,_ . . P0k, rkJi,, r t : . . . i k i,r~ 1)

(8)

The sequences of components' transitions ends when the system reaches the mission time. As an example, let us consider the case in which the transition times of the jth component obey a Weibull distribution characterized by the parameters c~i, Ai." For this case we have

Si(t ) = e a,,,*,: ~(t) = 1 - e ~i,",; fi(t) = aiA,t'"-'e a,,,,,.

The sequence of transitions in the system life-history then continues with the second pair (i2, 1-2). The unbiased bivariate pdf that the first transition is of component il at time 1-~ and the second is of component i2 at time 1-2 is

P(i2,

(7)

(4)

Biasing of the transitions may be performed by modifying any of the parameters a i, A, j = 1. . . . . n, into the desired values 6i, "~i. j = 1. . . . . n. Correspondingly we have L(t)=e

f~(t)

L"~': ~ ( t ) = l - e

L";':

ajXjt ~i 'e L"-~'.

=

Substituting these expressions in eqns (1)-(8), the general weight update after the kth transition (ik, l-k) becomes Ogik I~i k 1-~ Ik

Wk = W(ik, l-k) = Wk-l &i, Aik r~'*

k=3

However, the pdf which plays the fundamental role in determining the two-transition-ordered sequence (il, 1-~,), (i2, 1-2) is the conditional pdf of having the

× exp - 2

m=k

(hi,,, - Mml(rk ~ '~'° ' - rk'-'l)~'

(9/

In the simple case in which a / = ~j = 1, j = 1 , . . . , n,

61

Biasing transition probabilities in direct M o n t e Carlo

the distribution reduces to the exponential; correspondingly, the general weight update after the kth transition (i~, r~) is

with m e a n AGO), for which the cdf and the pdf may be written as F(t) = 1 - e x p [ - A g ( t ) ]

Ai~ w~ = w(i~, r~) = w~_, = -

× exp - ~

(12)

f ( t ) = dg_~(t)exp[-Ag(t)]A ot

(A~., - A~)(r~ - r~ ,).

(10)

and assume g(t) = f~ ?,(r) dr, where ?,(t) has the usual meaning of transition rate, typical of the reliability problems. As a general case, we assume that ?,(t) may be written as the sum of m transition modes

m-k

Equations (9) and (10) can be easily recognized as the well known expressions pertaining to the biased analog Monte Carlo procedure, for Weibull and exponential distributions respectively.

y(t) = ~

y,(t).

(13)

i=l

Sampling random variables from an inhomogeneous univariate distribution by biasing techniques

Here, each m o d e is assumed to be represented by a t w o - p a r a m e t e r Weibull distribution with transition rate given by

Monte Carlo sampling of a r a n d o m variable t from a given univariate cumulative distribution function (cdf) F(t) may be p e r f o r m e d by solving for t the equation

yi(t) =

i = 1. . . . .

m.

(14)

By allowing several different values of the p a r a m e t e r s in eqn (14), eqn (13) m a y reasonably represent m a n y of the time-dependent transition rates of interest in the reliability assessments. The cdf F ( t ) may be shown to be a m i n i m u m extreme value distribution, ~° in which the parent distributions are the individual Weibull distributions

(11)

= F(t)

otiAi t~'-I

where ~: ~ U(O, 1) is a uniform r a n d o m n u m b e r in the interval (0,1). For simple distributions, such as uniform and exponential, this general inversion technique is quite straightforward. However, in m o r e complex cases the solution of eqn (11) might b e c o m e a difficult task since one has either to solve a transcendent equation or to use ad hoc M o n t e Carlo methods or specific algorithms. Let us consider an inhomogeneous Poisson process

F~(t) = 1 - e x p [ - A : m ] (15) f ( t ) = a i A : m-~ e x p [ - A : " ' ]

Therefore, in order to obtain the values of the r a n d o m

Direct Monte Carlo of a three mode Weibull distribution

xl0 -5

cpu time = 0.1658 s n. trials = 100000

4

~ ; " '

~, ~.~ ~ ~ ' ' "~' ', i

1

I

,

~fl w

I~1 i

., m I

i

iq

ill

II'~l~ i

Ill

l ~

q

rl

l s¢

I

r '

,i

.

i i

i

i

i

j

I

,

t (~

~Ii _tJ

)

d

0

100

200

,

,

300

400

,

,

500

,

I it

600

time Fig. 1. Comparison of the direct Monte Carlo curve (dashed line) to the analytical (solid line) three-mode Weibull distribution.

E. Z i o

62

Biased Monte Carlo of a three mode Weibull distribution

x10-5 7

cpu time = 0.0136 s n. trials = 10000

4

,

/* .II

I

0

0

L

I

I

L

I

100

200

300

400

500

600

time Fig. 2. Comparison of the Monte Carlo curve (dashed line) to the analytical (solid line) three-mode Weibull distribution, in the case for which the distribution is forced to be exponential.

variable t, we m a y s a m p l e from each of the F,(t) by using the g e n e r a l i n v e r s i o n t e c h n i q u e , a n d t h e n take the m i n i m u m . T h e g e n e r a l s c h e m e p r e s e n t e d in the p r e v i o u s s e c t i o n for the biasing of the t r a n s i t i o n s in direct

M o n t e C a r l o m a y be e x p l o i t e d here for the s a m p l i n g of r a n d o m variables. I n d e e d , the d i s t r i b u t i o n from which we sample m a y be forced to be o n e for which a s o l u t i o n to e q n (11) m a y be easily o b t a i n e d . Let us, for e x a m p l e , m o d i f y the p a r a m e t e r s of the W e i b u l l

Biased Monte Carlo of a three mode WeibuU distribution

x10-5 7

cpu time = 0.0104 s n. trials = 10000

r

f

~

-. ,

I

t)

0

-, . . . . . . . . . .

,

,-,'

,

,

-_

i

i

L

i

I

100

200

300

400

500

.,,.,p,;?,,

600

time Fig. 3. Comparison of the Monte Carlo curve (dashed line) to the analytical (solid line) three-mode Weibull distribution, in the case for which the distribution is forced to be uniform over the range (0, 600).

Biasing transition probabilities in direct Monte Carlo distribution modes exponential:

so

as

to

force

them

to

be

63

a figure of merit of the form 1

hi --> Xi

We then have, from eqn (14), that the biased transition rates of the different modes are all constant •/,(t) = ,~,.

(17)

By so doing, the biased distribution becomes P ( t ) = 1 - exp ( - ~ )~,t) i=1

o.2T

where T is the time required for a single Monte Carlo sampling and 0 -2 is the average of the squared deviations between the Monte Carlo and the analytical values. In the three cases here investigated we found q = 2.75 x 1016 for the direct Monte Carlo, q = 12-43 x 1016 for the biased-to-exponential Monte Carlo and q = 31.35 x 10 ~v for the biased-to-uniform case.

(18)

f ( t ) = ~ Ai e x p ( - ~ ftit) i:l

q

(16)

ai---~ ~i = 1.

ACKNOWLEDGMENTS

i=1

thus greatly simplifying the Monte Carlo sampling (indeed, the inversion of eqn (11) is now straightforward). The weight associated with each individual Monte Carlo sampling, initialized to unity, must be appropriately updated by multiplying each by the ratio of the natural pdf of eqn (12) to the biased one of eqn (18), as shown in eqn (8). As a simple example, we have considered a three-mode Weibull distribution, of p a r a m e t e r s A~ = 4 X 10 ~:

A2 =

2 x 10 ~;

al : 0 " 7 ;

O~2 : 0"6;

A~ = 1 X 10 s: O~3 : 0"5.

Figure 1 shows the results of l0 t direct Monte Carlo samplings c o m p a r e d to the analytical distribution, in the range of (0, 600). We then biased each of the three modes of the Weibull distribution to be exponential (a, = 1, i = 1, 2, 3), with the A~ values being 10 times the natural ones. Figure 2 shows the results of this biased Monte Carlo with 104 trials: the i m p r o v e m e n t in the efficiency of the sampling is evident and also the c o m p u t e r time required is reduced. Finally, we biased the natural distribution to a uniform distribution over the interval (0,600) and the results for 10 4 Monte Carlo samplings are reported in Fig. 3. As expected, in this case the efficiency is further improved, since the natural distribution does not differ much from a uniform distribution over the range of interest. The efficiencies of the direct, the biased-to-exponential and the biased-to-uniform methods may be summarized in

The author wishes to thank Professor Marzio Marseguerra for his contribution to this work.

REFERENCES 1. Goertzel, G. & Kalos, M. H., Monte Carlo methods in transport problems. Progress in Nuclear Energy, Vol. II, Series II. Pergamon Press, New York, 1958. 2. Rubinstein, R. Y., Simulation and the Monte Carlo Method, John Wiley & Sons, New York, 1981. 3. Siu, N. O. & Deoss, D. L., A simulation model for dynamic system availability analysis, MITNE-287, Massachussetts Institute of Technology, Cambridge (1989). 4. Lewis, E. E. & Bohm, F., Monte Carlo simulation of Markov unreliability models. Nucl. Engng Design, 77, (1984) 49-62. 5. Zhuguo, T. & Lewis, E. E., Component dependency models in Markov Monte Carlo simulation. Reliab. Engng System Safety, 13 (1985) 45-61. 6. Lewis, E. E. & Zhuguo, T., Monte Carlo reliability modeling by inhomogeneous Markov processes. Reliab. Engng System Safety, 16 (1986) 277-96. 7. Wu, Y. F. & Lewins, J. D., Monte Carlo studies of engineering system reliability. Annals Nucl. Energy, 19 (10-12) (1992) 825-59. 8. Marseguerra, M. & Zio, E., Nonlinear Monte Carlo reliability analysis with biasing towards top event. Reliab. Engng System Safety, 40 (1993) 31-42. 9. Kendall, M., Stuart, A. & Ord, J. K., Kendall's Advanced Theory of Statistics, Vol. I, Fifth edition, Charles Griffin & Co., London, 1987. 10. Gumbel, J. E., Statistics of Extremes, Columbia University Press, New York, 1958.