Filtering of Markov Processes in Discrete Time by Monte-Carlo Method

CopniRht © 1f..\C Stochastir Co ntrol \ ·illlius. Lithuanian SSR. L·SS R. 19 ~ 1i

FILTERING OF MARKOV PROCESSES IN DISCRETE TIME BY MONTE-CARLO METHOD V. B. Svetnik, V. S. Zaritsky and V. S. Posdnikin Lpllillgrali Agrirl//II/ra/ f l/.llill/I", Ll'llillgmli. L'SSR

Abstract. For approximate calculation of optimal estimation of non-observable component of partially-observable Markov proc·ess Mo nte-Carlo method is used. The advantages of this method are the simplicity of computer programming of multidimensi onal problems, the possibility to evaluate the accurac y of approximati o n while solving the problem, universality. The use of the method for filtering processes both of determinated and random structures is considered. The "disorder" problem serves as an example. The application of the method for other estimation problems is discussed. Keywords. Nonlinear filtering; discrete systems; Markov processes; numerical methods; Monte-Carlo method; modelling.

INTR ODUCTION xiII!ation ini tial model 01' the process:X: t by model close to it and for which integration in (3) - (6) may be done explicitly. Various modifications of procedures of gaussian approximation: extended Kalman filter, se cond order filters, multiple linearisation filters and others all belong to this kind of methods. A detailed survey of these procedures one c an find in the work of Dmitriev and Shime levich (1977). Di Mazi and Runggaldier (1982, 1984) used approximation of the initial and transient dens ity function s by the sequence of conver g ent to them functions for each of which integrals in (3) - (6) are calculated explicitlv. In Kushner (1985) Markov process is approximated by a Markov chain with a finite set 01' states, that's why integrati on in (3) - (6) i s re p lac ed by summation.

Consider a partially observable Markov process X1;€Rn.+·'U'1;E:R""),-t;€R s in discrete time i-O,t, •••• The problem of optimal filtering is to calculate for each conditional expectation.

J1'\,t=Mt"'t/s!}

(1)

Where m t € Rn. is the vector of optimal estimations of nonobservable component

~

s~= (so' . . . , ~t) is the sample of observable component. The accuracy of fil tering is determined by n, x n,. covariance matrix of estimation errors

rt=Mt(~-mt)(~-t-J1'\,t)/S:}. (2) There are many works on the solving the optimal nonlinear filtering problem,includ ing funda mental works of R.L.Stratonovitch, R. Sh. Liptser and A.N. Shiryaev, G. Kushner and others. As mt

=

{n't}t ftC""t It:

)dV-t ,

0)

rt=SR"("'t-n'\)('l\-tnt)Tft("'-t/.s~ )d~t,

(4) The considered methods posse s not only und~~bted advanta ges , .but also.a num~er of olsadvantages:.a hlgh com~l~cacY . ln computer pro grammlng o~ ~u~tldlmenslonal pr o blems, ~ack.of p osslblllty.tO cont:ol the approxlmatlon accuracy whlle solvlng

then the filtering problem is reduced to calculating a posteriori density function and the further n-dimensional integration to obtain m-t;, r~ . For density function we have a well known difference equation 1

where

ft "' 1 C1Jt+Js:+ ) ={Jlj}t/5:).ft"" /t('l\..,,5t+J~t,5Jd:&t/~n.(n.u.merator) d.1Jt + 1 , . f (1J t j 1J- t \ the problem and to change thlS t+1/ t

tt-I'

t+1

(5) accuracy, nonuniversality at solving different problems ana consequently the difficulty in making universal computer programs. In a sense, method of filterin6 based on calculation of integrals by Monte-Carlo method, has no such disadvantages. Apparantly, such solution was first considered by Handshin (1970), then by Zaritsky, Svetnik and f himelevich (1975~1976), by Dmitri~v and Shimelevich (1977), by Svethnik (1984, 1985). The main idea of the present work is the following. For the vector of optimal estimations J1'\,t and the covariance matrix

t't/

is transient density function of the process xt-. The initial condition to (5) is the density function \=f(~~)(S f(19::~)d~)-1 (6) o 0 ~o/ 0 o'~O RI\, 0 0' 0 0, where fo(~0,5o) is initial density function o f the process x t . The explicit calculation of integrals (3) - (6) is possible in particular cases only, for example, when density functions fo and tt+1/t are gaussian, that's why to find J1'\,t and ~ one has to use numerical methods. As a rule, known methods are based on appro-

rC-&/l:

\:; I

\', B, S\t:tnik. \', S, Zaritsh and \ ', S. I'osc\nikin ~ their representations as ratios of functionals expectations are true. (Here pro ofs are not given). Functionals are given on the trajectories of an introduceCl in a special way auxiliary Markov process (see (9), (10) in section 1). It is the replacement of exact values of expectations in these representations by their sample mean calculated on simulated samples of the auxiliary process that gives formulae for the approximate calculation of m t and ft, In section 1 this approach is used in the case, when the process x t submits to the system of difference stochastic equations with a determinated structure. The last term means that in the proce s s of estimation the model (the ri ght parts of the equations (7» doesn 't change or changes in a determinated way. In section 2 :.lonte-Carlo method is applied t o the equation with random structure (20) (Kasakov (1977), Svetnik (1984». Specifically, this method can be used for the so called "disorder" problem (Kliegene and Telksnis (1983». I<:onte-Carlo method allowing in this case not only to calculate a posteri ori probability of disorder but also to estimate the moment of its appearance. In section 3 MonteCarlo method i s used to s olve problems of interpolation, filtering with a tolerance control (Dmitriev and Zaritsky (1974) and tolerance control (Zaritsky, Svetnik and Shimelevich (1 976), Svetnik

(1 985) •

-1

It+l= a.(lt,1 t ,t)+ H(lV5t,t)(C(lt,St,t) x x (§,t+1- A(cV5t,i»)+FClt,St;t) W(t+1) ;

(ll)

Q)(St+l,St,lt,t)=(5 H1 - Ml.t,St,t'Y x x (C(i"t' St,t)

r (5

t +1 - A(lt' '§,V 1::);

.ri+t p'\;(d.et(CC~t,St,t»)) xQ,)

(St+l' ~

t'

-1/2

EXp(-O.5 X

(12)

Zt,i)).

At t =0 Zo is a rando m vect or with density function

Fo

"0

:fo(lo)=fo('lJo)

\'lJ-o-l

o.

( 13)

is the function of the random vector and sample of observation So :

Po = .f0 ( S %

0 ) \

~

o

__

~

0

(14 )

Here " " is redenot at i on of argument in (13) and repl a cement of a non-random argument in (14) by a random initial vector co . That's why to(zo) has the sense of densit y function, and 9 0 is a random value. In (ll) uJ (t) E: Rnis a gaussian process indepe ndent with Ze. for all l=O , I, ... ,i-l; H=6,6; +6l6~; 1 p= (t\6; + 6z.6~- HC- W )'/2.. In (13) and (14) corresponJingly

t(~o)= SfC'OQ ,n)dh; Ifo

~(1o~)=t(~'So)/toC'd"o)· (15)

The auxiliary process Zt and functional of it P't have "a posteriori" char acter , because sample of observations so'~ l ' ••• ) 1. THE PROCESS ,iITH DETERi· INATED ~'t is one of inputs of the ri ght parts STRUCTURE (11) and (12) as well a5 in the initial condition (14). Let the process x t submit at t=1,2, ... The procedure of calculation of m t and to the system of difference sto c hastic r't consists in the followin g . According equations to the equation (11) and the initial ~+;' C1(Xt ,t)+ £1(Xt ,t)c1(t+l) +{ilx:t ,t)t-2.(t +1) , £'t+l = A(xvt )+B/x t ' t)E 1Ct+1) +l\(xt,t) £2,(t+1), (7) where a., A,6 1 ,u"B"El2, (ar guments are omitted) are correspondingly ~xl,5x1, n-xm, rt-xp , 5x m, 5'P determinated vector - and matrix functions; t. 1 (t) Eo R''', c 2.(t) E RP are gaussian processes independe nt with Xl, for all e :0,1, . . . ,t - 1;

M{t1(t)}=-M {c/t)} = M{E./t)E;Ct)}=O; MtE./t)(Cr)r= 100t't; M{E}t)E.:C't)}= le,\>,"; "C'=1,2, ... ;

condition (13) be ginning from t=O and further on for t = 1,2, . . . samples l.t[o(..], 0(.. = Cl.. of the process l.t are simulated in parallel. On these samples according to (12) and (14) samples Pt[~] of the funct ional Ft are calculated. The obtained samples are substituted then into the approximate formul a e L

mtc;::

L

n\,=(fl lfol] 1\[01..])/];; Pt[ol]

(16)

L L -1 I is a unit matrix of the corresponding r~ r=(L:(l [ 0 (8) of expectations in (9) and (10) by sample for all =et and i . At t=O Xo is a random mean. Let ~t[~], m t [ l], mt [ L] be vector with a densit y function fo(x o )' el,!;ments of vector lot, rn t , m. t an~ rt(L,j) It can be proved that for mt and rt (L ,J) be elements of matrices rt ' r t ; formulae j, L- ~, Then with the a~curacy to the m t = M {~t.J\}( M{ J\}f1, (9)

r

= M {Ci\-tnt)(lt-mt)Tpt}(Ml9tn-1 (10) t are true. Here Zt E: Rn. is an auxiliary members of the order L-~ the following f.\a rkov process, and Pt E: R is a funcformulae, characterizing the method accutional of it submitted to the system of racy take place difference stochastic equations

153

Filtering of Markm Processes

(18)

(9)

If in the right parts of this formulae 1'l1--t(L) _ ~s replaced by nt,-t.(i.) , Ft(~,P by rt (~,d) and expectations by sample mean then we shall obtain approximate formulae, which can be used in the solving of filtering problem. In a particular case (Handshin, 1970) in the equation (7) Ba =0, 1\ = 0 and functions CL(l\,t), 81 ({}-t,t) do not depend on s-t. That is why is a Markov process and 5t has a sense of measurement of function A(X-t, -t) with an error B2.(X-t,t)E.2.(t+l). In this case equation (11) coincides with the first equation (7) and modelling of the auxiliary procese is done consequently according to a priori description.

"'-t;

Lt l.-t i.t "t+t=CL (X t ,t)+i>, (Xt ,t)E/t+l) + 62. (xt ,t)E-/t+l),

t

= l/lt =K,X-,v+1 ,X}). t

t+1

..P-t+1

~ 5\ (dd (Cj~(lV5i>,t )f1l2"

"EXP (-0.5

cph(~t+1,.5i;,~t,t));

(23)

cpJ (S1o+1' E,t, li-,t)= (~tH-A.tt(~t,5 t , t) )Tx t

1

"(C4t Cl t ,St,t)r x (St+1- Ait(lt,~t,t )). Here it= 1,l is a process with the same matrix as that of the process ~t (24)

V<

of transient The process ·tJJ (t) in addition to the innumerated characteristics does not dependonjt. atf=O?l, .... ,1;-I; HJt=O~t(B~t)T + 6~1o(I)~tf;

pjt=(6:t(6~-t;)T +5~t(6!t)T_ Hjt(C~t)hjt)T)1/: At t= O

. ;:'0 is a random vector with ~ ~ ;t)+B 1 (x.-t,t)C f (t+l) + Bl. (::ct ,t)E../t+1). t

~ tt+fA (X

~t 6 Lt a't '''' .... Here CL, "o2,B" Bz . are elements with number ~t of the given family of vector and matrix functions.In addition to the above mentioned properties of t.,(t) and t-2.(t) we assume their independence with Lt for all .{. = 0,1, .. . , t -1 • We assume as well that for. all Lt the. condition (d) r~corded for B~t(:X:-t.,t) Bz... (:c.t,t) and c't>(X-t,-t)takes place. Values of the process Lt have the sense of numbers of the descriptions (structures) in accordance to which process x,t is evolving in time. Random transitions of x t from one value to another reflect random structure changes (Kasakov,1 977). Let determine process [.1; for t = 1,2., . " by the matriX of conditional transient probabilities

,X)== (Pr{l

4t

"'t -z ... conditional probabilities.

Let (X t ' l.-t) = ( trt , S t ' Lt) be a Markov process of a mixed type, in which L= 1.I Process x,t for t = 1,2., . . . submits to the system of difference stochastic equations

t~l

xC}' t+1 -A Clt' S t't))+f'h(i~ 1; , ~ ·V t)0(t+1)",

:Pt+1/1;(x10+1 /Xt ) \ ,.Vt+1- l

2. THE PROCESS WITH A RANDOM STRUCTURE

'Y-\:+1 / +(X

One can prove that with the assumptions made for the process x t formulae (9)and (10) are true for m t and rt • But the auxiliary process ~-t and functional 9t are submitted at t = 1,2., .• . by a system of difference equations with a random structure h 4 4 )~ It,,=a {It'.~t;t)+H (It,S1o,tXC (It,St,t) x (22)

(21)

Here P~(/) is a conditional probability. For -\;=0 the dens i t y function J-o(Xo) of a random veotor Xo and 1 ><1 vector of initial conditional probabilities are given. Elements of matrix J\+-I It and vector Po determine in probability sense the causes of random structures changes. So if Pr(~=L/::C J= PrLi.o=L\; P,,[Lt+1=L/Lt=K,X-t+1 ,xt l= loo I =Pl'l'~ t'<=i./i. t =K} then structure changes will be initiated by an independent Markov chain. On the contrary if p~ t Lt+1 = l/i.t=K,X t +1' x'\;\ = P,,{ ~t+1 / X H1 } then the structure changes will be determined by the process ::et'

the

(20) density function (13); 9Q is a random variable (14); to - is a random variable with probabilities equal to the corresponding vector elements ;Po(lo'~o) = 3\xo) l~o=lo (25) The left parts of (24) and (25) ~O, ... ,Si; are samples of observations. The procedure of apprOXi mate calculation ef In-\> and r t is the following. By the equation (23), the initial condition (13) in accordance with the initial and transient probaDilities (25), (24) samples Z t (0<. J and it (o<..J of processes Zt and Jt are modelled in parallel. The sequence of modelling is the followin g : l,ioi..] -- Jt[o(.] ..... It[o(.] - - h"'1[o<.] - . . . On these samples by (23) and (14) samplesFt(o<.] are calculated. Then the obtained samples are substituted to (16)~ and (17) by which approximate values m t and r t are calculated. The approximate calculation of a posteriori probabilities (26) used in a number of problems for drawing a conclusion on the state of the system, submitting to the equation (20) is of interest to the processes with random structure.

(27)

\' . B. S,elllik, \ '. S. Zaritsb' and \ '. S. Posdnikin

154

where j ·Ch) = 1 i f h = i. and iL(A)= 0 i f ho# ~ ". The accuracy of calculation according t o (26) is determined by the formula (we i gnore members of the order L- 2 )

M{Crr_/~)-nt(~))2/.s~} =L- 1 " "M { (j ~

(J

t )

-1C -t ( l ) /- F~

(28)

} ( M {Ft })~

It is often required fo r the process with random structure not only to evaluate the structure's number but also the moment of the first appearance of the structure with such a number. ~e s hall explain it on the exa ~ ple of the s o called problem of " disorder" (Klie ge ne and Telksnis

(1983».

Let pro cess ~t = 1 ,2. at t = 0 have s tructure 1 and then in a r ando m moment A=1,2 •.. _,X change it for structur e 2 in wh i ch it stays furth er on during the entire ti me of observation . As s ume that there is a sample to, , , . , t:K of fix e d number X of observations . It i s required firstl y to calculate a poster iori pro babili ty

Tt:x(~) = PI' {~x = 2/s~ },

3. INTERPOLATION AND TOLERANCE CONTROL

Let it be required to calculate in case of random or determinated structures ,X

Int ,x. = M {t\l1 o } ;

rt.X fCt

•

= M {(1J- i

x =Pt'Ft =

m

t,X

)(1J- t

m

T

t.'x

) 1'1:

X

:00

}. '

~ lg~1 .

These characteristics are analogues of (1), (2) and (26) in case of interpolation. One can use for their approximate calculation f ormulae L..

-1

)

"", int."" ",,=(2:. ~-t[oCJ.Pt[oCJ ill-x 0£. =1

In t.X

rt,X~rt,J\-~/=ct(i~.L[9(..]-m-,.J(~t~J-rr\ .""jpJC[o(.~~ 0(.=1 " "',,,," -1

L

at-t X CL)=

,

(F-, .t ~ (J 10[0(.] p:J{.[o(,] )

%$L

which in the ~ i v en pr obl em is a posteri ori probability of " disorder " during X times, and secondly if nx (2. ) = 0 to estimate the moment of " diso r der" that is to find

$'JC = ~ Fx [o(..J

.

Formulae (18), (1 9) and (28 ) are true for calculation accuracy analysis if in them Ft is replaced by jJ'X and index t of m, r , :It , in- , r,:it by tte do uble index -t , X.

Dmitriev and Zaritsky (1974) considered the problem of calculation of estimations Q

m1;

=

M {,\

t t "t+1 Iso, -C-o Eo l::lt:. },

where where 0 < }J- (X t +., J:)t ) <:. 1 is a probability of t ransit i on fr om str uctur e 1 t o structure 2 . The soluti on of the "di so r der" probl em by Monte- Carl o method is reduce d to c a lculation_by (27) at ~=2. appr oxima te val ue TeJ( (2) and , i f ar,,(2) 0# 0 t o the applicati on of the appr oximate form ula

where 71[::<. ] is the moment of th e first t r ansi ti on h L,{..] fr om 1 to 2; (X) means summation only fo r those ~ for which A [oI..J ~:Ko ' If it's requir ed to find a post eri ori di stributi on of the "di sor der" moment one can use the ap proximate fo rmula

Pt' L (el

l;\,= {. 15~. ~x=2.} ~ L (X)

~(~, .Px[oC])(~

-1

9x l9LJ)

,t=1,:K..

where (t) means summation only f or those 0(,. fo r whi ch P:K [cl.] = e.. Mark, that the ca l c ulati on of a posteriori characteristics of " disorder" is do ne in a universal way, because the mode lling pr ocedure of processes ~ t , ~t is the same for all characteristics.

t +1 '1>+1

~

Q =Q~ ...)(Q;

Q

n.

t

eR ;~=(~'''''~).

Q

The estimation ~t is the solution of the optima l filt ering problem in the presence of add iti onal data tha t '\J-~ Q'l;' for all 't' = ~ For the approximate calculation of estimation m,Q. it' s posible to use formula ~

where (~) means summati on onl y fo r those 0/... fo r which i!:.~[o(JE:Q , It's often of intetest t o determine a posteriory probability

(Zar itsky ,Svetnik and Shimelevich (1976) . For calculating by Mo nte-Carlo method one can use formula

In all the considered problem the procedure of mode lling of proc e ss es Z t , ~ t and calculating of Pt remains the same.

Filtering of \Iarkm' Processes

CONCL USI ON S

On the basis of formulae of the (9) and (10) type for the solving of filtering problems a procedure of approximate cal~ulation of . . rn t , r t , :flt(L) by MonteCarlo methoa lS suggested. The main part of the procedure is the modelling of auxiliary processes L t ,jt and calculating on their samples the samples of the functiona19tO Formulae for the accuracy analysis of approximate calculations were obtained. These formulae may be used in the process of solving problems. The universality of the suggested approach to the solving of filtering problems and some other problems is manifested in the fact that modelling part remains the same (This is clearly seen in the example of the problem of "disorder"). Standard modellin5 programs of random processes available at present make computer programming by using Monte-Carlo method quite si~ple, and parallel character of modelling procedure will, apparantly, be the cause of its effective realization for parallel processors. However, one should mind the low accuracy of the method, the order of which is, roughly speaking t... -112 , which requires a great number of simulated samples and consequently the necessity to use highspeed computers or to apply one or another method of accuracy raising, but t~ey in their turn may essentially compllcate computer programmin!:':.

REFERENCES Di Masi, G.B. and Runggaldier, W.J.(1982). Appoximations and bounds for discreteti :.. e nonlinear filtering. Lect. Notes Contr. and Inf. Sci, 44. Springer Verlag, Berlin, pp. 191-202. Di Masi, G.B. and Runggaldier, W.J.(1984). Approximations for discrete-time partially observable stochastic control problems. Lect. Notes Contr. and Inf. Sci, 61. Springer-Verlag, Berlin, pp. 19·1-202. Dmitriev, S.P. and Zaritsky, V.S. (1974). Optimal proceSSing of measurements with using the results of tolerance control. Avtomatica i telemechanica g, 31-36 (in Russian). ' Dmitriev, S.P. and Shimelevich, L.I.(1977). Konlinear problems of navigation data rocessin. CNII "Rumb", Leningrad in Russian). Handshin, I.E. (1970). r.lonte-Carlo techniques for prediction and filtering of nonlinear stochastic processes. Automatica, 6, 555-563. Kasakov, 1.1'.,. (1'9"77). Statistical dynamics of systems with changing structure. Nauka, Moscow (in Russian). Kliegene, H. and Telksnys, L. (1983). 1.1ethods of the detection of change points in the properties of random processes. Avtomatica i telemechanica, 10, 5-56 (in Russian). Kushner, H.J. (1985) Probability methods for approximations in stochastic control and for elliptic equations. Nauka, Moscow (in rtussian).

155

Svetnik, V.B. (1984). Optimal estimation of phase coordinates of systems with random structures. Using of microcomputer and microprocessors for technological processes automatization. Leningradsky selskochoziastvenniY institut, Leningrad, pp. 35-42 (in Russian). Svetnik, V,B. (1985). Using of Monte-Carlo metnod for computering of multydimensional integrals taking place in nonlinear filtering problem YII Allunion conference "Monte-Carlo method in computing mathematics and mathematical phisics". Sibirian otdelenie AN SSSR, Novosibirsk, pp. 113-116 (in Russian). Zaritsky, V,S., Svetnik, V.B. and Shimelevich, L.I. (1975). The Monte-Carlo method in problems of optimal data processing. Avtomatica i telemechani~, 12, 95-103 (in Russian). Zaritsky, V.S., Svetnik, V.B. and Shimelevich, L.I. The calculation of a posteriori probabilities for decision rules of tolerance control. Izv. AN SBSR.Tethnicheskay Kibernetika: 1, 209-215 (in Russian).