Burkatovskaya, Yulia

Kabanova, Tatiana

Vorobeychikov, Sergey

English

Tomsk State University Repository

Филиал Кемеровского государственного университета  в г. Анжеро-Судженске Национальный исследовательский Томский государственный университет Кемеровский государственный университет Институт проблем управления им. В. А. Трапезникова РАН Институт вычислительных технологий СО РАН    ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ  И МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ (ИТММ–2015) Материалы  XIV Международной  конференции   имени  А .  Ф .  Терпугова  18–22 ноября 2015 г. Часть 1       Издательство Томского университета  2015 3 ВЕРОЯТНОСТНЫЕ И ЭКОНОМЕТРИЧЕСКИЕ МЕТОДЫ И МОДЕЛИ  CUSUM ALGORITHMS FOR PARAMETER ESTIMATION IN QUEUEING SYSTEMS WITH JUMP INTENSITY OF THE ARRIVAL PROCESS Yu. Burkatovskaya, T. Kabanova, S. Vorobeychikov  National Research Tomsk Polytechnic University, Tomsk, Russia  National Research Tomsk State University, Tomsk, Russia   1. Introduction Markovian arrival processes form a powerful class of stochastic processes introduced in [1] and [2] and thereafter they are widely used now as models for input flows to queueing systems where the rate of the arrival of customers depends on some external factors. MAP is a counting process whose arrival rate is governed by a continuous-time Markov chain. One of the problems connected with MAP is the estimation of intensity parameters by observing flow of events. A survey of estimation methods is given in [3]. Its emphasis is on maximum li-kelihood estimation and its implementation via the EM (expectation-maximization) algorithm. This approach is developed for different conditions in [4, 5], etc. The survey [6] with a huge bibliography is focused on matching mo-ment method which is also widely used for parameter estimation in MAP be-cause of its simplicity. This method is used, for example, in [7]. Bayesian ap-proach based on the a posteriori probability of the controlling chain state is de-veloped in [8]. Quality of those methods is typically examined via simulation.  In this paper we propose a different approach to MAP parameter estimation using the sequential analysis methods described in [9] and [10]. The key idea is to consider time intervals between arrivals as a stochastic process which parameters change in random instants. First we detect these points using sequential change point detection methods. Then we estimate the intensity parameters under the assumption that the intensity is constant between detected change points. 2. Problem statement We consider a Markov-modulated poisson process, i.e. a flow of events, controlled by a Markovian chain with a continuous time. The chain has two states, transition between the states happens at random instants. The time of sojourn of the chain in the l -th state is exponentially distributed with the parameter l , 1,2=l . The flow of events has the exponential distribution with the intensity parameter 1  or 2  subject to the state of the Marcovian chain. The parameters of the system 1 , 2  and the instants of switching between the states are supposed to be unknown. We also suppose that ll  , i.e. changes of the 4 controlling chain states occur more rarely than observed events. The sequence of instants of arriving events is observed. The problem is to estimate the parameters 1 , 2 , 1 , 2 . 3. Algorithm 1 Consider the process 1}{  ii , where 1=  iii tt  is the length of the i -th interval between arriving events in the observed flow. If the controlling chain is in the l -th state then the mean length between events is equal to l1/ . So at the first stage of our procedure we try to detect the instants of the chain transition from one state to another as the instants of change in the mean of the process 1}{  ii  using CUSUM procedures. Let the parameters 1 , 2  satisfy the condition   ;<<0 12      ,11 12 d  (1) where d  is a certain positive parameter. If the constant d  is unknown then it can be chosen as the minimal difference between the mean lengths of the intervals i  when the controlling chain is in different states that should be detected. Let n  be a lower bound of the mean number of events between switchings of the controlling chain states. For the model under consideration it means that lln  . We suppose that n is rather large (for example, n10), hence changes of the controlling chain states occur more rarely than observed events. This situation is typical for real processes such as call-center or http-server because one of the states can be interpreted as a "usual" state of the system and another state as a "peak-time" state and during each of these states several customers are supposed to arrive. Besides processes having this property are often used for simulation study of algorithms for processes with jump intensity of customer arrivals (for example, see [5, 8]).  Choose then an integer parameter 1k  describing the memory depth According to our analysis, a good choice of the parameter k  is 2nk  . The idea is to compare the values i  and ki , ki  . If there are no changes of the controlling chain state within the interval ],[ 1 iki tt   then the values i  and ki  have the identical exponential distribution with the mean 11/  or 21/ . If the chain state changes within the interval ],[ 1 iki tt   then the expectations of the values i  and ki  are different. On one hand the parameter k  should allow us to detect changes with minimal delay, on the other hand it should not be too large to contain more than one chain state change within the interval ],[ 1 iki tt  . As the initial state of the chain is unknown, we shall consider two CUSUM procedures simultaneously. The first procedure is set up to detect increase in the mean of the process and hence, decrease of the intensity, and the second procedure is set up to detect decrease in the mean and hence, increase of the intensity. For the first procedure we choose the positive parameter d  and introduce the sequence of the statistics  5  .,=(1) kiz kiii    (2) For the second procedure we introduce the sequence of the statistics   .,=(2) kiz ikii    (3) Consider then four hypothesis concerning the state of the controlling chain: – ),( 1 ikil ttH   – the intensity of the arrival process on the interval ],[ 1 iki tt   is constant and equal to l , 1,2=l ; – ),( 1,  ikiml ttH  – the intensity of the arrival process on the interval ],[ 1 iki tt  changed once from l  to m , 1=l , 2=m  or  2=l , 1=m ; The statistics )( jiz , {1,2}j  (2), (3) have the following properties:           .0>    1,2;=0,< 0;>    1,2;=0,< 2,1(2)(2) 1,2(1)(1)HzElHzEHzElHzEiliili  (4) So the means of statistics (2), (3) change from negative to positive values when the intensity of the process changes. These properties determine the construction of the procedures. We introduce positive constants 1h  and 2h  as the procedures thresholds and construct the cumulative sums (1)iS  and (2)iS  which are recalculated at the instants it . It is defined as follows   .0,=;>},{0,max=2,1    ;=)()()()(1)()(0llilililililhSifSkizSSlS  (5) Reaching the threshold lh  by the sum )(liS  results in a decision considering the parameters changes; 1l  indicates decision on increase of the mean length of the interval between the events and decrease of the process intensity, 2l  indicates the opposite. Let the sequence   0)(  mlm  be the sequence of the instants when the cumulative sum in the l -th procedure reaches the threshold lh , i.e.   0;=)(0l     .:>min= )()( 1)( lljlmjlm hSt    (6) Consider a sequence   0)( ilin  associated with the sequence   0)(  mlm  as follows   0;=)(0ln     .0=0,>:max= )( 1)()()( ljljlmjlm SStn   (7) Thus the instant )(lmn  is the first instant when the cumulative sum becomes positive to reach then the threshold. We consider the instants (1)in  )((2)in  as the estimators for the instants when the mean length of the interval between the events increases (decreases).  When implementing the procedure it is possible to encounter false alarm situations. We shall record all the exceeding the thresholds by either first or the 6 second cumulative sum. If the same sum reaches threshold several times in a row, we only record the first occurrence. Thus the procedure for estimation of instants of intensity switching is described as follows. Calculate two cumulative sums given by equations (5). Then construct the sequences  )(lm ,  )(lmn  defined by equations (6), (7). Let (2)1(1)1 < nn , then the initial value of the intensity is equal to 1 . Define the sequence      0.,>:min= 0;,>:min=0;=12(2)(2)222(1)(1)120lqnnqlqnnqqliilliil  (8) The values ,, 21 qq  are calculated using formula (8) while it is possible. If         =>:=>: 12(2)(1)2(1)(2) liilii qnnqnn  then we set Nq l =12   ( Nq l =22  ), where N  is the instant of the last occurrence. Here the odd instants 12 lq  are the estimators of the instants when the intensity changes from 1  to 2 , and the even instants 22 lq  are the estimators of the instants when the intensity changes from 2  to 1 . Define estimators for the parameters 1 , 2    ,=ˆ,=ˆ 222111 TNTN   (9) where 1N  is the total number of events occurred at the intervals ],[ 122 ll qq , Nq l 12  and 1T  is the total length of these intervals; 2N  is the total number of events occurred at the intervals ],[ 2212  ll qq , Nq l 22  and 2T  is the total length of these intervals; 0l . Define estimators for the parameters 1 , 2    ,=ˆ,=ˆ 222111 TLTL   (10) where 1L  is the total number of the switching points Nq l 12 , 2L  is the total number of the switching points Nq l 22 , 0l . The parameters   and ih  affect the characteristics of the CUSUM procedure, i.e., the mean delay and the mean time between false alarms (see [10]). If there are no additional conditions then the procedure is considered to be optimal when the probabilities of the false detection and the skip of the change are equal. It can be guaranteed by the choise of the parameter   as 2d . If the memory depth is equal to k then the sum )(liS  has to reach the threshold lh  in not more then k steps (while 0)1( lEz ). It can be provided by the choise of lh  as .4ndkhl   Note that the parameter lh  should not be significantly less than its upper bound because it can increase the number of false alarms. 4. Algorithm 2 Let we have a certain period of observation ][0,T  and N  is the number of 7 occurrences at the interval. First, we calculate the mean of the length between occurrences using the usual formula   .=ˆ NT  (11) The value ˆ  exceeds the mean length of the interval i  when the controlling chain is in the first state, and vice versa, the mean length of the interval i  exceeds the value ˆ  when the controlling chain is in the second state, i.e.   .1<ˆ<1 21  E  (12) We introduce the sequence of the statistics   .ˆ=(1)  iiz  .ˆ=(2)  iiz  (13) Using in Algorithm 1 statistics (13)  instead of (2), (3) we obtain Algorithm 2. The parameter lh  can be chosen as 4ndhl  , i.e. as at the first algorithm.  5. Numerical simulation The model for the considered flow and the suggested algorithms was implemented with varying parameters. The results of the simulation for Algorithm 2 are presented in the table below.   T 1  2  1  2  1h  2h  ˆ  1ˆ  2ˆ  1ˆ  2ˆ  1000 5 1 0,3 0,2 0,5 0,5 0,3883 5,2336 1,2276 0,3220 0,1732 1000 5 1 0,3 0,2 0,8 0,8 0,3929 5,0591 1,2370 0,2573 0,1322 1000 5 1 0,3 0,2 1 1 0,4355 4,7475 1,1972 0,2587 0,1144 1000 5 1 0,1 0,2 0,5 0,5 0,2668 5,6804 2,0053 0,2912 0,2604 1000 5 1 0,1 0,2 0,8 0,8 0,2924 5,1544 1,6825 0,1180 0,1880 1000 5 1 0,1 0,2 1 1 0,2501 5,2498 2,6283 0,1207 0,1297 1000 5 2 0,1 0,2 0,5 0,5 0,2351 6,1085 2,8854 0,3632 0,2656 1000 5 2 0,1 0,2 0,8 0,8 0,2564 5,1785 2,8092 0,1652 0,1289 1000 5 2 0,1 0,2 1 1 0,2486 5,1949 2,8831 0,1219 0,1162 10000 5 1 0,3 0,2 0,8 0,8 0,3830 4,8379 1,3439 0,2316 0,1318 10000 5 1 0,3 0,2 1 1 0,3766 4,6783 1,4326 0,1917 0,1157  First, the quality of the proposed algorithms on the threshold parameters lh  was studied. Increasing of lh  leads to decreasing of probability for the cumulative sums to reach the thresholds and hence an intensity change can be undetected. It causes increasing of error of the estimators iˆ  because of not correct estimation of the controlling chain current state. On the other hand, increasing of lh  leads to decreasing the total number of false alarms. These theoretical conclusions are supported by the simulation results. As the thresholds increase the estimators of the switching parameters lˆ  decrease because less switching points are detected on the first stage of the procedures. In the Table for 1== 21 hh  one can see that the estimators lˆ  considerably less the real 8 values of the parameters l . The best results are obtained for 0,8== 21 hh  for all intensity parameter values. Thus, choice of the algorithm parameters is a rather difficult problem requiring further theoretical investigations. Increasing of the simulation time from 1000 to 10000 does not influence significantly the estimators quality. This result stress the fact that the proposed algorithms can be used for a small sample size. Conclusion МАРs are used as models for real processes, particularly, for call-centers or http-server customers (see [3], [4]), healthcare systems (see [5]), etc. Input flow intensity estimation and pertinent model setup is necessary to develop dispatch-ing rule, to calculate optimal number of servers, etc. The suggested algorithms do not need the distribution function of the observing flow and, hence, can be applied to parameter estimation of other types of flows.   This paper is supported by The National Research Tomsk State University Academic D.I. Mendeleev Fund Program (NU 8.1.55.2015 L) in 2014–2015.  References 1. Neuts, M. F. A versatile Markovian point process. // J. Appl. Probab. 1979. No. 16. P. 764–774.  2. Lucantoni D. M., Meier-Hellstern K. S., Neuts M. F. A single server queue with server vacations and a class of non-renewal arrival processes // Adv. Appl. Prob., 1990. Vol. 22. P. 676–705.  3. Asmussen S. Phase-type distributions and related point processes: Fitting and recent advances // Chakravarthy and Alfa (eds.), Matrix-analytic methods in stochastic models, NY: Marcel Dekker. Lect. Notes Pure Appl. Math. 1997. Vol. 183. P. 137–149.  4. Breuer, L., Kume, A. An EM algorithm for Markovian arrival processes observed at discrete times // Measurement, Modelling, and Evaluation of Computing Systems and Dependability and Fault Tolerance Lecture Notes in Computer Science. 2010. Vol. 5987. P. 242–258.  5. Okamura H., Dohi T., Trivedi K. S. Markovian Arrival Process Parameter Estimation With Group Data // IEEE/ACM Transaction on Networking, 2009. Vol. 17, No. 4, P. 1326–1340.   6. Gerhardt I., Nelson B. L. On capturing dependence in point processes: Matching moments and other techniques // Technical report, Northwestern univ. 2009. 7. Duffie D., Glynn P. Estimation of Continuous-Time Markov Processes Sampled at Random Time Intervals // Econometrica. 2004. Vol. 72. P. 1773–1808.  8. Gortsev A. M., Zuevich V. L. Optimal estimation of parameters of an asynchronous doubly stochastic flow of events with arbitrary number of the states. // Tomsk State University Journal of Control and Computer Science. 2011. Vol. 4(17). P. 25–40.  9. Page E. S. Continuous inspection schemes // Biometrica. 1956. Vol. 42, No. 1. P. 100–115.  10. Lorden G. Procedures for reacting to a change in distribution // Annals. Math. Statist. 1971. No. 42. P. 1897–1971.    

CUSUM algorithms for parameter estimation in queueing systems with jump intensity of the arrival process

CUSUM Algorithms for Parameter Estimation in Queueing Systems with Jump Intensity of the Arrival Process

http://vital.lib.tsu.ru/vital/access/services/Download/vtls:000520206/SOURCE1

CUSUM algorithms for parameter estimation in queueing systems with jump intensity of the arrival process

Abstract

Similar works

Full text

Available Versions

Tomsk State University Repository

Tomsk State University Repository