Pseudo-random numbers for constructing discrete Markov chains by the Monte Carlo method

PSEUDO-RANDOMNUMBERSFORCONSTRU~T~G~ISCRETE MARKOVCHAINSBYTHEMONTECARLOMETHOD" I.M.SOBOL' Moscow (Received 14 LJecember 1972)

A METHOD of mode&g on the basis of returns is described for constructing discrete Markov chains; the method enables the coordinates of certain points of the infinite-dimensional unit cube to be used as pseudo-random numbers. In certain problems our method for constructing chains accelerates the convergence of the mean to the stationary values. rn~oducti~~. Consider a homogeneous Markov chain

with a finite or countable number of states oj, /=I, transition matrix

2,....We denote by (nil) the stochastic

Pij=P {l;k+i=wq&c=Oi), and assume that a stationary (final) distribution normalized solution of the system of equations

i, j-l,

2,. . . ,

c-9

pp, exists, representing the unique non-negative

We further assume that the initial state Ei =oio, is f=ed, recurring and positive. Finally, we assume that the functions f (oi) and g(oi, o’) are such that the mathematical expectations exist with respect to the stationary distribution M*f = 2 j (oi) pi*, *=1

M*g =

2 g (mit 6) pi*&. i, j==l

(4)

In certain problems we wish to evaluate the quantities (4) when the stationary distribution pj*. is unknown. For this, we can construct a sufficiently long chain (1): since the laws of large numbers hold for Markov chains [ 1-5 3 , for large N we have

The method of inverse functions is often used to construct the chain (1); in the present case the algorithm of this method is extremely simple. We divide the interval OGzc:(l into intervals *Zh. vjThis1, Mat. mat. Fiz., 14, 1, 36-44, 1974.

36

37

Constructing discrete Markov chains by the Monte Carlo method

Assume that jk=wi is known. To find Ek.A,, we choose a random number yk (i.e., the value of a random variable y, uniformly distributed in the interval (&I)) and check which of the intervals (6) contains the point s--y,: if

YkEAij,

then

!$+,=6+.

(7)

It was shown in f4] that, for chains (1) with a finite number of states, if we use, instead of yt, + ’ * , n;ik,. . . in (7), numbers y,, . . . , y,, . . . , forming a completely uniformly distributed

sequence, the relations (5) will hold in the “classical” sense: the convergence in probability is replaced by ordinary convergence. In the present paper the chain (1) will be constructed by using the coordinates of the basepoints X, of certain in~nite-Dimensions quadrature formulae 163 of the type

where X = (%, x2,. . .>, m

K”={O~a,~l;n.=1,2

,...

>,

This method may be called modelling on the basis of returns. Modelling on the basis of returns.Instead of y,, y2, . . . in (7), we shall use the coordinates of apoint X,=(5,,, xl.‘, . . .) , but only until the chain returns to the state wiO. If gr+, =CB’O,we have to use, instead of “fr+f, ~~~2,. . . the coordinates of a second point Xz= (G_~,xz2,. . .) until the chain returns a second time to the state w’0, etc. The chains thus constructed (given a suitable choice of the base-points X,) also enable us to evaluate M*f and M*g from (5).

Notice two features of the method. Firstly, with this method of construction, the chain (1) proves to be not “completely ergodic” (see [4]). For instance, if the function h(w’, e$, 0’) is such that the mathematical expectation

exists, then the relation

(9) will not in general be satisfied for the chain. This situation can be corrected, however: in order to ensure that relations of the type (9) hold, we can introduce modelling on the basis of composite returns (see below). The convergence when modeiling on the basis of returns can then prove to be faster than the convergence in probab~ity in (5).

I. M. Sob&

38

Secondly, it would seem that there are still no “good” (from the practical point of view) completely uniformly distributed sequences of numbers. But there are “good” quadrature formulae of the type (8): as the X, one can use the points of generalized LPr and Holton sequences [6]. Mathematical expectation of sums over a cycle. We shall describe the piece of the chain (1)

until the next return as a simple cycle. For instance, the first simple cycle is g,gkfmiO for k = 2, . . . , T, but &+,=co~@.

. . . -&,

where

Lemma I

Let all the conditions listed at the start of the introduction be satisfied. Consider the random variables cDf and Y 6, dependent on the simple cycle El- . . . -t !&:

cf T

m’r =

5

Y, =

&kh

2 g(Ek,

Ek+l)-

ti=i

K=i

For the mathematical expectations of these variables we have MY,/Mr=M*g.

M(D,/Ma-=M*f,

(11)

The proof will be omitted, since relations similar to (11) are used in the Kolmogorov-Doeblin method [2, 31. We merely remark that the laws of large numbers, corresponding to (5), are easily derived from Eqs. (11). In particular, if we consider a piece of the chain (l), consisting of m simple cycles, i.e., N=N,=rli. . . +L, then, as rn-+-,

Since here

k=l

S=l

we find, on dividing one of these relations by the second, that N* +-zf

(ck$

$$= M*f.

Reduction to injkite-dimensional integrals. Let the point X-

(xi, x2, . . .) d?. We shall agree to write XEIIiocrfar . . . atiO,if x1EAioa,, xz~A=,~,, . . . , ~t+t~&, while the xt+z, x~+s, - - are arbitrary. Consider the sum of such parallelepipeds over all possible or, . . . , at, different from iO:

It can easily be seen that if, instead of rl, rr, . . . we use in Eq. (7) the coordinates x1, 52, . . . of the point X, then the length of a simple cycle El- . . . -El is equal to c if and only if X=Gt. The measure (volume) of the set Gt is easily evaluated: since the random point I’= (ri, ‘fz, . . .) is uniformly distributed in K”, we have

Constructing discrete Markov chains by the Monte Carlo method

P(I’EG~)

VGt =

39

= P(z = t).

Consequently, (13)

and OD

G,=K"-

ct=1

G,

is a set of measure zero. We define in K” the function

(14) k=l

Obviously,

cs m

s ’

(D(X)dX=

Cp(X) d>Y.

(15)

f=1 Gt

Km

Since the random variable CD f is equal to 0 (I’), its mathematical expectation is equal to the integral Mm, = MQ (I‘) =

s Q(X) dX.

(16)

Km

It is easily shown in the same way that if, for X=G*, l
‘r (x) = t

b?(Ek, gkil),

we define functions

z(X)=t,

k=l

the mathematical expectations of the random variables Y 8 and T are M’& = MY (I’) =

5 Y (X)dX, K-

Mr=Mr(I’)=

s z(X)dX.

(17)

Km

If quadrature formulae of the type (8) are used to evaluate approximately the integrals (16) and (17), we obtain precisely the algorithm for modelling on the basis of returns. Unfortunately, the convergence of the quadratures for the functions (Ii (X) and Y! (X) does not follow directly from the results of [6], since these functions are not bounded. These functions are extremely simple, however: they are constant in each of the parallelepipeds IIioala, ario. It is therefore quite possible that the quadratures will usually be convergent. Finite-dimensional approximation. Convergence. Consider the chain Ei+ * - * -fb-f ..-, which is constructed according to the same rules as the Markov chain (l), starting with gi=oio, except that, if ~A#oio fork = 2,3,. . . , nfl, we set En+2=aio.In other words, if z)n, we set z=n+l. The subsequent cycles will be constructed in the same way.

40

I. M. Sobol’

The function corresponding to this chain is t

if

‘E

f (Ek),

@l(X) =

l<.f
XcGt,

k=l nf1

isf

(18)

x.

(&) at other points

k=l

Obviously, the CD,(X) depend only on the first n coordinates x1, . . . ,G, of the point X. Denote by K”={O
P=(q, ...,z,)the projection of the point X onto the n-dimensional cube k=l,

2,. . . ) n}.

Theorem

We select in K" a sequence of points Xi, X,, . . . such that, for every n, their projections P,,P:,...on to K" are uniformly distributed in K”. If 1f (cd) ( GC,then

lim lim -YL&zJn(X,,= n--r00m-m m s=1

1 @(X)dX.

(19)

Km

Pro05 Thefunction@‘,(X)=@.(P)isboundedin K",((@,(X) IGC(n+l)),andis It is therefore Riemann-integrable, and it follows constant in each of the parallelepipedsIL,,,,,. cs.tioe by Weyl’s well known theorem that

lim A- &3,(X,)=

m-m

1

S=l

cD,(P)dP = 1 @,(X)dX.

c33

Km

If”

Now consider

s

@(X)dX

-

5 CD,(X)dX

1[@LX) - @n (X)1dx-

t =n+2 Gt

Ii”

KC=

2

=

Recalling Eqs. (14) and (18), we find that

2 sI@-@,,ldXs t=n-J-z

Gt

s2

g t =n+2

Gt

k=n+2

The last sum is the remainder of the convergent series II

c

tVGt = MT

1

If&)IdX<'

2 t =n+z

t&y

41

Constructing discrete Markov chains by the Monte Carlo method

and tends to zero as n-+m. Hence Eq. (19) follows from Eq. (20) as n-+00. The theorem is proved. No matter what the point X, of K”, matrices (pij), exist such that X,=G,. For such problems, an attempt to evaluate CP(X,) leads to the program “cycling”, while the @,,(X,) can be computed without any complexities. We shall take the case when it is possible to avoid “cycling”, without introducing m,(X). > >O in the matrix (2) for all i. We renumber the states in Assume that the probabilities Piio/a such a way that i,, = 1. If X,= (zSi, . . . , zsk, . . .) and lim

zsk

=

(21)

0,

R+rn

then X,@G,. For, since it turns out that xsh’
Gl.

E z

t=*

It is clear that in the case considered, rDn(X,)=cD(X8) for all G-k’. Condition (21) is satisfied e.g., for a generalized Holton sequence. Incidentally, any points X,= (Q, Q,, . . .) in (8) can be “cut out” and replaced by points X,‘= (z*~,. . . , tsksr0, 0,. . .), where k, + 00 when s +m. If the projections of the points Xi, X2,. . . onto P are uniformly distributed in Kn for every n, the projections of the points Xi’, X2’,. . . . will have the same property. The finite-dimensional approximations Y n (X) and ‘t,,(X) are defined by just the same sort of expressions as (18). Notice that any piece of Markov chain (1) constructed in practice can be regarded as a realization of a finite-dimensional chain with fairly large n. Modelling on the basis of composite returns. We fur two initial states Ei=oio, %2=oil, of the chain (l), such that pi,i,>O. We define a composite cycle as a piece of the chain Et-+. . .-+&t, if gT’+i=O ‘0, ETr+z=oij and both these equations are satisfied simultaneously the first time (~‘2 1).

Instead of 72, Y3,. . . in (7), we shall use the coordinates of a point X1= (xii, zi2, . . .) , until the composite cycle ends. Then, instead of T%,+~,Y~~+~, . . . , we use the coordinates of a second point Xz= (x2,, xz2, . . .) , until the second composite cycle ends, etc. Lemma 2

Let ail the conditions listed at the start of the introduction be satisfied, and let M*h exist. Consider the random variable

(22) dependent on the composite cycle E,-+&expectation of this variable we have

. . . -Err,

&f+i=E,,

E,r,+~=b

For the mathematical

42

I. M. Sabol’

scheme of the proofi Consider the extended E.51Markov chain ni --f . . . -+ qlk.-+ . . , , where qk==@A,&+f ). The states Sar of this chain are the pairs (u’, ~9). A simple cycle for the extended chain is the same as a composite cycle for the chain (1). If we introduce the function V($P, Qt) = h(fJP, 09, 07) a,, with R’= to=, we), S-J’=(wx, a?), where (s8,is the Kronecker delta, @hwill be equal to Y,,, for the extended chain and by Lemma 1,

On evaluating the transition probabilities rcif=&gx~Ly and the stationary probabilities ,T~‘=P~*P~~,for the extended chain, we find that

&ph’=

r,

h’fR’,

Qf)n,‘n,t=M.*h.

r,t

Consider the sums over the composite cycle

Then,

MB,/Mz'=M*f, MB$Mz'=M*g. These equations follow from Eq. f23),, since, if fa =f(o’) Wf or M*g respectively.

or h=g(w’,

(24) oj), then M*h becomes

Equations (24) show that chains constructed by modelling on the basis of composite returns enable M*f and M*g to be evaluated as well as M*h. In general, the mathematical expectations with respect to the stationary distribution for any functions of s successive states, with s~+I-l, can be evaluated with respect to the chain (I), which is constructed on the basis of returns of order 4, if the end of a eompo~te cycle is defined by the condi~ons &+k=& for k = 1,2, . . . f 4. By means of invite-Dimensions points X1 , X2, . . . , we can thus construct chains (1) which, though not “completely ergodic”‘, have “ergodicity of arbitrary order (1”. ~~~e~caI exurnpie. We considered a chain with three states 1,2,3 and with the matrix 0.4 0.3 0.6 (Irij :=

0.2

it.2

0.6

i?.Z 0.4

0.4

for which P;= P/Bz~r/z%, ‘L/z2). We constructed three realizations &-rg,+ . , . -+&, of this chain: A, starting with .$=S on the basis of standard pseud~~ndom numbers; B, starting with 5,=3 on the basis of simple returns; and C, starting with g;-3, &=2, on the basis of composite returns. In cases B and C, the coordinates of the same points Xt , X,, . . . , forming an L&-sequence [6] , were used.

Constructing discrete Markov chains by the Monte Carlo method

In Tables l-3

we quote the empirical frequencies

43

of the different states (digits), the pairs and

triples in each of the chains A, B, and C, and also the corresponding

stationary

probabilities

(*). All

the frequencies and probabilities are multiplied by 1000. The closeness of the empirical frequencies the stationary probabilities were estimated by means of the usual x2 criterion. The corresponding probabilities

P = P(x')are quoted in the last columns of the tables.

Even with a chain as short as this, it is clear that the chains B and C enable the limiting distribution

to be approached

more rapidly than does the chain A. It is also clear that modelling

the basis of simple returns is not sufficient when we want to estimate the distribution TABLE 2

TABLE 1

170 172 173

A

B c +

)

3.X 325 334

on

of the triples.

480 503 493

0.38 0.87 0.76

15 40 11 65 10 h5

A

B C

182 1 318 t 500 1 -

l

65 60 65

1181 54164

93 63 62

85 116 191 221 173 0.24 102 97 201 196 202 0.95 98 98 204 204 194 0.92

1 6-i ’ ( 100 ( 109 1191 1200 ( 200 ( -4

TABLE 3

j lil~iIt~2+il/ A

0

B C

0 0123

L

121516 1

A B C *

?,

+/ 5

5

12 5 7

12 5 0

5 5 12

0 i 11 1 11 (

122 /

2?i /

27,?.j 2i? j

131 /

$12 j

123 j

213

IO 16 14

10 18 12

32 13 7

12 23 19

17 0 17

25 37 34

28 44 38

50 ii

13

1 13

19

1 22

1 30

( 33

1 38

11

(

(

?i31 1 223132il it?/ 331)

132 1

133 1

313 1

252 1

45 0 1:

53 39 48

34

55

91

1

On the order

50 53 50 23 37 39 34102 46 43 43 38

1 38) 401 40) 401 44

ofconvergence. The

233 /

332 1

333 /

323 1

P

1

usual method of modelling the chain (1) in which “existing”

Ye, y?, . . . , are used in (7) may be treated as modelling on the basis of returns with random points X,=I’,, Xz=rz, . . . , where all the r, are independent and uniformly distributed in K".This approach enables an expression to be obtained for the error of the approximations (5). random numbers

Let N=N,,=T~+ . . . +r, be a piece of chain, consisting of m simple cycles. If m is large, then, by the central limit theorem,

to

I. M. Sobol’

44

where the random variables cm and c,’ are approximately normal with (0,l). Recalling Eqs. (11) and (16), we can easily find that

whence it follows that, with high probability,

The variances DO! and Dr can be estimated empirically. (In the case of a finite chain (1) the expression for Dt appears in [SJ ). If, for modelling on the basis of returns, we use points X, which satisfy the conditions of our theorem, the order of convergence can prove to be better. For, if the chain (1) is finite and all the p,) are binary rational, the function CI,,, (X) will take a finite number of values and will be a finite Haar sum [6]. But then, CDn (P) ES,, and if the points X1, X2,. . . , form a generalized LP,-sequence, the convergence rate in (20) will be 0(1/m). Notes. In a sense, the best approach is to choose an initial state (I)‘~~, corresponding to maxllflum p,“*. Then,X,, X2, . . .is a minimum and we can expect that it will be possible to manage with a finite-dimensional approximation with fewer than n dimensions. The replacement of evaluation of M@ (I’) by the evaluation of MO n (I?) implies the approximation of the Monte Carlo algorithm with infinite constructive dimensionality (c.d. = w) by an algorithm with c.d. = n (see [7]). Translated by D. E. Brown REFERENCES 1.

KOLMOGOROV, A. N., Markov chains with a countable number of possible states, By&Z. MGU 1,3, 1-6, 1937.

2.

KOLMOGOROV, A. N., A local limit theorem for classical Markov chains, Izv. Akad. Nuuk SSSR, Ser. matem, 13,4, 281-300, 1949.

3.

SARYMSAKOV, T. A., Foundations of the theory of Markov processes (Osnovy teorii protsessov Markova), Gostekhizdat, Moscow, 1954.

4.

CHENTSOV, N. N., Pseudo-random numbers for modelling Markov chains, Zh. vjchisl. Mat. mat. Fiz., 7,3, 632-643, 1967.

5.

KEMENY, J. G., and SNELL, J. L., Finite Mwkov chains, Van Nostrand, 1959.

Constructing discrete Markov chains by the Monte Carlo method 6.

SOBOL’, 1. M., Multi-dimensionalquadrature formulae and Haar functions (Mnogomemye kvadraturnye formuly i funktsii Khaara), Nauka, Moscow, 1969.

I.

SOBOL’, I. M., On the dimensionality of Monte Carlo algorithms, in: Monte Carlo Methods and their Applications (Metody Monte-Karl0 i ikh urimeneniya), Novosibirsk, VTs SO AN SSSR, 138-l 39, 197 1.

45