A regenerative modification of the multilevel splitting Alexandra Borodina and Evsey Morozov Institute of Applied Mathematical Research, Russian Academy of Sciences and Petrozavodsk University, Russia Supported by Russian Foundation for Basic Research, Grant 07-07-00088. 1
1 Introduction The multilevel splitting (MS): huge simulation time reduction vs crude Monte Carlo (CMC); statistical properties of the estimators are not studied in detail. Known: – Markov queues: the overflow probability estimate is unbiased [Garvels]; – consistency and asymptotic normality in adaptive MS [Cerou, Guyader, 2005]. 2
For a regenerative queue: splitted processes as a family of (dependent) regenerative cycles ⇒ regenerative MS . We estimate the probability of overflow within a cycle for: I.1) workload process in M/G/1 (and GI/G/1); I.2) non-Markov queue using an embedded Markov chain (MC); II. A stationary probability to exceed a level. As a result: a randomization of the thresholds and a change of algorithm. 3
2 Preliminaries I. A Markov queueing process X = { X n , n ≥ 1 } : state space E ; MS is (typically) used to estimate the probability γ c = P ( X hits set A before 0) 1. Given thresholds 0 < L 1 < · · · < L M ≡ L ; 2. Partition of E : A = C M ⊂ · · · ⊂ C 1 : C i = { x ∈ E : x ≥ L i } ; 3. R i copies upon reaching L i ; termination at L . 4
Then γ c = p 1 p 2 · · · p M with unbiased estimate M N � γ ( m ) = p i = ˆ ; (1) mR 1 · · · · · R M i =1 estimate ˆ p i of p i = P ( C i | C i − 1 ) is evaluated by CMC; N = # { reaching A = [ L, ∞ ) } ; m = the number of i.i.d. runs starting at 0. (More details: [5, 6, 7, 16].) II. A process X n , n ≥ 1 , with regeneration instants : T (0) = 0 , T ( n + 1) = min( k > T ( n ) : X k = 0) , n ≥ 0; (2) 5
the same distribution of X n + T ( k ) , n ≥ 0 , for each k ≥ 1 and independent of the X n , n < T ( k ) ; i.i.d. regeneration cycles G n = ( X k , T ( n − 1) ≤ k < T ( n )) ; i.i.d. cycle periods β n = T ( n ) − T ( n − 1) . We estimate the cycle overflow probability γ c = P ( max 1 ≤ n<β X n ∈ A ) , (3) 6
and, if weak limit X n ⇒ X exists, the stationary probability γ s = P ( X ∈ A ) , (4) for queue size and for workload in M/G/1 queue. (A generic element with no index.) I. Estimating γ c (under stability condition). Consistency: i.i.d indicators I k = I ( T ( k ) ≤ n<T ( k +1) X n ∈ A ) ; max n γ c ( n ) = 1 � I k → E I = γ c ; (5) n k =1 asymptotical normality: by regenerative CLT. 7
II. Estimating γ s : positive recurrence, E β < ∞ ; (dependent) indicators I k = I ( X k ∈ A ) and i.i.d. variables T ( i ) − 1 � Y i = I k , i ≥ 1 . (6) k = T ( i − 1) Since Y ≤ β , by regenerative theory, n γ s ( n ) = 1 I k → E Y � E β = γ s , n → ∞ , (7) n k =1 the time average E Y/ E β = P ( X ∈ A ) if weak limit X exists (e.g. if β is aperiodic). 8
Confidence 100(1 − δ )% interval for γ s (by a regenerative CLT): � � � � γ n − z δ v s ( n ) , γ n + z δ v s ( n ) √ n √ n , (8) P ( N (0 , 1) ≤ z δ ) = 1 − δ/ 2 ; where empirical variance � n i =1 ( Y i − γ s ( n ) β i ) 2 v s ( n ) = 1 ⇒ σ 2 s ≡ V ar ( Y − βγ ) , (9) 2 n β n ( β n = sample mean cycle length.) A minimal condition for the convergence: E ( Y − γβ ) 2 < ∞ . 9
(Under E Y 2 < ∞ , E β 2 < ∞ , v n is strongly consistent [8].) Applicability to dependent cycles [1, 12, 15]: below D -groups of dependent cycles with D = R 1 × · · · × R M . III. M/G/1 queue: input rate λ , service time S , distribution F , stability: ρ = λ E S < 1 . The queue size at the departure instants: a positive recurrent (aperiodic) Markov chain (MC): ν n +1 = ( ν n − 1) + + ∆ n , n ≥ 0 , ( ν 0 = 0) , (10) 10
∆ n = # { arrivals } during service of customer n + 1 . Waiting time (workload) Markov chain: W n , n ≥ 1 ; W n ⇒ W exists. If distribution of W is analytically available [Ross,Seshadri] [13]: 1 ≤ n<β W n ≥ L ) = 1 − P ( W + S ≤ L ) γ c = P ( max . (11) P ( W ≤ L ) M/M/1: ρ = λ/µ < 1 , P ( W ≤ x ) = 1 − ρe − ( λ − µ ) x : γ c = (1 − ρ ) e − ( µ − λ ) L 1 − ρe − ( µ − λ ) L , (12) 11
while (applying technique from the birth-death processes) 1 ≤ n<β ν n ≥ L ) = ρ L − 1 − ρ L P ( max 1 − ρ L . (13) NOTE 1. M/M/1 has been used to justify the quality of MS [5, 6, 7, 9, 10, 14]. If service time has subexponential integrated-tail distribution � ∞ F e ( x ) = 1 F ( y ) dy, ( F = 1 − F ) E S x and F ( xe y/ √ x ) lim = 1 , locally uniformly in y , ∈ R , (14) F ( x ) x →∞ 12
(holds for Pareto distribution) then stationary queue ν satisfies [2] 1 − ρF e ( L − 1 ρ γ s = P ( ν ≥ L ) ∼ ) , L → ∞ . (15) λ 3. Splitting and regenerations in M/G/ 1 by splitting: dependent regeneration (sub)cycles: ν ( i ) = { ν ( i ) n , τ i − 1 ≤ n < τ i } , τ 0 = 0 , i = 1 , . . . , mD ( D = R 1 × · · · × R M ) with periods α i = τ i − τ i − 1 . D -group: ( ν ( i ) : kD < i ≤ ( k +1) D ) , k = 0 , . . . , m − 1 . 13
Fig. 1: splitted processes with termination instants b, c, d (case a)) and cycles with regenerations τ n pasted together with common pre-history, with no idle-time, case b). � (t) Set�A L 2 a) L 1 L 0 a b c t d � (t) b) Set�A L 2 L 1 L 0 � � ���� t 0 � � � � 14
1. Estimating γ c = P (max 1 ≤ n<β ν n ∈ A ) : 1.a) Markov queue: standard MS, termination upon reaching A with (at most D - dependent) indicators τ i − 1 ≤ n<τ i ν ( i ) I i = I ( max n ∈ A ) , i = 1 , . . . , mD. (16) To compensate overestimation if a process starting at L i reaches 0 before L i +1 : take extra R i +1 R i +2 · · · R M (virtual) cycles, i = 1 , . . . , M − 1 . 1.b) Non-Markov M/G/1 queue, to keep independence after splitting: use MC 15
(10): makes transitions between regions G i = [ L i , L i +1 ) ⇒ : A randomization of the thresholds and modification of the algorithm: Condition A : Under transition y ∈ G i → x ∈ G k , k > i (crossing the thresholds L i +1 , . . . , L k ), generate R i +1 · · · R k processes at state x, i = 0 , . . . , M − 1; k = i + 1 , . . . , M ( R M = 1) . 2. Estimating γ s = P ( ν ≥ L ) (weak limit ν n ⇒ ν exists): evaluate the time (=number of arrivals) the process is in A using indicators I i,n = I ( ν ( i ) n ∈ A ) . Apply condition A evaluating the number of arrivals in all processes. 16
For W n , n ≥ 1 with jumps X n ( � = ± 1) : the same randomization. NOTE 2. By PASTA: estimating overflow probability for continuous time. Randomization is widely used, [3, 11] . 4. Consistency and asymptotic normality 1. Consistency: estimate γ ( m ) = regenerative estimate (5) based on n = mD dependent indicators (16), so: n I i → E � D i =1 I i 1 � = γ c . (17) n D i =1 17
NOTE 2. No need n = mD for the convergence . 2. Estimating γ s = P ( ν ≥ L ) : cycle structure, Fig 2. � 1 � 2 � m .�.�. � 1 � 2 � 3 � D � D+1 � mD ... ... .. . ... L 0 t .�.�. � 2 � D =T � mD =T � 1 � 1 m Variables τ i − 1 I ( ν ( i ) � Z i = k ≥ L ) , (18) k = τ i − 1 18
belonging to the same D -group are dependent; i.i.d variables iD � Y i = Z j , i = 1 , . . . , m ; (19) j =( i − 1) D +1 the length of the i th D -group is iD � β i = α j , j =( i − 1) D +1 so β i = T ( i ) − T ( i − 1) , i = 1 , . . . , m . Since Y ≤ β = st α 1 + · · · + α D (see Fig. 2) and E α < ∞ by ρ < 1 , then E Y < ∞ , 19
and (by regenerative theory) the estimate � n i =1 Z i → E Y γ s ( n ) = E β = P ( ν ≥ L ) ≡ γ s , n → ∞ (20) � n i =1 α i is consistent. The limit ν n ⇒ ν exists in the system M/G/1. NOTE 4. To estimate γ s in discrete-time setting one can choose any discrete- event scale because, for each cycle, the number of arrivals equals the number of departures. An asymptotic normality: by regenerative interpretation. 20
Applicability of CLT for estimating γ s = P ( ν ≥ L ) : assumption E S 2 < ∞ implies E α 2 < ∞ , Wolff [20]. Since β = st α 1 + · · · + α D then E β 2 < ∞ . Similarly, consistency and asymptotic normality of the estimates for the workload. In confidence interval (8): variance includes Y i depending on target probability. 21
5. Simulation results Fig. 3 estimate of γ s = P ( ν ≥ L ) in M/Pareto/ 1 with service time distribution F ( x ) = 1 /x 4 , x ≥ 1(= 1 , x ≤ 1) , (21) input rate λ = 0 . 45 ( ρ = 0 . 6 ), N M +1 = 10 6 . 22
9e-06 Asymptotic Splitting (regenerative) 8e-06 Splitting (regenerative, EMC) 7e-06 6e-06 5e-06 4e-06 3e-06 2e-06 1e-06 0 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 L 23
The simulation stopping rule: a given number N of the target events to be achieved. For instance, for L = 50 we take L i = { 7 , 14 , 21 , 28 , 35 , 42 } (that is M = 6 ), N = 10000 and R i ≡ 10 . More details in Table 1: S=standard MS, EMC= embedded MC. The estimate γ s ( n ) is (20) , and ”asymp” = (15). (In all experiments the number of replications m = nD is widely varied in [10 3 , 10 4 ].) 24
Recommend
More recommend
