Application of the interacting particle system method to piecewise deterministic Markov processes used in reliability
AApplication of the interacting particle system method to piecewisedeterministic Markov processes used in reliability
H. CHRAIBI a , A. DUTFOY a , T. GALTIER a,b, ∗ , J. GARNIER c a EDF R&D - D´epartement PERICLES, 7 Boulevard Gaspard Monge, 91120 Palaiseau, France b Universit´e Paris-Diderot - Laboratoire de Probabilit´es Statistique et Mod´elisation, 75205 Paris Cedex 13,France c Ecole Polytechnique - Centre de Math´ematiques Appliqu´ees, 91128 Palaiseau Cedex, France
Abstract
Variance reduction methods are often needed for the reliability assessment of complex in-dustrial systems, we focus on one variance reduction method in a given context, that is theinteracting particle system method (IPS) used on piecewise deterministic Markov processes(PDMP) for reliability assessment . The PDMPs are a very large class of processes whichbenefit from high modeling capacities, they can model almost any Markovian phenomenonthat does not include diffusion. In reliability assessment, the PDMPs modeling industrialsystems generally involve low jump rates and jump kernels favoring one safe arrival, we callsuch model a ”concentrated PDMP”.Used on such concentrated PDMPs, the IPS is inefficient and does not always providea variance reduction. Indeed, the efficiency of the IPS method relies on simulating manydifferent trajectories during its propagation steps, but unfortunately concentrated PDMPsare likely to generate the same deterministic trajectories over and over. We propose anadaptation of the IPS method called IPS+M that reduces this phenomenon. The IPS+Mconsists in modifying the propagation steps of the IPS, by conditioning the propagation toavoid generating the same trajectories multiple times. We prove that, compared to the IPS,the IPS+M method always provides an estimator with a lower variance. We also carry outa quick simulation study on a two-components system that confirms this result.
Keywords:
Rare event – Reliability assessment – PDMP or PDP – Dynamical hybridsystem – PycATSHOO – Sequential Monte-Carlo Samplers – Feymann-Kac particle filters ∗ Corresponding author
Email address: [email protected] ( T. GALTIER ) Preprint submitted to Chaos: An Interdisciplinary Journal of Nonlinear Science May 23, 2019 a r X i v : . [ s t a t . C O ] M a y . Introduction For both safety and regulation issues, the reliability of industrial systems has to be as-sessed. The considered systems (dams or nuclear power plants, for instance) are complexdynamic hybrid systems, so only simulation methods can be reasonably considered to as-sess their reliability. The failure of such a dynamic hybrid system generally correspondsto a physical variable of the system (temperature, pressure, water level) entering a criticalregion. The simulation of such a system requires to accurately model the trajectory of thephysical variables. The evolution of these physical variables are generally determined bysimple differential equations derived from the laws of physics. As the physical context giv-ing the differential equations generally depends on the statuses of the multiple componentsof the systems (on, off, or failed), the differential equations can change whenever there isa change in the components statuses. To encounter for this hybrid interplay between thediscrete process of components’ statuses and the continuous evolution of the physical vari-ables, we model the evolution of the state of a system by a piecewise deterministic Markovprocess (PDMP) [5, 6, 7, 3]. PDMPs are meant to represent a large class of Markovianprocesses that do not include diffusion, as such, they benefit from high modeling capacity:they can model many industrial systems. For instance, EDF has recently developed thePyCATSHOO toolbox [2], which allows the modeling of dynamic hybrid systems, the mainoption within this toolbox is to evaluate the dependability criteria of the studied systemsby Monte Carlo simulations.As the industrial systems are highly reliable, their failure is a rare event, and the naiveMonte-Carlo simulation method (MC) is computationally too intensive in this context. Theobjective of our work is to set up new algorithms to accelerate the reliability assessment ofsuch industrial systems. To do so we want to use faster methods, such as variance reductionmethods. A variance reduction method is a method that yields an unbiased estimator witha smaller variance than the Monte-Carlo estimator. The estimation being more accurate,we need less simulation runs to reach the desired accuracy, thus we save computationaltime. The variance reduction is generally achieved by altering the simulation process, andthen correcting the bias induced by this modification of simulation process by appropri-ately weighting each simulation. Several acceleration methods for variance reduction canbe proposed, such as importance sampling methods or particle filter methods (also calledsubset simulations). In this article we focus on one example of particle filter method: theinteracting particle system (IPS) method [9].Unlike in importance sampling methods, in the IPS method we keep simulating accordingto the original system. Another difference is that we do not simulate directly the trajec-tories on the entire observation duration, the observation duration is subdivided into smallintervals of time and we simulate the trajectories sequentially one interval of time after theother one. These sequential simulations of the trajectories consist alternating between anexploration step and a selection step. During the exploration step we simulate the trajecto-ries on a small time interval, therefore exploring the most probable trajectories on a short2orizon of time. Then we apply a selection step on these trajectories replicating the trajec-tories which seem ”close” to failure and giving up the less ”promising” trajectories. At thenext exploration step, only replicated trajectories are continued, before the next selection,and the next exploration and so on... until each selected trajectory reaches the end of theobservation time window. This way the effort of simulation is concentrated on selected tra-jectories, which have higher chance of becoming a failing trajectory before the end of theobservation time window. In the end we get more failing trajectories to fuel our estimation,and, if the selection was well done, the IPS method yields an unbiased estimator with asmaller variance than the MC estimator.When we try to apply the IPS method in order to estimate the failure probability ofsuch reliable complex hybrid systems, the IPS methods turn out to be inefficient. Indeedthe estimation provided by the IPS method often has, in this case, a higher variance thanthe one of the Monte-Carlo estimator. IPS does not perform well in this context because theapplication case (a reliable complex hybrid system) makes it hard to conduct the explorationsteps of the IPS method efficiently. Indeed, an industrial system is often modeled by whatwe call a ”concentrated PDMP”, which is a PDMP with low jump rates and concentratedjump kernels on boundaries. The typical jump rate is low because it is the sum of thefailure rates of the components in working condition and of the repair rates of the failedcomponents which are all very low. These failure and repair rates are very low because thecomponents are reliable and their repairs are slow. The typical jump kernel is concentratedbecause most of the probability mass of a jump is concentrated on one (safe) output. Indeedthe jump kernel on boundaries model the automatic control mechanisms within the system.During such a control mechanism there is a small probability that some component(s) fail ondemand but the most likely output is that the system jumps on the safe state aimed by thecontrol mechanism, consequently the probability mass of the jump kernel is concentratedon this output. Due to these characteristics of the model there is a high probability that nocomponent failure or repair occurs during the short exploration time, and with a PDMP itmeans that all the trajectories are likely to follow the same deterministic paths. So when weexplore the trajectory space most simulated trajectories end up being the same one, hencelimiting our exploration of the trajectory space. To avoid this pitfall, a particular filterwas proposed in [20] that enhances the occurrence of random jumps (failure or repairs) ormodifies the occurrence time of the last jump. However the proposed method is limited toa different case of PDMP. This class of PDMP does not include concentrated PDMP, as itcontains only PDMP without boundaries which allows continuous jumps kernels, i.e. it doesnot allow to model automatic control mechanisms in components. Moreover, an a prioribound on the number of jumps in a time interval is required, but in our case we do not havesuch information.In order to adapt the IPS to concentrated PDMPS, we propose instead to use an approachbased on the memorization method developed in [16]. The idea is to start the exploration byfinding the most likely trajectories continuing each batch of replicated trajectories. Then, wecondition the rest of the exploration to avoid these trajectories. As a result, the simulatedtrajectories have much more chances to differ which improves the quality of the exploration3nd reduce the variance of the estimator. To correct the bias induced by this modificationof the simulation process, we have to modify the weight of each trajectory. We call ouradaptation of the IPS to PDMP the IPS+M for sequential Monte-Carlo sampler with mem-orization method.The rest of the paper paper is organized as follows: Section 2 is dedicated to the presen-tation of our model of the system, Section 3 presents the IPS method and introduces theoptimal potential functions for this method, Section 5 presents the IPS+M method whichadapts the IPS algorithm for PDMPs with low jump rates, Section 6 explains how to forcethe differentiation of the trajectories using the memorization method, and finally in Section7 we illustrate the better efficiency of the IPS+M method on a toy example.
2. A model of the system based on a PDMP with discrete jump kernel
We denote by Z t the state of the system at time t . Z t is the combination of the physicalvariables of the system, noted X t , and of the statuses of all the components within thesystem, noted M t : Z t = ( X t , M t ). We consider that X t ∈ R d , and that M t ∈ M = { On, Of f, F } N c where F corresponds to a failed status, and N c is the number of componentsin the system. The value of M t is sometimes referred as the mode of the system. Here weconsider only three categories for the status of a component : On, Of f, F , but it is possibleto include more categories as long as the set of the possible modes M stays countable.The process Z t is piecewise continuous, and each discontinuity is called a jump. Betweentwo jumps there is no change in the components’ statuses, and the dynamics of the physicalvariables can be expressed thanks to an ordinary differential equation derived from the lawof physics: d X t dt = F M t ( X t ) . We note φ ( x,m ) ( t ) the solution of this equation when X = x and M = m . Then for anytime s > , if T is the time separating s from the next jump time, we have ∀ t ∈ [0 , T ) , Z s + t = (cid:0) X s + t , M s (cid:1) = (cid:0) φ ( X s ,M s ) ( t ) , M s (cid:1) . Similarly, a flow function on the states can be defined. If z = ( x, m ), then we defineΦ z ( t ) = (cid:0) φ ( x,m ) ( t ) , m (cid:1) , and so ∀ t ∈ [0 , T ) , Z t + s = Φ Z s ( t ) . (1)As the physical variables are often continuous, the jumps are essentially used to modelchanges in the statuses of the components. These jumps can occur for two reasons.Firstly, a jump can correspond to an automatic control mechanism (See Figure 1). In agiven mode m , such mechanism is typically triggered when the physical variables crosssome threshold. For each mode m , we define an open and connected set Ω m , so that thesethresholds determine its boundary ∂ Ω m . So when M t = m the value of the physical variables4 igure 1: A jump at boundary. X t is restricted to the set Ω m ⊂ R d . Letting E m = { ( x, m ) , x ∈ Ω m } be the set of the possiblestates with mode m , in terms of state a jump associated to a control mechanism is triggeredwhenever the state Z s + t hits the boundary of E m . The set of possible states is thereforedefined by: E = (cid:91) m ∈ M E m . Secondly, jumps can correspond to a spontaneous failure or repair. In such case the jumpoccurs before Z s + t hits the boundary of E m . The occurrence time of a spontaneous jump is Figure 2: A spontaneous jump. modeled by a jump rate λ ( Z s + t ). For instance, when we consider that components fail orare repaired one at a time, this jump rate is the sum of the failure rates of the components5n working condition and of the repair rates of the broken components. In a more generalcase, each transition from a mode m to a mode m + has its own rate λ m → m + ( X s + t ), and wehave λ ( Z s + t ) = (cid:88) m + ∈ M λ M s → m + ( X s + t ) (2)We define the cumulative jump rate Λ z ( t ) by Λ z ( t ) = (cid:82) t λ (cid:0) Φ z ( u ) (cid:1) du , so that ∀ t ∈ [0 , T ) , Λ Z s ( t ) = (cid:82) t λ (cid:0) Z s + u (cid:1) du . Eventually the cumulative distribution function (cdf) of T (the time untilthe next jump starting form a state Z s = z ) takes the form: P ( T ≤ t | Z s = z ) = (cid:26) − exp [ − Λ z ( t )] if t < t ∗ z , t ≥ t ∗ z . (3)Here t ∗ z = inf { t > , Φ z ( t ) ∈ ∂E m } is the time until the flow hits the boundary starting froma state z = ( x, m ). When there is no boundary i.e. { s > , Φ z ( s ) / ∈ E m } = ∅ , we take theconvention that t ∗ z = + ∞ . With this definition, the law of T has a continuous part associatedto spontaneous jumps (spontaneous failures and repairs) and a discrete part associated toforced jumps (control mechanisms). A possible reference measure for T knowing Z s = z isthen ∀ B ∈ B ( R + ) , µ z ( B ) = leb ( B ∩ (0 , t ∗ z )) + t ∗ z < ∞ δ t ∗ z ( B ) , (4)where leb ( . ) corresponds to the Lebesgue measure.If a jump occurs at time S , then the distribution of the destination of the jump is expressedby a transition kernel K Z - S where Z - S is the departure state of the jump. If E is the closureof E , we have Z - S ∈ E . Let Z + S ∈ E be the arrival state, and B ( E ) be the Borelian σ -algebraon E . Then, the law of a jump from a departure state z - is defined by: ∀ B ∈ B ( E ) , P (cid:0) Z + T ∈ B | Z - T = z - (cid:1) = K z - ( B ) . (5)For any departure state z - ∈ E , we define x a ( z - , m + ) as the arrival for the vector of phys-ical variables when the transition m - → m + is triggered from a state z - . We define by z a ( z - , m + ) = (cid:0) x a ( z - , m + ) , m + (cid:1) the arrival state when the transition m - → m + is triggeredfrom a state z - . Using this definition, for any departure state z - ∈ E , we define the measure ν z - on (cid:0) E, B ( E ) (cid:1) by: ∀ B ∈ B ( E ) , ν z - ( B ) = (cid:88) m + ∈ M δ z a ( z - ,m + ) ( B ) . (6)The kernel K z - is absolutely continuous with respect to the measure ν z - . Denoting by K z - its density with respect to ν z - we have: P (cid:0) Z + T ∈ B | Z - T = z - (cid:1) = (cid:90) B K z - ( z + ) dν z - ( z + ) . (7)6ote that we only consider systems for which ν z - is discrete, whatever the departure state z - ∈ E may be. This hypothesis of discreteness of ν z - is mandatory to apply the methodpresented in Section 5.If ν z has a Dirac contribution at point z , then the kernel must satisfy K z ( z ) = 0 , so itis not possible to jump on the departure state. In some applications the physical variablesare continuous, so that K z - ( z ) is zero whenever x − (cid:54) = x . In some cases, one might want themodel to include renewable and aging components. Then the vector X t should include thetime since the last renewal of such component, and so the vector X t can be discontinuous atthe time of a renewal. In both situations, when z - is not on the boundary of E m - , the jumpkernel has this form: ∀ z - ∈ E, B ∈ B ( E ) , K z - ( B ) = (cid:88) m + ∈ M λ m - → m + ( x - ) λ ( z - ) δ z a ( z - ,m + ) ( B ) , (8)and K z - ( z + ) = λ m - → m + ( x - ) λ ( z - ) x + = x a ( z - ,m + ) . (9)When z - ∈ ∂E m - , a control mechanism is triggered : some components are required to turnon, or to turn off, so that the system reaches a desired state z c . The transition to this stateis usually very likely, so when z - ∈ ∂E m - , the jump kernel K z - tends to concentrate all ora big portion of the probability mass on this state z c . In industrial systems there can becomponents which have a small probability to fail when they are asked to turn on. Thisphenomenon is referred as a failure on demand. We denote by κ i ( z − ) the probability thatthe i th component fails on demand when the control is triggered. We denote by f od ( i, m - , m )the indicator being equal to one if the i th component fails on demand during the transitionfrom m - to m , and to zero otherwise. We also define ask ( z - , z + ) as the set gathering theindices of the components supposed to turn on during the control mechanism triggered from z - to z + . Finally, when z - ∈ ∂E m - , we have that: K z - ( z + ) = x + = x a ( z - ,m + ) (cid:89) i ∈ ask ( z - ,z + ) (cid:0) κ i ( z − ) (cid:1) fod ( i,m - ,m + ) (cid:0) − κ i ( z − ) (cid:1) (1 − fod ( i,m - ,m + )) . (10)To generate Z t = ( Z s ) s ∈ [0 ,t ] a trajectory of the states of the system, one can repeat thefollowing steps: Starting with s = 0,1. Given a starting state Z s = z , generate T the time until the next jump using (3),2. Follow the flow Φ until s + T using (1), and set the departure state of the next jumpas being Z - T = Φ z ( T ),3. generate Z s + T the arrival state of the next jump using K Z - T
4. repeat starting with s = s + T until one gets a trajectory of size t .Defined in this way, the process Z t is Markovian [4].7 .2. Law of trajectory Let us denote by n ( Z t ) the number of jumps in the trajectory Z t , and by S k the timeof the k th jump in the trajectory Z t (with the convention that S = 0 and S n ( Z t )+1 = t ). T k = S k +1 − S k denotes the time between two consecutive jumps. We define the σ -algebra S t as the σ -algebra generated by the sets in (cid:83) n ∈ N ∗ B (cid:16)(cid:110)(cid:0) z k , t k (cid:1) k ≤ n ∈ ( E × R ∗ + ) n , n (cid:80) i =0 t i = t (cid:111)(cid:17) ,where B ( . ) indicates the Borelians of a set. Letting Θ t : Z t → (( Z k , T k )) ≤ k ≤ n ( Z t ) be theapplication giving the skeleton of a trajectory, we can define the law of trajectories as animage law through Θ t . We have that for B ∈ S t : P z o (cid:16) Z t ∈ Θ − t ( B ) (cid:17) = (cid:90) B n (cid:89) k =0 (cid:16) λ z k ( t k ) (cid:17) tk
0, are called thepreponderant trajectories. For example, assuming that λ is positive, if the trajectory z t involves no jump, it is preponderant because we have P (cid:0) Z t = z t (cid:1) = 1 − exp (cid:0) − Λ z ( t ) (cid:1) > ν z - are discrete, there can be an other typeof preponderant trajectories, which are the trajectories that only jump on the boundaries ∂E m . Indeed, such trajectories z t would verify P (cid:0) Z t = z t (cid:1) = n (cid:89) k =0 exp (cid:104) − Λ z sk ( t ∗ z sk ) (cid:105) n (cid:89) k =1 K z − sk ( z s k ) > , (12)where s k is the time of the k th jump in z t . Conversely, some realizations can be consideredas negligible as they verify P (cid:0) Z t = z t (cid:1) = 0. These trajectories are the ones that involve aspontaneous jump, i.e. a jump starting from the interior of a set E m .This can be better understood by looking at the equation (11) which shows that the measure ζ z ,t defined by : ∀ B ∈ S t , ζ z ,t (Θ − ( B )) = (cid:90) dδ t ∗ n ( t n ) dν z − n ( z n ) ( z k ,t k ) k ≤ n ∈ B dµ z n − ( t n − ) ... dν z − ( z ) dµ z o ( t ) (13)where z − j = Φ z j − ( t j − ), and t ∗ n = t − (cid:80) n − i =0 t i is a reference measure for the law of trajecto-ries. This reference measure highlights the fact that the trajectories with no jump or with8nly jumps on boundaries can concentrate some probability mass on their own, becausethey refer to a Dirac contribution of the measure ζ z ,t . Indeed for such trajectories the timesbetween two consecutive jumps verify t k = t ∗ z sk ∀ k < n , and so, they always refer to thediscrete part of the measures µ z k − ∀ k < n , therefore such trajectories are related to thediscrete part of ζ z ,t . Equation (13) also shows that the remaining probability mass is dis-tributed continuously among the trajectories with at least one spontaneous jump. Indeed,in a negligible trajectory, if for example the k + 1 th jump is a spontaneous jump, then t k relates to the continuous part of the reference measure µ z k − . For reliability assessment of a highly reliable system one often models the system by aPDMP with low jump rates and concentrated jump kernels on the boundaries. Indeed, thecomponents of the system are often reliable and their repair takes time, hence the low jumprates, and as failures on demand during a control mechanism are unlikely the jump kernelson boundaries are concentrated on one safe arrival state (i.e. the state aimed by the controlmechanism). We call this kind of PDMP a concentrated PDMP, mainly because the law ofone trajectory concentrates a big part of its probability mass.This trajectory happens to be the trajectory with no failure and no repair. As jumprates are low the probability of not having a spontaneous jump is close to one. For instanceat a k + 1-th jump this probability verifies P z sk (cid:16) T k = t ∗ z sk (cid:17) = exp (cid:104) − Λ z sk ( t ∗ z sk ) (cid:105) (cid:39) . So only jumps on boundaries are likely, and when the process hits a boundary ∂E m , thearrival state aimed by the control mechanism is very likely. Denoting by z s k this state for a k -th jump we have : K z − sk ( z s k ) (cid:39) . So if z t is a trajectory with no failure and no repair of reasonable size we have: P (cid:0) Z t = z t (cid:1) = n (cid:89) k =0 exp (cid:104) − Λ z sk ( t ∗ z sk ) (cid:105) n (cid:89) k =1 K z − sk ( z s k ) (cid:39) . (14)In this article we adapt the IPS method for concentrated PDMPs. It is interesting to makethe connection between our work and the modified particle filter developed in [12] whichcan be applied to the related context of discrete time processes with discrete spaces. Let t f be an observation time. We denote by E t f the set of trajectories of size t f thatverify (1), and h ∈ M ( E t f ). The methods presented in this article can be used for theestimation of any quantity p h defined by: p h = E [ h ( Z t f )] .
9e are interested in estimating the probability, noted p D , that the system fails before thefinal observation time t f knowing it was initiated in a state z . Letting D be the set oftrajectories of length t f which pass through the critical region D ⊂ E , we have p D = P z (cid:0) Z t f ∈ D (cid:1) = E z (cid:2) D ( Z t f ) (cid:3) . Although our application relates to the case where h = D , the IPS method will be presentedwith an arbitrary (bounded) function h .
3. The IPS method
For a measure η and a bounded measurable function f we note η ( f ) = (cid:82) f dη . For a mea-sure η and a kernel V , ηV denotes the measure such that ηV ( h ) = (cid:82) (cid:82) h ( y ) V ( dy | x ) dη ( x ),and for a bounded measurable function f , V ( f ) is the function such that V ( f )( x ) = (cid:82) h ( y ) V ( dy | x ).Consider a subdivision of the interval [0 , t f ] into n sub-intervals of equal lengths, noted[ τ k , τ k +1 ), and such that 0 = τ < τ < · · · < τ n -1 < τ n = t f . Let V k be the Markoviantransition measure extending a trajectory of size τ k into a trajectory of size τ k +1 , such that ∀ z τ k ∈ E τ k , B ∈ E τ k +1 , P (cid:0) Z τ k +1 ∈ B | Z τ k = z τ k (cid:1) = V k ( B | z τ k ) . (15)For each k < n we denote G k the potential function on E τ k , such that: ∀ z k ∈ E k , G k ( z k ) ≥ . (16)These potential functions are the main inputs of the method. The choice of the potentialfunctions is important, because it will ultimately determine the variance of the estimator of p h provided by the IPS method. A good choice depends on the system and on the targetfunction h . The potential functions are used to define the target probability measures ˜ η k foreach k ≤ n , such that: ˜ η k ( d z τ k ) ∝ k (cid:89) s =0 G s ( z τ s ) k − (cid:89) s =0 V s ( d z τ s +1 | z τ s ) , (17)or equivalently ∀ B ∈ E τ k , ˜ η k ( B ) = E (cid:104) B ( Z τ k ) (cid:81) ks =0 G s ( Z τ s ) (cid:105) E (cid:104)(cid:81) ks =0 G s ( Z τ s ) (cid:105) . (18)Originally the IPS method comes from filtering methods. Filtering methods aim at estimat-ing the target measures and, in these methods, the potential functions are chosen so thatthe η k match the target measures. But in our cases we have no interest in estimating thetarget measures, we only want to estimate p h . So the potential functions can be chosen10ore freely. They are used to propose a probabilistic representation of p h in terms of aselection+mutation dynamics, which makes it possible to build an estimator of p h with areduced variance.We define the propagated target measures η k such that η = ˜ η and for k ≥ η k +1 =˜ η k V k . We have : η k +1 ( d z τ k +1 ) ∝ k (cid:89) s =0 G s ( z τ s ) k (cid:89) s =0 V s ( d z τ s +1 | Z τ s ) , (19)or equivalently ∀ B ∈ E τ k +1 , η k +1 ( B ) = E (cid:104) B ( Z τ k +1 ) (cid:81) ks =0 G s ( Z τ s ) (cid:105) E (cid:104)(cid:81) ks =0 G s ( Z τ s ) (cid:105) . (20)For k = 0 we consider that η = δ , but the methods would still be valid if we had η (cid:54) = δ .We define Q k such that for f ∈ M ( E τ k +1 ), Q k ( f )( Z τ k ) = (cid:90) E τk +1 f ( z τ k +1 ) V k ( d z τ k +1 | Z τ k ) G k ( Z τ k )and set Q k,n = Q k Q k +1 . . . Q n . Let Ψ k be the application that transforms a measure η defined on E τ k into a measure Ψ k ( η ) defined on E τ k as follows:Ψ k ( η )( f ) = (cid:82) G k ( z ) f ( z ) dη ( z ) η ( G k ) . (21)We say that Ψ k ( η ) gives the selection of η through the potential G k . Notice that ˜ η k is theselection of η k as ˜ η k = Ψ k ( η k ). The target distributions can therefore be built according tothe following pattern, in which a propagation step follows a selection step: η k Ψ k −−−−→ ˜ η k .V k −−−−→ η k +1 . We also define the associated unnormalized measures ˜ γ k and γ k +1 , such that for f ∈ M ( E τ k ):˜ γ k ( f ) = E (cid:34) f ( Z τ k ) k (cid:89) s =0 G s ( Z τ s ) (cid:35) and ˜ η k ( f ) = ˜ γ k ( f )˜ γ k (1) , (22)and for f ∈ M ( E τ k +1 ): γ k +1 ( f ) = E (cid:34) f ( Z τ k +1 ) k (cid:89) s =0 G s ( Z τ s ) (cid:35) and η k +1 ( f ) = γ k +1 ( f ) γ k +1 (1) . (23)Denoting f h ( Z τ n ) = h ( Z τn ) (cid:81) n − s =0 G s ( Z τs ) , notice that we have: p h = γ n ( f h ) = η n ( f h ) n − (cid:89) k =0 η k (cid:0) G k (cid:1) . (24)11 .2. The IPS algorithm and its estimators The IPS method provides an algorithm to generate weighted samples which approximatethe probability measures η k and ˜ η k respectively for each step k . For the sample approxi-mating η k , we denote Z jτ k the j th trajectory and W jk its weight. Respectively, for the sampleapproximating ˜ η k , we denote ∼ Z jτ k the j th trajectory and ∼ W jk its associated weight. For sim-plicity reasons, in this paper, we consider that the samples all contain N trajectories, butit is possible to modify the sample size at each step, as illustrated in [17]. The empiricalapproximations of η k and ˜ η k are denoted by η Nk and ˜ η Nk and are defined by:˜ η Nk = N (cid:88) i =1 ∼ W ik δ ∼ Z iτk and η Nk = N (cid:88) i =1 W ik δ Z iτk . (25)So for all k ≤ n and f ∈ M ( E τ k ),˜ η Nk ( f ) = N (cid:88) i =1 ∼ W ik f (cid:0) ∼ Z iτ k (cid:1) and η Nk ( f ) = N (cid:88) i =1 W ik f (cid:0) Z iτ k (cid:1) . (26)By plugging these estimations into equations (22) and (22), we get estimations for theunnormalized distributions. Denoting by ˜ γ Nk and γ Nk these estimations, for all k ≤ n and f ∈ M ( E τ k ), we have:˜ γ Nk ( f ) = ˜ η Nk ( f ) k − (cid:89) s =0 η Ns ( G s ) and γ Nk ( f ) = η Nk ( f ) k − (cid:89) s =0 η Ns ( G s ) . (27)Plugging the estimations η Nk into equation (24), we get an estimator ˆ p h of p h defined by:ˆ p h = η Nn ( f h ) n − (cid:89) k =0 η Nk (cid:0) G k (cid:1) . (28)The algorithm builds the samples sequentially, alternating between a selection step and apropagation step. The k th selection step transforms the sample ( Z jk , W jk ) j ≤ N into the sample( ∼ Z jk , ∼ W jk ) j ≤ N . This transformation is done with a multinomial resampling scheme: the ∼ Z jk ’sare drawn with replacement from the sample ( Z jk ) j ≤ N , each trajectory Z jk having a proba-bility W jk G k ( Z jτk ) (cid:80) Ni =1 W ik G k ( Z iτk ) to be drawn each time. We denote by A jk the ancestor index of the j th trajectory in the selected sample, such that ∼ Z jk = Z A jk τ k . We let ˜ N jk = card { i, A ik = j } be thenumber of times the particle Z jk is replicated in the sample ( ∼ Z jk , ∼ W jk ) j , so N = (cid:80) Nj =1 ˜ N jk . Af-ter this resampling the weights ∼ W jk are set to N . The interest of this selection by resamplingis that it discards low potential trajectories and replicates high potential trajectories. So theselected sample focuses on trajectories that will have a greater impact on the estimations ofthe next distributions once extended. 12 nitialization : k = 0 , ∀ j = 1 ..N, Z j i.i.d. ∼ η and W j = N , and ∼ W j = G ( Z j ) (cid:80) s G ( Z s ) while k < n doSelection: Sample ( ˜ N jk ) j =1 ..N ∼ M ult (cid:0) N, ( ∼ W jk ) j =1 ..N (cid:1) ∀ j := 1 ..N, ∼ W jk := N Propagation :for j := 1 ..N do Sample Z jτ k +1 from V k +1 ( . | ∼ Z jτ k )set W jk +1 = ∼ W jk for i := 1 ..N do Set ∼ W ik +1 = W ik +1 G k +1 ( Z iτk +1 ) (cid:80) j W jk +1 G k +1 ( Z jτk +1 ) if ∀ j, ∼ W jk +1 = 0 then ∀ q > k, set η Nq = ˜ η Nq = 0 and Stop else k := k + 1 Figure 3: IPS algorithm
If one specifies potential functions that are not positive, there can be a possibility that ata step k we get ∀ j, G k ( Z jτ k ) = 0 , and so the probability for resampling cannot be defined.When this is the case, the algorithm stops and we consider that ∀ s ≥ k the measures ˜ η Ns and η Ns are equal to the null measure.Then the k th propagation step transforms the sample ( ∼ Z jk , ∼ W jk ) j ≤ N , into the sample( Z jk +1 , W jk +1 ) j ≤ N . Each trajectory Z jk +1 is obtained by extending the trajectory ∼ Z jk on theinterval [ τ k , τ k +1 ) using the transition kernel V k . The weights satisfy W jk +1 = ∼ W jk , ∀ j . Thenthe procedure is iterated until the step n . The full algorithm to build the samples is dis-played in Figure 3.For the sake of simplicity, we will make the following assumption: ∃ ε , ε ∈ R + suchthat ∀ Z τ k ∈ E τ k : ε > G k ( Z τ k ) > ε > , ((G)) Theorem 1.
When (G) is verified the estimator (28) is unbiased and strongly consistent.
The proof of theorem 1 follows from Theorems 7.4.2 and 7.4.3 in [8].
Theorem 2.
When (G) is verified: √ N (ˆ p h − p h ) d −→ N →∞ N (cid:0) , σ IP S,G (cid:1) , (29)13 here σ IP S,G = n − (cid:88) k =0 γ k (1) η k (cid:16)(cid:2) Q k,n ( f h ) − η k Q k,n ( f h ) (cid:3) (cid:17) (30)= n − (cid:88) k =0 (cid:40) E z (cid:20) k − (cid:89) i =0 G i ( Z τ i ) (cid:21) E z (cid:20) E [ h ( Z τ n ) | Z τ k ] k − (cid:89) s =0 G − s ( Z τ s ) (cid:21) − p h (cid:41) . (31)A proof of this CLT can be found in [8] chapter 9 at the theorem 9.3.1 . For the estimationof the variance σ IP S,G we refer the reader to [17].
We have seen that the resampling steps have the advantage of replicating high potentialtrajectories and discarding low potential trajectories. However the resampling steps alsointroduce some additional fluctuations into the estimation (see (31)). So we would liketo trigger them only when it is judicious. Typically, not when the potentials of all thetrajectories are similar, as in this case there is not point in discarding or replicating sometrajectories over others. In order to avoid pointless resampling, one can trigger the selectionstep only when the weights are unbalanced. This is done in the Sequential Monte Carlo(SMC) algorithm with adaptive resampling presented in Figure 4. In this algorithm, theheterogeneity of the weights is quantified using the effective sample size. At the k th step theeffective sample size is defined by: ESS k = (cid:16)(cid:80) Nj =0 W jk G k ( Z jτ k ) (cid:17) (cid:80) Ni =0 (cid:0) W ik G k ( Z iτ k ) (cid:1) . (32)Its value is between 1 and N and it measures the homogeneity in the candidate weights W ik G k ( Z iτk ) (cid:80) j W jk G k ( Z jτk ) : when ESS k = N the weights are perfectly balanced and are all equal to N ,and conversely when ESS k = 1 all the weights are null except one, which concentratesthe totality of the mass. Therefore, one considers the weights are too unbalanced when ESS k ≤ eN where e ∈ [0 ,
1] is a tuning parameter.Note that in the presented algorithm one can use alternative strategies to select highpotential trajectories. Here, the presented algorithms include a standard multinomial re-sampling procedure, but one can also use residual resampling or stratified resampling withoutaltering the properties of the estimator. Empirical results suggest that these alternative re-sampling schemes yield estimations with smaller variances [11, 14]. There are also recenttheoretical results on the performance of different resampling schemes [13]. MCMC stepswith invariant distribution ˜ η k can also be included in the algorithm after the resamplingstep. Some adaptations of the algorithm for parallel implementations have also been stud-ied in [19] for instance. 14 nitialization : k = 0 , ∀ i := 1 ..N, Z i = ( z ) and W i = N , and ∼ W i = G ( Z i ) (cid:80) j G ( Z j ) while k < n doSelection:if ESS k ≤ eN then Sample ( ˜ N jk ) j =1 ..N ∼ M ult (cid:0) N, ( ∼ W jk ) j =1 ..N (cid:1) and set ∀ i = 1 ..n, ∼ W ik := N elsefor i := 1 ..N do set ∼ Z iτ k := Z iτ k Propagation :for i := 1 ..N do Sample Z jτ k +1 from V k +1 ( . | ∼ Z jτ k )set W ik +1 = ∼ W ik for i := 1 ..N do Set ∼ W ik +1 = W ik +1 G k +1 ( Z iτk +1 ) (cid:80) j W jk +1 G k +1 ( Z jτk +1 ) if ∀ j, ∼ W jk +1 = 0 then ∀ q > k, set η Nq = ˜ η Nq = 0 and Stop else k := k + 1 Figure 4: SMC algorithm with adaptive resampling steps
4. Choice of the potential functions
Note that the variance of ˆ p h depend on the number of subdivisions and on the choiceof the potential functions. We display here an important hint on how to select potentialfunctions that yield a small variance. Indeed, the theoretical expressions of the potentialfunctions that minimize the asymptotic variance of the IPS estimator are known [1]. Theorem 3.
For k ≥ , let G ∗ k be defined by: G ∗ k ( z τ k ) = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) E (cid:104) E (cid:2) h ( Z τ n ) (cid:12)(cid:12) Z τ k +1 (cid:3) (cid:12)(cid:12) Z τ k = z τ k (cid:105) E (cid:104) E (cid:2) h ( Z τ n ) (cid:12)(cid:12) Z τ k (cid:3) (cid:12)(cid:12) Z τ k − = z τ k − (cid:105) (33) if E (cid:104) E (cid:2) h ( Z τ n ) (cid:12)(cid:12) Z τ k (cid:3) (cid:12)(cid:12) Z τ k − = z τ k − (cid:105) (cid:54) = 0 , and G ∗ k ( z τ k ) = 0 otherwise. For k = 0 , we define G ∗ ( z τ ) = (cid:114) E (cid:104) E (cid:2) h ( Z τ n ) (cid:12)(cid:12) Z τ (cid:3) (cid:12)(cid:12) Z τ = z τ (cid:105) . (34) The potential functions minimizing σ IP S,G are the ones that are proportional to the G ∗ k ’s k ≤ n . The optimal variance of the IPS method with n steps is then σ IP S,G ∗ = E (cid:104) E (cid:2) h ( Z τ n ) (cid:12)(cid:12) Z τ (cid:3) (cid:105) − p h + n (cid:88) k =1 E (cid:34)(cid:114) E (cid:104) E (cid:2) h ( Z τ n ) (cid:12)(cid:12) Z τ k (cid:3) (cid:12)(cid:12) Z τ k − (cid:105)(cid:35) − p h . (35)For the reliability assessment case, the optimal target distributions have the form:˜ η ∗ k ( d z τ k ) ∝ (cid:113) E (cid:2) E [ D ( Z τ n ) | Z τ k +1 ] (cid:12)(cid:12) Z τ k = z τ k (cid:3) k (cid:89) s =1 V s ( d z τ s | z τ s -1 ) , (36)and the potential functions would be defined by (33) and (34), taking h = D . As we donot have closed-form expressions of the functions z τ k → (cid:113) E (cid:2) E [ D ( Z τ n ) | Z τ s +1 ] | (cid:12)(cid:12) Z τ s = z τ k (cid:3) ,we propose to use instead a parametric approximation of these expectations based on ourknowledge of the system and denoted by U α ( Z τ s ) so that we take˜ η k ( d Z τ k ) ∝ U α ( Z τ k ) k (cid:89) s =1 V s ( d Z τ s | Z τ s -1 ) (37)and ∀ s > , G s ( Z τ s ) = U α ( Z τ s ) U α ( Z τ s -1 ) and G ( Z τ ) = U α ( Z τ ) . (38)For a system including similar components in parallel redundancy, we propose to set U α ( Z τ s ) = 1 when Z τ s has already reached the failure region once, and to set U α ( Z τ s ) =exp (cid:2) − α ( b ( Z τ s ) + 1) (cid:3) L ( τ s ) otherwise, where L is a positive function, and b ( Z ) indicatesthe number of working components within a state Z . Here α is a parameter tuning thestrength of the selection.
5. The IPS+M method for concentrated PDMPs
When it is used on a reliable system and therefore on a concentrated PDMP (see Section2.4), the IPS method tends to loose in efficiency. This efficiency loss can be attributed to theexploration steps. Remember that an exploration step comes after a selection step: it buildsa sample ( Z jτ k +1 , W jk +1 ) j ≤ N by extending the trajectories of a selected sample ( ∼ Z jτ k , ∼ W jk ) j ≤ ˜ N .This newly built sample ( Z jτ k +1 , W jk +1 ) j ≤ N fulfills two goals : 1) It contributes to the empir-ical approximation η Nk +1 of η k +1 . 2) It is used as a candidate sample for the next selection.But this second goal is often poorly achieved with a concentrated PDMP. Indeed, in orderto get a good approximation ˜ η Nk +1 of ˜ η k +1 , it is preferable that the candidate sample to se-lection ( Z jτ k +1 , W jk +1 ) j ≤ N contains as many different trajectories as possible, along with highpotential trajectories. Unfortunately, with this kind of PDMP, it is generally not the case:the candidate sample often contains several replicates of the same trajectories, and no high16otential trajectory. Therefore each distribution ˜ η k is poorly represented, and so is eachtarget distribution η k +1 , which eventually deteriorates the quality of the estimator ˆ p D or ˆ p h .To understand why the exploration steps are not likely to generate many different trajec-tories with a concentrated PDMP, we have to come back at the beginning of the propagationstep. At that point, the sample ( ∼ Z jτ k , ∼ W jk ) ≤ j ≤ N is naturally clustered because of the previousselection step, each of the clusters containing several replicates of the same trajectories. Wecan rewrite (25) in the following way˜ η Nk = N (cid:88) i =1 ∼ W ik δ ∼ Z iτk = 1 N N (cid:88) j =1 ˜ N jk δ Z jτk (39)where (cid:80) Nj =1 ˜ N jk = N . In practice many of the ˜ N jk are null and only a few are positive andthe N resampled trajectories ( ∼ Z jτ k ) j ≤ N are concentrated on a few trajectories. Then, each ofthe ˜ N jk trajectories of the j -th cluster is extended by using the same distribution V k ( . | Z jτ k ).(For all index i such that A ik = j the trajectory ∼ Z iτ k is extended with the kernel V k ( . | Z jτ k )).As the kernel V k ( . | Z jτ k ) corresponds to a concentrated PDMP, it is likely to extend all thetrajectories of a cluster in the same manner. The trajectory a k,jτ k +1 which extends Z jτ k until τ k +1 without spontaneous jump or failure concentrates the mass of the kernel V k ( . | Z jτ k ).Indeed, at this point we have : V k ( a k,jτ k +1 | Z jτ k ) = P (cid:0) Z τ k +1 = a k,jτ k +1 (cid:12)(cid:12) Z τ k = Z jτ k (cid:1) (cid:39) . (40)Therefore each of the trajectories ∼ Z iτ k = Z jτ k in a cluster tends to be extended in a k,jτ k +1 . Thus,the trajectories within a cluster are likely to stay clumped together during the propagation,and the propagated sample ( Z jτ k +1 , W jk +1 ) j ≤ N is very likely to be clustered too. When thepreponderant trajectories a k,jτ k +1 have low potential values, the sample is not likely to containhigh potential trajectories. Consequently the selection step having no good candidates andtoo few candidates, it tends to yield an inaccurate estimation of the distributions ˜ η k .This situation is typical of reliability assessment. In that context, a well constructed poten-tial function is close to G ∗ k wherein h = D . So the potential of a trajectory G k +1 ( Z τ k +1 )should be high if its final state Z τ k +1 is more degraded than the state Z τ k . This generallyimplies that Z τ k +1 includes at least one component failure between τ k and τ k +1 . As thepreponderant trajectories a k,jτ k +1 do not contain failure between τ k and τ k +1 they generallyare associated with low potential values.The segment of a k,jτ k +1 on ( τ k , τ k +1 ] relates to a Dirac contribution of the measure ζ Z jτk ,τ k +1 − τ k .We can, therefore, decompose the expected propagation of the trajectory Z jτ k in this way: δ Z jτk V k ( f ) = f (cid:0) a k,jτ k +1 (cid:1) V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1) + (cid:90) E τk +1 \{ a k,jτk +1 } f (cid:0) z τ k +1 (cid:1) V k (cid:0) d z τ k +1 | Z jτ k (cid:1) , (41)17here f ∈ M ( E τ k +1 ). And the expected propagation of η Nk would be:˜ η ˜ Nk V k ( f ) = N (cid:88) j =1 ˜ N jk N δ Z jk V k ( f )= N (cid:88) j =1 ˜ N jk N f (cid:0) a k,jτ k +1 (cid:1) V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1) + ˜ N jk N (cid:90) E τk +1 \{ a k,jτk +1 } f (cid:0) z τ k +1 (cid:1) V k (cid:0) d z τ k +1 | Z jτ k (cid:1) (cid:39) N (cid:88) j =1 ˜ N jk N f (cid:0) a k,jτ k +1 (cid:1) . (42) In order to diversify the simulated trajectories, and to increase the precision of theestimation, we propose to modify the way we extend the selected trajectories. Here weconsider that the size of the propagated sample can differ from the size of the previousselected sample. We now denote ˜ N k the size of the k th selected sample, and N k +1 the sizeof the k th propagated sample, with the convention N = N . We stressed out, in section 5.1,that the propagation step aims at providing an estimation of η k +1 = ˜ η k V k using the selectedsample ( ∼ Z jk , ∼ W jk ) j ≤ ˜ N k . In other words, the selection step aims at providing a propagatedweighted sample ( Z jτ k +1 , W jk +1 ) j ≤ N k +1 to estimate the distribution ˜ η ˜ N k k V k defined by:˜ η ˜ N k k V k ( f ) = ˜ N k (cid:88) j =1 ∼ W jk δ ∼ Z jk V k ( f ) = N k (cid:88) j =1 ˜ N jk ˜ N k δ Z jk V k ( f ) , (43)where f ∈ M ( E τ k +1 ). We denote by ¯ V k the Markovian kernel from E τ k to E τ k +1 such that,for any trajectory Z jτ k , ¯ V k ( . | Z jτ k ) is the conditioning of V k ( . | Z jτ k ) to E k \{ a k,jτ k +1 } . ¯ V k ’s densityverifies: ¯ V k ( z τ k +1 | Z jτ k ) = V k ( z τ k +1 | Z jτ k )1 − V k (cid:0) a k,jτ k +1 (cid:12)(cid:12) Z jτ k (cid:1) z τk +1 (cid:54) = a k,jτk +1 . (44)Using (43) we can decompose ˜ η ˜ N k k V k as follows:˜ η ˜ N k k V k ( f ) = N k (cid:88) j =1 ˜ N jk ˜ N k (cid:34) V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1) f (cid:0) a k,jτ k +1 (cid:1) + (cid:16) − V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1)(cid:17) δ Z jk ¯ V k ( f ) (cid:35) . (45)In the original IPS algorithm, the sample approximating ˜ η ˜ N k k V k is built by directly extend-ing each trajectory in the selected sample. When we extend the replicates of a cluster, inaverage a proportion V k ( a k,jτ k +1 | Z jτ k ) of the replicates are extended in a k,jτ k +1 . This proportionof trajectories extended in a k,jτ k +1 then serves as an estimation of V k ( a k,jτ k +1 | Z jτ k ). But it is notnecessary to misspend all these replicates to estimate the probability of the preponderant18rajectory. If we use equation (45), we would need to generate the trajectory a k,jτ k +1 only onceto assess its contribution to the propagation of the cluster. Also, a k,jτ k +1 is easy to get. Togenerate it, it generally suffices to run the simulation process starting from the state Z jτ k until time τ k +1 , while setting the jumps rates and the probability of failure on demand to zero.Therefore, for each cluster, we propose to use an additional replicate to generate a k,jτ k +1 and compute exactly its contribution. So for any j ∈ { , . . . N k } , we will extend the selectedtrajectory Z jk , N jk times, where N jk = ˜ N jk + ˜ N jk > . We denote j i the index of the i th replicates of Z jk , and consider the added replicate has index 0 such that for i ∈ { , . . . , ˜ N jk } we have ∼ Z j i τ k = Z jk . The additional replicate is deterministically extended to the preponderanttrajectory, so we have Z j τ k +1 = a k,jτ k +1 , and we set its weight to W j k +1 = V k ( a k,jτ k +1 | Z jτ k ) ˜ N jk ˜ N k , sothat it carries all the mass associated to the preponderant trajectories of a cluster. Then wecan use all the remaining ˜ N jk trajectories in the cluster to estimate the non preponderantpart of the cluster’s propagation (the 1 st term in the right hand side of equation (41)). For i >
0, we condition the extensions to avoid a k,jτ k +1 generating the Z j i τ k +1 according to thekernel ¯ V k ( . | Z jτ k ) and set W j i k +1 = − V k ( a k,jτk +1 | Z jτk )˜ N jk ˜ N jk ˜ N k .Usually, the simulations of a restricted law are carried out using a rejection algorithm, butin our case a rejection algorithm would perform poorly. The rate of rejection would be toohigh, as it would be equal to V k ( a k,jτ k +1 | Z jτ k ) which is typically close to 1. For PDMPs, suchsimulations, conditioned to avoid a preponderant trajectory, can be efficiently carried outusing the memorization method. This method, introduced in [16], shares similarities withthe inverse method. It therefore benefits from not using any rejection, and so it is well suitedto our applications. The memorization method is presented in section 6.The target distributions ˜ η k are still estimated with ˜ η ˜ N k k , using equation (39), but for k = 0 to n − η k +1 , the propagation of a target distribution, is now estimated by : η N k +1 k +1 = N k +1 (cid:88) i =1 W ik +1 δ Z ik +1 = N k (cid:88) j =1 , ˜ N jk > N jk (cid:88) i =0 W j i k +1 δ Z jik +1 = N k (cid:88) j =1 , ˜ N jk > ˜ N jk ˜ N k (cid:34) V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1) δ a k,jτk +1 + (cid:16) − V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1)(cid:17) ˜ N jk ˜ N jk (cid:88) i =1 δ Z jik +1 (cid:35) (46)Let ˜ N k = ( N , ˜ N , N , ˜ N , . . . , ˜ N k ) and N k = ( N , ˜ N , N , ˜ N , . . . , N k ) We now note ˜ γ ˜ N k k and γ N k k the estimations of the unnormalized distributions, and for all k ≤ n and f ∈ M ( E τ k ),we define them by:˜ γ ˜ N k k ( f ) = ˜ η ˜ N k k ( f ) k − (cid:89) s =0 η N s s ( G s ) and γ N k k ( f ) = η N k k ( f ) k − (cid:89) s =0 η N s s ( G s ) . (47)19n the end, p h is estimated using the equation :ˆ p h = η N n n ( f h ) n − (cid:89) k =0 η N k k (cid:0) G k (cid:1) . (48)The full modified version of the algorithm is presented in Figure 5. We call this modifiedversion of the IPS algorithm the IPS+M algorithm.Throughout the rest of the paper, the notation E M will indicate that the expectation isassociated to the IPS+M method and E will still denote the expectation for the original IPSmethod. Initialization : k = 0 , ∀ j = 1 ..N, Z j = ( z ) and W j = N , and ∼ W j = G ( Z j ) (cid:80) s G ( Z s ) while k < n doSelection: ˜ N k = N , and sample ( ˜ N jk ) j =1 ..N k ∼ M ult (cid:0) N, ( ∼ W jk ) j =1 ..N k (cid:1) ∀ i = 1 .. ˜ N k , ∼ W ik := N ∀ i = 1 .. ˜ N k , set N ik = ˜ N ik + ˜ N ik > set N k := (cid:80) N k i =1 N ik Propagation :for j := 1 ..N k doif N jk > then set Z j τ k +1 = a k,jτ k +1 and W j k +1 = V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1) (cid:80) ˜ N jk i =1 ∼ W j i k for j = 1 .. ˜ N jk do Sample Z j i τ k +1 from ¯ V k ( . | Z jτ k ) and set W j i k +1 = (1 − V k (cid:0) a k,jτk +1 | Z jτk (cid:1) )˜ N jk (cid:80) ˜ N jk i =1 ∼ W j i k for i := 1 ..N k do ∼ W ik +1 = W ik +1 G k +1 ( Z iτk +1 ) (cid:80) j W jk +1 G k +1 ( Z jτk +1 ) if ∀ j, ∼ W jk +1 = 0 then ∀ q > k, set η N q q = ˜ η ˜ N q q = 0 and Stop else k := k + 1 Figure 5: IPS+M algorithm
In this section we show that the estimator ˆ p h obtained with the IPS+M method havethe same basic properties as the IPS estimator. With the IPS+M method, ˆ p h convergesalmost surely to p h , it is unbiased, and it satisfies a CLT. The proofs that we provide inthis section follow the reasoning of the proofs in [8]. We present how to adjust the originalproofs to take into account that the extensions of the trajectories within a cluster are no20onger identically distributed. Finally we show that the asymptotic variance of the CLT isreduced with the IPS+M method. With the IPS+M method, we assumed that ∀ k, ˜ N k = N . For p ≤ n, we define F p thefiltration associated to the sequence of the p first random samples built with the IPS+Malgorithm: (cid:0) ( Z jτ ) j ≤ N , ( ∼ Z jτ ) j ≤ N , ( Z jτ ) j ≤ N , . . . , (cid:1) . So when p is an even number such that p = 2 k , F p is adapted to the vector (cid:0) ( Z jτ ) j ≤ N , ( ∼ Z jτ ) j ≤ N , . . . , ( Z jτ k ) j ≤ N k (cid:1) . For an odd num-ber p = 2 k + 1, F p is adapted to the vector (cid:0) ( Z jτ ) j ≤ N , ( ∼ Z jτ ) j ≤ N , . . . , ( Z jτ k ) j ≤ N k , ( ∼ Z jτ k ) j ≤ N (cid:1) .For f ∈ M ( E τ n ) we let Γ N p, n ( h ) be define byΓ N k, n ( f ) = γ N k k ( Q k,n ( f )) − γ k ( Q k,n ( f ))= γ N k k ( Q k,n ( f )) − γ n ( f ) (49)and Γ N k +1 , n ( f ) = ˜ γ ˜ N k k ( V k Q k +1 ,n ( f )) − ˜ γ k ( V k Q k +1 ,n ( f ))= ˜ γ ˜ N k k ( V k Q k +1 ,n ( f )) − γ n ( f ) . (50)Using a telescopic argument we getΓ N p, n ( f ) = (cid:98) p (cid:99) (cid:88) k =0 γ N k k ( Q k,n ( f )) − ˜ γ ˜ N k − k − ( V k − Q k,n ( f ))+ p> (cid:98) p +12 (cid:99) (cid:88) k =1 ˜ γ ˜ N k − k − ( V k − Q k,n ( f )) − γ N k − k − ( Q k − ,n ( f )) , (51)with the convention for k = 0 , ˜ γ ˜ N − − ( V − Q ,n ( f )) = γ n ( f ).Noticing that γ N k k (1) = ˜ γ ˜ N k − k − (1) = γ N k − k − ( G k − ), we can rewritte (51) asΓ N p, n ( f ) = (cid:98) p (cid:99) (cid:88) k =0 γ N k k (1) (cid:16) η N k k ( Q k,n ( f )) − ˜ η ˜ N k − k − V k − ( Q k,n ( f )) (cid:17) + p> (cid:98) p +1 (cid:99) (cid:88) k =1 ˜ γ ˜ N k − k − (1) (cid:16) ˜ η ˜ N k − k − ( V k − Q k,n ( f )) − Ψ k − ( η N k − k − )( V k − Q k,n ( f )) (cid:17) , (52)where for k = 0, we use the convention γ N (1)˜ η ˜ N − − ( V − Q ,n ( f )) = γ n ( f ). The benefit of thisdecomposition is that it distinguishes the errors associated to the propagation steps and theerrors associated to the selection steps. For the propagation steps, using (41) we easily getthat for any f ∈ M ( E τ k +1 ): E M (cid:104) η N k k ( f ) (cid:12)(cid:12) F k − (cid:105) = ˜ η ˜ N k − k − V k − ( f ) . (53)21or the selection steps, as the resampling schemes are the same ones as for the IPS algorithm,we still have for any f ∈ M ( E τ k ): E M (cid:104) ˜ η ˜ N k k ( f ) (cid:12)(cid:12) F k (cid:105) = Ψ k ( η Nk )( f ) . (54)Thus, each selection step and propagation step is conditionally unbiased. Note that γ N k k (1)is F k − -measurable and γ ˜ N k k (1) is F k -measurable, so, when the samples are generated withthe IPS+M algorithm, (Γ N p, n ( h )) p ≤ n is a F p -martingale. Therefore, ˆ p h stays unbiased withthe IPS+M method, because E M [Γ N n, n ( f h )] = E M [ˆ p h − p h ] = 0 . Thanks to this martingale decomposition, we can use the same arguments as in the proofin the chapter 7 in [8]. The hypotheses of Theorems 7.4.2 and 7.4.3 [8, page 239 and 241]are satisfied with the IPS+M method too, which yields the following theorem:
Theorem 4.
For any h ∈ M ( E τ n ) , ˆ p h converges almost surely to p h , and, for any f ∈ M ( E τ k ) , η N k k ( f ) converges to η k ( f ) almost surely, γ N k k ( f ) converges to γ k ( f ) almost surely.5.3.3. A Central Limit Theorem Theorem 5.
If the potential functions verify (G) and the samples are generated with theIPS+M algorithm ,then we have the following convergence in distribution: √ N (ˆ p h − p h ) −→ N →∞ N (cid:0) , σ M,G (cid:1) , where σ M,G = η (cid:18)(cid:104) Q ,n ( f h ) − η Q ,n ( f h ) (cid:105) (cid:19) + n (cid:88) k =1 γ k (1) ˜ η k - (cid:18)(cid:0) − V k − ( a τ k | Z τ k - ) (cid:1) ¯ V k - (cid:104) Q k,n ( f h ) − ¯ V k - Q k,n ( f h ) (cid:105) (cid:19) + n (cid:88) k =1 ˜ γ k (1) ˜ η k - (cid:18)(cid:104) V k - Q k,n ( f h ) − ˜ η k - V k - Q k,n ( f h ) (cid:105) (cid:19) . (55) Proof.
This proof is very similar to what is done in Chapter 9 of [8]. In order to provethat ˆ p h satisfies a CLT, we begin by proving that the errors associated to the selection andpropagation steps are normally distributed using Lindeberg’s theorem.For a sequence of function ( f k ) k ≤ n such that f k and f k +1 are in M ( E τ k ), we define thesum of errors until the p th selection and propagation by: M N p, n ( f ) = (cid:98) p (cid:99) (cid:88) k =0 η N k k ( f k ) − ˜ η ˜ N k − k − V k − ( f k )+ p> (cid:98) p +12 (cid:99) (cid:88) k =1 ˜ η ˜ N k − k − ( f k − ) − Ψ k − ( η N k − k − )( f k − ) . (56)22or j ∈ { , . . . N } we let U N (2 k +1) N + j ( f ) = 1 √ N (cid:16) f k +1 ( ∼ Z jk ) − Ψ k ( η N k k )( f k +1 ) (cid:17) . (57)For k ≥ j ∈ { , . . . , N k } and i ∈ { , . . . , N j } , we consider that the indices j i are orderedin such way that j > N and j i < N when i >
0. With such indexing ∀ s ∈ { , . . . , N } , ∃ j ∈ { , . . . , N k } and i ∈ { , . . . , N j } such that s = j i , and for such s we let U N kN + s ( f ) = 1 − V k (cid:0) a k,jτ k +1 | Z jτ k (cid:1) √ N (cid:0) f k +1) ( Z j i k +1 ) − ¯ V k ( f k +1) )( Z jk ) (cid:1) . (58)For j ∈ { , . . . , N } , let U Nj ( f ) = 1 √ N (cid:0) f ( Z j ) − η ( f ) (cid:1) . (59)Thus, √ N M N p, n ( f ) = ( p +1) N (cid:88) k =0 U Nk ( f ) . (60)Noting P Nk a filtration adapted to the k first trajectories generated in the IPS+M algorithm.Note that we have that E (cid:2) U Nk ( f ) |P Nk − (cid:3) = 0, and E (cid:2) U Nk ( f ) |P Nk − (cid:3) < ∞ , and | U Nk ( f ) | < √ N sup k ≤ n, Z τk ∈ E τk {| f k ( Z τ k ) |∧| f k +1 ( Z τ k ) |} , so the Lindeberg condition is clearly satisfied. Then,we have that (cid:104)√ N M N p, n ( f ) (cid:105) p = ( p +1) N (cid:88) k =0 E (cid:2) U Nk ( f ) |P Nk − (cid:3) = η N (cid:18)(cid:104) f − η N ( f ) (cid:105) (cid:19) + (cid:98) p (cid:99) (cid:88) k =1 ˜ η Nk -1 (cid:18)(cid:0) − V k − ( a τ k | Z τ k -1 ) (cid:1) ¯ V k -1 (cid:104) f k − ¯ V k -1 f k (cid:105) (cid:17)(cid:19) + (cid:98) p +1 (cid:99) (cid:88) k =1 ˜ η Nk -1 (cid:18)(cid:104) f k − ( Z τ k -1 ) − Ψ k -1 ( η N k -1 k -1 ) f k − (cid:105) (cid:19) . (61)As η N k k and ˜ η Nk converge almost surely to η k and ˜ η k , (cid:104)√ N M N p, n ( f ) (cid:105) n converge in probability23o σ p ( f ) = η (cid:18)(cid:104) f − η ( f ) (cid:105) (cid:19) + (cid:98) p (cid:99) (cid:88) k =1 ˜ η k -1 (cid:18)(cid:0) − V k − ( a τ k | Z τ k -1 ) (cid:1) ¯ V k -1 (cid:104) f k − ¯ V k -1 f k (cid:105) (cid:17)(cid:19) + (cid:98) p +12 (cid:99) (cid:88) k =1 ˜ η k -1 (cid:18)(cid:104) f k − ( Z τ k -1 ) − ˜ η k -1 f k − (cid:105) (cid:19) . (62)By application of the Lindeberg’s theorem for triangular array (see for instance Theorem 4on page 543 in [18]), we get that √ N M N p, n ( f ) converges in law to a centered Gaussian ofvariance σ p ( f ). As a corollary, if for p (cid:54) = 2 k we take f p = 0 and for p = 2 k f k = Q k,n ( f h ),we get that √ N (cid:18) η N k k Q k,n ( f h ) − ˜ η Nk − V k − Q k,n ( f h ) (cid:19) −→ N →∞ ˜ η k -1 (cid:18)(cid:0) − V k − ( a τ k | Z τ k -1 ) (cid:1) ¯ V k -1 (cid:104) Q k,n ( f h ) − ¯ V k -1 Q k,n ( f h ) (cid:105) (cid:17)(cid:19) and if for p (cid:54) = 2 k − f p = 0 and for p = 2 k − f k − = V k − Q k,n ( f h ), we get that √ N (cid:18) ˜ η Nk − ( V k − Q k,n ( f h )) − Ψ k − ( η N k − k − )( V k − Q k,n ( f h )) (cid:19) −→ N →∞ ˜ η k -1 (cid:18)(cid:104) V k − Q k,n ( f h ) − ˜ η k -1 V k − Q k,n ( f h ) (cid:105) (cid:19) . Following from Theorem 4, γ N k k (1) and ˜ γ Nk (1) converges almost surely to γ k (1) and ˜ γ k (1) , by an application of Slutsky’s Lemma, we get that √ N Γ N N, n ( f h ) converges in law to acentered Gaussian with variance σ M,G = γ (1) η (cid:18)(cid:104) Q ,n ( f h ) − η Q ,n ( f h ) (cid:105) (cid:19) + n (cid:88) k =1 γ k (1) ˜ η k -1 (cid:18)(cid:0) − V k − ( a τ k | Z τ k -1 ) (cid:1) ¯ V k -1 (cid:104) Q k,n ( f h ) − ¯ V k -1 Q k,n ( f h ) (cid:105) (cid:17)(cid:19) + n (cid:88) k =1 ˜ γ k -1 (1) ˜ η k -1 (cid:18)(cid:104) V k -1 Q k,n ( f h ) − ˜ η k -1 V k -1 Q k,n ( f h ) (cid:105) (cid:19) . (63)24 .3.4. Variance reduction Theorem 6.
The variance of the original IPS can be decomposed as follows: σ IP S,G = σ M,G + n (cid:88) k =1 γ k (1) ˜ η k - (cid:18) v k ( Z τ k - ) ¯ V k - (cid:16)(cid:104) Q k,n ( f h )( a τ k ) − Q k,n ( f h )( Z τ k ) (cid:105) (cid:17)(cid:19) , (64) where v k ( Z τ k - ) = V k - ( a τ k | Z τ k - ) (cid:0) − V k - ( a τ k | Z τ k - ) (cid:1) . Therefore we have σ M,G ≤ σ IP S,G .Proof. σ IP S,G = n (cid:88) k =0 γ k (1) η k (cid:16)(cid:2) Q k,n ( f h ) − η k Q k,n ( f h ) (cid:3) (cid:17) (65)= η (cid:18)(cid:104) Q ,n (( f h ) − η Q ,n (( f h ) (cid:105) (cid:19) + n (cid:88) k =1 γ k (1) ˜ η k − V k − (cid:16)(cid:2) Q k,n ( f h ) − V k − Q k,n ( f h ) + V k − Q k,n ( f h ) − η k Q k,n ( f h ) (cid:3) (cid:17) = η (cid:18)(cid:104) Q ,n (( f h ) − η Q ,n (( f h ) (cid:105) (cid:19) + n (cid:88) k =1 γ k (1) ˜ η k − V k − (cid:18)(cid:104) Q k,n ( f h ) − V k − Q k,n ( f h ) (cid:105) (cid:19) + n (cid:88) k =1 ˜ γ k − (1) ˜ η k − (cid:18)(cid:104) V k − Q k,n ( f h ) − ˜ η k − V k − Q k,n ( f h ) (cid:105) (cid:19) (66)Temporarily using the notation V k -1 ( a τ k | Z τ k -1 ) = p k , for any f ∈ M ( E τ k ), we get V k − (cid:18)(cid:104) f ( Z τ k ) − V k − f (cid:105) (cid:19) = V k − (cid:18)(cid:104) f ( Z τ k ) − p k f ( a τ k ) − (1 − p k ) ¯ V k − ( f ) (cid:105) (cid:19) = V k − (cid:0) f ( Z τ k ) − p k f ( Z τ k ) f ( a τ k ) − − p k ) f ( Z τ k ) ¯ V k − f + p k f ( a τ k ) + 2 p k f ( a τ k ) ¯ V k − f + (1 − p k ) (cid:0) ¯ V k − f (cid:1) (cid:17) = p k f ( a τ k ) + (1 − p k ) ¯ V k − ( f ) − p k f ( a τ k ) − p k (1 − p k ) f ( a τ k ) ¯ V k − ( f ) − p k (1 − p k ) f ( a τ k ) ¯ V k − f − − p k ) f ¯ V k − ( f ) + p k f ( a τ k ) + 2 p k f ( a τ k ) ¯ V k − f + (1 − p k ) ¯ V k − ( f ) = p k (1 − p k ) (cid:2) f ( a τ k ) − f ( a τ k ) ¯ V k − f + ¯ V k − ( f ) (cid:3) + (1 − p k ) (cid:0) ¯ V k − ( f ) − ¯ V k − ( f ) (cid:1) = p k (1 − p k ) ¯ V k − (cid:16)(cid:2) f ( a τ k ) − f ( Z τ k ) (cid:3) (cid:17) + (1 − p k ) (cid:16) ¯ V k − ( f ) − ¯ V k − ( f ) (cid:17) . (67)25n particular, for f = Q k,n ( f h ) we get V k − (cid:18)(cid:104) Q k,n ( f h ) − V k − Q k,n ( f h ) (cid:105) (cid:19) = p k (1 − p k ) ¯ V k − (cid:16)(cid:2) Q k , n ( f h )( a τ k ) − Q k , n ( f h ) (cid:3) (cid:17) + (1 − p k ) (cid:16) ¯ V k − ( f ) − ¯ V k − ( f ) (cid:17) (68)Plugging (68) into the second line of (66), yields (64).
6. Efficient extension of the trajectories using the Memorization method
This section presents the memorization method that was first introduced in [15]. Remem-ber that we considered that a trajectory a t is preponderant whenever p a t = P ( Z t = a t ) > a t , the memorization method allows togenerate a trajectory Z t which differs from this preponderant trajectory a t . The interest of the method, compared to a rejection algorithm, is that we generate atrajectory Z t (cid:54) = a t in one shot, whereas a rejection algorithm may generate several timesthe preponderant trajectory a t before generating a trajectory different from a t . This isespecially interesting when the probability p a t = P ( Z t = a t ) is close to 1, as, with a rejectionalgorithm, the average number of tries to get a trajectory different from a t would be − p a t which is then very high. Therefore with a rejection algorithm much computational effortwould be wasted generating a t over and over. Note that the IPS+M, greatly unbalances the weights of the propagated samples. Con-sequently it is useless to consider an algorithm which triggers a resampling step, when thevalue of on the effective sample size is below a threshold. Indeed, the weights are so unbal-anced that the effective sample size would always be very low and would trigger a resamplingat each step.
The key idea of the memorization method is to consider the stopping time τ defined suchthat: ∀ s < τ, Z s = a s and Z τ (cid:54) = a τ . (69)This time τ is the time at which the trajectory Z t differentiates itself from a t . So, to generate Z t knowing τ ≤ t is equivalent to generate Z t knowing it differs from a t . In order to simulatea trajectory Z t avoiding a t , one can follow these three steps:1. generate τ knowing τ ≤ t , and set Z τ − = a τ − ,2. generate Z τ knowing Z τ (cid:54) = a τ ,3. generate the rest of the trajectory normally until t .These steps are not difficult to realize, except for the first one.26 .3.2. Generate τ knowing τ ≤ t To achieve this first step, the authors in [15] propose to generate τ knowing τ ≤ t byusing a method equivalent to the inverse transform sampling method. We present hereafterthe theoretical foundation for this method. We denote by F the cdf of τ knowing τ ≤ t : F ( v ) = P (cid:0) τ < v | τ ≤ t (cid:1) , (70)and we denote by F − its generalized inverse defined by F − ( x ) = inf v> { v | F ( v ) ≥ x } . (71)We also denote by ˜ F the function defined by˜ F ( v ) = P ( Z v − = a v − ) = n ( a v ) (cid:89) k =0 exp (cid:104) − Λ a sk ( t k ) (cid:105) n ( a v ) (cid:89) k =1 (cid:16) K a − sk ( a s k ) (cid:17) tk> (72)where Θ v ( a v ) = (cid:0) ( a s k , t k ) (cid:1) ≤ k ≤ n ( a v ) . Note that ˜ F is discontinuous in each jump times s k where K z − k ( z k ) (cid:54) = 1, so the inverse of ˜ F is not necessarily defined everywhere on [ p a t , F − , the generalized inverse of ˜ F defined by˜ F − ( x ) = sup v> { v | ˜ F ( v ) ≤ x } . (73)˜ F − extends the inverse of ˜ F constantly where it is not defined, this extension being donefrom the left so that ˜ F − is right continuous.The inverse transform sampling method consists in generating U ∼ U nif (0 ,
1) and taking F − ( U ) as a realization of τ | τ ≤ t which is a truncated random variable. The simulation ofsuch random variables is also presented in [10]. Note that the expression of the cdf F canbe related to ˜ F , indeed we have: ∀ v < t, F ( v ) = P ( τ < v ) P ( τ ≤ t ) = 1 − P ( τ ≥ v )1 − P ( τ > t ) = 1 − P ( Z v − = a v − )1 − P ( Z t = a t ) = 1 − ˜ F ( v )1 − p a t . Consequently we have that F − ( U ) = inf v> { v | F ( v ) ≥ U } = sup v> (cid:8) v | ˜ F ( v ) ≤ − U (cid:0) − p a t ) (cid:1)(cid:9) , = ˜ F (cid:16) − U (cid:0) − p a t (cid:1)(cid:17) . (74)Also, as U has uniform distribution on [0 , U = 1 − U (cid:0) − p a t ) is a uniform on [ p a t , U ∼ U nif ( p a t ,
1) and taking ˜ F − ( ˜ U ) as a realization of τ | τ ≤ t .27ssuming we first generate the trajectory a t and generate ˜ U according to a uniformdistribution on ( p a t , F − ( ˜ U ). We consider that duringthe generation of a t , we computed and memorized P ( Z s − k = a s − k ) and P ( Z s k = a s k ) for eachjump in the trajectory, and also for P ( Z t = a t ). Then we distinguish two cases: eitherthere exists k ≤ n ( a t ) such that P ( Z s − k = a s − k ) ≥ ˜ U > P ( Z s k = a s k ) < , either there exists k ≤ n ( a t ) such that P ( Z s k = a s k ) ≥ ˜ U > P ( Z s k = a s − k +1 ) where we take the convention that s − n ( a t )+1 = t . The first case is quite simple as by definition of ˜ F − we get ˜ F − ( ˜ U ) = s k . In thesecond case, ˜ F being continuous and strictly decreasing on [ s k , s k +1 ), it is inversible on thisinterval, and ˜ F − corresponds to F ’s inverse on ( ˜ F ( s − k +1 ) , ˜ F ( s k )]. So ˜ F − ( ˜ U ) ∈ [ s k , s k +1 )and ˜ F ( ˜ F − ( ˜ U )) = ˜ U . Notice that ∀ v ∈ [ s k , s k +1 ) , ˜ F ( v ) = ˜ F ( s k ) × exp (cid:2) − Λ a sk ( v − s k ) (cid:3) . (75)So in particular, for v = ˜ F − ( ˜ U ), we have :˜ U = ˜ F ( ˜ F − ( ˜ U )) = ˜ F ( s k ) × exp (cid:104) − Λ a sk ( ˜ F − ( ˜ U ) − s k ) (cid:105) , (76)or equivalently log (cid:32) ˜ F ( s k )˜ U (cid:33) = (cid:90) ˜ F − ( ˜ U ) − s k λ a sk ( u ) du. (77)To determine ˜ F − ( ˜ U ) we look for the value s such that the integral (cid:82) s λ a sk ( u ) du is equal tolog (cid:32) ˜ F ( s k )˜ U (cid:33) by dichotomy, then we set ˜ F − ( ˜ U ) = s k + s .To sum up the generation of a realization of τ | τ ≤ t we proceed as follow:1. Generate ˜ U ∼ U nif ( p a t , k = 02. If P ( Z s k = a s k ) ≥ ˜ U > P ( Z s − k +1 = a s − k +1 ), then we find s ∈ [0 , s k +1 − s k ) such thatlog (cid:32) ˜ F ( s k )˜ U (cid:33) = (cid:90) s λ a sk ( u ) du, and we set τ = s k + s .3. If P ( Z s − k +1 = a s − k +1 ) ≥ ˜ U > P ( Z s k +1 = a s k +1 ), then τ = s k +1
4. If the condition above were not satisfied, set k = k + 1, if k ≤ n ( a s ) repeat the steps2 to 4 In the IPS+M algorithm, assuming we apply the memorization method on an interval( τ k -1 , τ k ] and for the i th cluster, we apply it knowing that Z τ k -1 = a ik -1 . Note that this28rajectory a ik -1 is not necessarily preponderant but when we extend it into a ik , the piece oftrajectory a i ( τ k -1 ,τ k ] ) is preponderant because : P ( Z ( τ k -1 ,τ k ] = a i ( τ k -1 ,τ k ] (cid:12)(cid:12) Z τ k -1 = a ik -1 ) = P ( Z τ k = a iτ k (cid:12)(cid:12) Z τ k -1 = a ik -1 ) > . (78)So, in the IPS we try to generate trajectories of the cluster that verify Z τ k -1 = a ik -1 , butavoid the piece of trajectory a i ( τ k -1 ,τ k ] .
7. Numerical illustrations
In order to confirm our results numerically, we have applied the IPS method and theIPS+M method to two two-components system.
The first system is a room heated by two heaters in passive redundancy. Heaters are pro-grammed to maintain the temperature of the room above negative values, turning on whenthe temperature drops below some positive threshold and turning off when the temperaturecrosses a high threshold. The second heater can activate only when the first one is failed.The system fails when the temperature falls below zero. X t represents the temperature of the room at time t . M t represents the status of theheaters at time t . Heaters can be on, off, or out-of-order, so M = { ON, OF F, F } . Thestate of the system is Z t = ( X t , M t ).The differential equation rule the temperature can be derived from the physics. x e is theexterior temperature. β is the rate of the heat transition with the exterior. β is the heatingpower of each heater. The differential equation giving the evolution of the the temperatureof the room has the following form: d X t dt = β ( x e − X t ) + β M t or M t = ON . The heaters are programmed to maintain the temperature within an interval( x min , x max ) where x e < < x min . We consider that the two heaters are in passive re-dundancy in the sense that: when X ≤ x min the second heater activates only if the first oneis failed. When a repair of a heater occurs, if X ≤ x min and the other heater is failed, thenthe heater status is set to ON , else the heater status is set to OF F . To handle the pro-gramming of the heaters, we set Ω m = ( −∞ , x max ) when all the heaters are failed m = ( F, F )or when at least one is activated, otherwise we set Ω m = ( x min , x max ) .Due to the continuity of the temperature, the reference measure for the Kernel is ∀ B ∈ B ( E ) , ν ( x,m ) ( B ) = (cid:80) m + ∈ M \{ m } δ ( x,m + ) ( B ). On the top boundary in x max , heaters turnoff with probability 1. On the bottom boundary in x min , when a heater is supposedto turn on, there is a probability γ = 0 .
01 that the heater will fail on demand. So,for instance, if z − = (cid:0) x min , ( OF F, OF F ) (cid:1) , we have K z − (cid:0) x min , ( ON, OF F ) (cid:1) = 1 − γ ,and K z − (cid:0) x min , ( F, ON ) (cid:1) = γ (1 − γ ), and K z − (cid:0) x min , ( F, F ) (cid:1) = γ .Let j be a transition from m to m + . For the spontaneous jumps that happen outside bound-aries, if the transition j corresponds to the failure of a heater, then:29 j ( x, m ) = 0 . . × x and, if the transition corresponds to a repair, then λ j ( x, m ) =0 . M j = F . Here the system failure occurs when the temperature of the room fallsbelow zero, so D = { ( x, m ) ∈ E, x < } . A possible trajectory of the state of this systemis plotted in figure 6. The probability of failure p was estimated to 2 . × − thanks to amassive Monte-Carlo of 10 simulations. Figure 6: A possible trajectory of the heated-room system(the mode is represented with colors)
The results of the simulation study for the heated-room system are displayed in table1. Here we have used the potential functions proposed in section 4. The value of α wasset to 1.1. We have tried different values of α between 0 . . .
1. Thevalue of α = 1 . n = 10, it reduces the variance bya factor 2 . . N = 10 the IPS+M is about 2 . We have seen that it is possible to improve the IPS method to make it similar to theSMC method. We may, therefore, think that the IPS+M algorithm could be improved byadding adaptive optional re-sampling steps in order to get a SMC+M algorithm. In practice,however, it is not beneficial to add these adaptive optional re-sampling steps. Indeed wehave noticed that, as we greatly modify the propagation process, the weights are greatly30C IPS IPS+Mˆ p . × − ˆ σ . × − n = 5 ˆ p . × − . × − ˆ σ . × − . × − n = 10 ˆ p . × − . × − ˆ σ . × − . × − Table 1: Empirical means and empirical variances on 100 runs with N = 10 for the MC, the IPS and theIPS+M methods imbalanced and the effective-sample-size ends up being extremely small, which would triggerthe re-sampling each time. Therefore adding adaptive optional re-sampling to the IPS+Mhas no effect, and in practice the IPS+M methods and the SMC+M methods are the same. The second system models a dam subjected to an incoming water flow. The physicalvariable of interest is the water level in the dam denoted by X t . The failure of the systemoccurs when the water level exceeds a security threshold x lim = 10 before time t f = 50. Theinitial level is set to X = 0. The water flow is characterized by the input debit Q = 10.The dam has two evacuation valves with output debit Q . Each valve can be either open,close or stuck closed. So M = { Open, Closed, Stuckclosed } . The valves are programmed inpassive redundancy, so if the valves are in functioning order there is always one valve openand one valve closed. Though, the valve can get stuck closed and this happens at randomtimes with exponential distribution with intensity λ = 0 . µ = 0 .
1. When both valves are stuck closed the reservoir of the dam startsfilling up according to the equation dX t dt = Q/S , where S = 10 is the surface of the reservoir. The results of the simulation study for the dam system are displayed in table 2. Herewe have used the potential functions: ∀ k < n, G ( Z τ k ) = exp (cid:2) α ( x lim − X τ k ) + α ( b ( Z τ k ) + 1) (cid:3) . (79)The value of α was set to − . α was set to − p . × − ˆ σ . × − n = 5 ˆ p . × − . × − ˆ σ . × − . × − Table 2: Empirical means and empirical variances on 50 runs with N = 10 for the MC, the IPS and theIPS+M methods estimator. In terms of computational cost, on this example the IPS+M method was 3.6times slower than the IPS, and 11.8 times slower then than the Monte-Carlo method. Sothe efficiency of the IPS+M is about 40 lower than the Monte-Carlo method. Clearly, theimplementation of the IPS+M method requires a careful choice of the form of the potentialfunctions and of their parameters.
8. Conclusion
This paper investigates the application of the IPS method to PDMPs. As the IPS methoddoes not perform well when it is used on a concentrated PDMP, we introduce and analyzethe IPS+M method, that is a modified version of the IPS that performs better with concen-trated PDMP. The IPS+M method is similar to the IPS but has different propagation steps.Its propagation steps focus on clusters of identical particles rather then on particles individ-ually. For each cluster a memorization method is used to get an empirical approximationof the distribution of the propagated cluster, which allows to greatly improve the accuracyof the method. We have shown that the proposed algorithm yields a strongly consistentestimation, and that this estimation satisfies a TCL. We prove that the asymptotic varianceof the IPS+M estimator is always smaller than the asymptotic variance of the IPS estima-tor. Simulations also confirm these results, showing that the IPS+M can yield a variancereduction when the IPS cannot. In terms of computational cost, our implementations of theIPS+M method give approximately the same efficiency as the Monte-Carlo method in theexamples considered in this paper, where the goal is to estimate a probability of the order of10 − for rather simple toy models. The numerical implementations certainly deserve morecareful attention. We also believe that there are ways to improve the efficiency of the IPS+Mmethod by finding a better class of potential functions. Another interesting improvementto the IPS+M method would be to propose an estimator of the variance. We believe thatit should be possible to adapt one of the estimators proposed in [17] for the IPS method inorder to get an estimator of the variance for the IPS+M estimator.32 eferences [1] H. Chraibi, A. Dutfoy, T. Galtier, and J. Garnier. Optimal input potential functions in the interactingparticle system method. Archive ouverte HAL , (hal-01922264), 2018.[2] H. Chraibi, J.-C. Houbedine, and A. Sibler. Pycatshoo: Toward a new platform dedicated to dynamicreliability assessments of hybrid systems. PSAM congress, 2016.[3] B. Cloez, R. Dessalles, A. Genadot, F. Malrieu, A. Marguet, and R. Yvinec. Probabilistic and piecewisedeterministic models in biology.
ESAIM: Proceedings and Surveys , 60:225–245, 2017.[4] M.H.A. Davis. Piecewise-deterministic markov processes: A general class of non-diffusion stochasticmodels.
Journal of the Royal Statistical Society. Series B (Methodological) , 46(3):353–388, 1984.[5] M.H.A. Davis.
Markov Models and Optimization . Chapman and Hall, Boca Raton, 1993.[6] B. de Saporta, F. Dufour, and H. Zhang.
Numerical Methods for Simulation and Optimization ofPiecewise Deterministic Markov Processes . Wiley, Hoboken, 2015.[7] B. De Saporta, F. Dufour, and H. Zhang.
Numerical methods for simulation and optimization ofpiecewise deterministic Markov processes: application to reliability . John Wiley & Sons, 2015.[8] P. Del Moral.
Feynman-Kac Formulae, Genealogical and Interacting Particle Systems with Applications .Springer, New York, 2004.[9] P. Del Moral and J. Garnier. Genealogical particle analysis of rare events.
The Annals of AppliedProbability , 15(4):2496–2534, 2005.[10] L. Devroye. Sample-based non-uniform random variate generation. In
Proceedings of the 18th conferenceon Winter simulation , pages 260–265. ACM, 1986.[11] R. Douc and O. Capp´e. Comparison of resampling schemes for particle filtering. In
Image and SignalProcessing and Analysis, 2005. ISPA 2005. Proceedings of the 4th International Symposium on , pages64–69. IEEE, 2005.[12] P. Fearnhead and P. Clifford. On-line inference for hidden markov models via particle filters.
Journalof the Royal Statistical Society: Series B (Statistical Methodology) , 65(4):887–899, 2003.[13] M. Gerber, N. Chopin, and N. Whiteley. Negative association, ordering and convergence of resamplingmethods. arXiv preprint arXiv:1707.01845 , 2017.[14] J.D. Hol, T.B. Schon, and F. Gustafsson. On resampling algorithms for particle filters. In
NonlinearStatistical Signal Processing Workshop, 2006 IEEE , pages 79–82. IEEE, 2006.[15] P.-E. Labeau. A Monte-Carlo estimation of the marginal distributions in a problem of probabilisticdynamics.
Reliability Engineering & System Safety , 52(1):65–75, 1996.[16] P.-E. Labeau. Probabilistic dynamics: estimation of generalized unreliability through efficient Monte-Carlo simulation.
Annals of Nuclear Energy , 23(17):1355–1369, 1996.[17] A. Lee and N. Whiteley. Variance estimation in the particle filter.
Biometrika , 105(3):609–625, 2018.[18] A.N. Shiryaev.
Probability, second edition . Springer-Verlag, volume 95 in graduate texts in mathemat-ics., New York, 1996.[19] C. Verg´e.
Mod`ele d’ˆılots de particules et application en fiabilit´e . PhD thesis, Ecole Polytechnique, 2015.[20] N. Whiteley, A.M. Johansen, and S. Godsill. Monte Carlo filtering of piecewise deterministic processes.
Journal of Computational and Graphical Statistics , 20(1):119–139, 2011., 20(1):119–139, 2011.