13—
Computer Studies of Some History-Dependent Random Processes:
With W. A. Beyer and R. G. Schrandt (LA-4246, October 28,1969)

This report is about a novel method for the study of several non-independent probability schemata with rather curious results on patterns of growth, iteration processes, and dependent random walks. (Author's note).

Abstract

Various history-dependent random processes are investigated by computer and a few cases are investigated theoretically. These processes include historydependent random walks, a combination of a birth process and a self-avoiding random walk, historydependent randomly-generated increasing integer sequences, and randomly-generated integer sequences which might have prime-like densities. A possible random ergodic theorem for history-dependent processes is discussed.

I—
Introduction

In this paper we consider some examples of random processes in which the probabilities of the outcome of the n^th step depend upon the entire past history of the process. This, of course, means that they

are non-Markovian processes. In contrast to Markovian processes, very little is known theoretically about these history-dependent processes. They are much more difficult to analyze. However, the real world abounds with examples of the latter.

Most of the examples are results of computer studies. In some instances, theoretical results are known and are reported here.

The computations were performed on a CDC-6600 machine. The random-number generator needed had the form x -5 ¹⁵ xi_1 (mod2³⁵ with xo = 515). The scales on the figures are linear.

II—
Random Walk Examples

A—
Self-Avoiding Random Walk

A process which cannot, in any reasonable way, be made into a Markov process is a self-avoiding walk, i.e., a walk on a lattice starting at some fixed point that is not permitted to visit any point more than once. A survey of this topic will be published by Domb.¹ Another version of a self-avoiding walk is a scattering process in which scattering is not permitted at a point that previously had been a scattering center. The physical idea is that the particle, which had been the scattering center, is moved by the scattering process.

B—
History-Dependent Walk On the Line (Pólya)

An old example of a history-dependent process is the P6lya urn scheme described by Feller.² This has been used as a model of phenomena, such as contagious diseases, where an occurrence of a disease increases the probability of further occurrences.

As a special case of the P6lya scheme on two symbols, consider the following example. Let ao = 0, al = 1 and an+l = ai (i = 0, 1, ..., n) where ai means that one of the ai (i = 0, 1, ..., n) is selected at random with uniform probability. Let Yn- i=o ai Ti=n

Then Prob. [Yn = k/n] = 1/n ^- 1 for n > 1 and k = 1, 2, ..., n- 1. As discussed by Feller in Ref. 3, p. 237, Y = limnooYn exists with probability 1, since [Yn] is a martingale and the distribution of Y is uniform on [0,1]; i.e., Prob. [0 <a < Y < b< 1] = b-a.

The process Yn can be interpreted as a random walk on the horizontal line of integers with 0 interpreted as a step to the left and 1 as a step to the right.

C—
History-Dependent Walk On a Plane Lattice

A two-dimensional version of the Polya scheme on two symbols is a walk on the lattice of plane points with integer coordinates starting at the origin and is executed according to the following rule. The decision to execute either a horizontal step or a vertical step is made independently each time with Prob. [horizontal step (H)] = Prob.[vertical step (V)] = 1/2. After the decision H or V is made, a second decision is made. In the case of H, a decision is made to take the step right or left in accordance with the P6lya two-symbol game discussed previously. In a similar way, a decision is made to take a step up or down in the case of V. Figure 1 shows the terminal points of 10,000 walks of 64 steps each for walks made by these rules.

Fig. 1. End points of 10,000 random walks of 64 steps on the plane quadratic lattice starting at the origin in which the decision to execute a horizontal or vertical step is made independently with equal probability. The steps horizontally and vertically are made in accordance with the Polya two-symbol game.

For comparison, the distribution of 10,000 walks of 64 steps in the case of classical P6lya walks is shown in Fig. 2. Classical P6lya walk is the walk on the points with integer coordinates starting at the origin and selecting one of the nearest neighbors with equal probability. The formula for the probability that the nth step takes the particle to (x, y) is I r7 r₄ ,2₂ n j (cos a + cos 3)ⁿ cos xa cos y/dadd* .

It seemed easier to generate the distribution of end points by a Monte Carlo procedure than to use this formula.

Fig. 2. End points of 1000 random walks of 64 steps on the plane quadratic lattice starting at the origin in the classical Polya case.

A second two-dimensional version of the Polya scheme on two symbols is as follows. The probability of a horizontal or vertical step is itself a Polya scheme on two symbols. The remaining decisions are made as before and the resulting distribution is shown in Fig. 3. It is peaked about the coordinate axes in an approximately hyperbolic manner.

D—
A History-Dependent Explosion

A plane configuration on the plane lattice of points with integer coordinates is generated. Let the origin be the first generation. Assuming that certain of the lattice points have been occupied by generations

* See ref. 2 p. 371.

Fig.3. End points of 1000 random walks of 64 steps on the plane quadratic lattice starting at the origin in which the sequence of decisions to execute a horizontal or vertical step forms a P61ya two-symbol game. The steps horizontally and vertically are made as in Fig. 1.

up to and including the nth generation, let the (n + 1)th generation be determined as follows. For each point of the nth generation, two random numbers are selected to determine two neighbor points from four possible neighbor points for possible positions to be occupied in the (nth + i)th generation. These two neighbor points are to be occupied provided they have not been previously occupied. As an example (see figure below), suppose the point (0) is the initial point. Assume that positions (1) and (2) are chosen by the random numbers as the new positions. Then, in the first generation, the walk is made from point (0) to points (1) and (2), because they are unoccupied. Now there are two terminal points, (1) and (2), from which to walk in the second generation. If the directions (1) -— (0) or (2) -4 (0) are now chosen as

one of the new directions, this walk is not executed, because point (0) is occupied. If both directions (1) -- (3) and (2) -* (3) are chosen, the first walk is executed, but the second is not. The walk terminates if all possible directions chosen for all terminal points lead to points that are occupied.

Two computer examples are given for this type of walk. In the first case, a maximum of two random numbers were chosen for the possible new directions of the walk for each terminal point. The same direction could occur twice if the two random numbers were in the same interval kk+l-<rl, r2 <-,k = 0, 1, 2, or 3 4 - 4

In this case, only one direction for the walk was allowed for this terminal point. In this run, the walk terminated after 108 generations, with 656 total points. This is plotted in Fig. 4. Fig. 4. A history-dependent explosion. The rules for generation of this explosion are explained in the text.

In the second example, this situation was not allowed because the second random number was rejected and a new number was generated until two distinct new directions were obtained for each terminal point. This walk apparently does not stop. In Fig. 5 the walk is plotted through 65 generations, with 3126 points occupied. There are 100 active terminal points in the sixty-fifth generation.

The explosion is in some respects similar to a self-avoiding random walker, except the walker multiplies. Put in other terms, the explosion is in some respects similar to a branched polymer.⁴

Fig. 5. A second history-dependent explosion with different rules for generation as explained in the text.

III—
Integer Sequences Generated by History-Dependent Random Processes

In this section, sequences of integers obtained from the following three random processes are discussed. (a) an+I=an+ ai, i =,..., n (b) an+i=ai + ai, i,j = 1,..., n (c) an+1 = ai - aj, i,j 1,..., n, (i > j)

where a' has the meaning given above. In (a), al is given. In (b) and (c), al, a2 are given.

For (a), we will only mention that Kac and Ulam have shown that the expected value of an, E(an), is asymptotic for large n to eV. In case (b), it is easy to show by induction that 1 E(an) = - (al + a ₂)n for n > 2

One thousand sample sequences of case (b) were obtained by a Monte Carlo sampling with al = -1, a₂ = 1. Let a_n be the nth member of the kth sequence with 1 <n< 100 and 1 <k < 1000. The averages 1000 bn 1000 E aⁿk=l

are plotted in Fig. 6 as a function of n. In this case, E(an) = 0, n = 3, 4,.... It is seen that the averages increase with n. Fig. 6. The averages over 1000 sequences generated by an+l = ai +aj, ij = 1,..., n with al = -1 and a2 as a function of n.

In Fig. 7, a second example is given with al = 0 and a₂ = I and 1 <k< 5000. Then E(an) = 1/3n. However, we plot the deviation 1 5000 \~ I ki/ -n 5000 S a n3 1 -c=l as a function of n. The quantity apparently increases linearly with n.

In the case (c), it can be shown that the expected value of an, En = E(an), satisfies the recurrence relation 2n n-² - n+2En+=En+ En (n = 1, 2, ...) n+ln(n+ 1)

The asymptotic value of En is discussed in the Appendix and it is shown that En = 0(1/n). Figures 8 and 9 show a graph of En for al, a₂ = 0 and 1 and al, a2, = 1, 0, respectively.

Fig. 7. Graph of the function g(n) = 1/5000 E-5 I a- E(an) I for the case I/irvvv .k-1 an --=1 where an+i = ai + aj with al = 0, a2 = 1. {an}k=1,...,5000 is a random sampling of 5000 sequences.

Fig. 8. Graph of the function En defined by En+₂ = [2n/(n + )]En+l - [(n² n + 2)/n(n + 1)] En which is the expected value of an where an+l = ai - aj,, i,j = 1, 2, ..., n. Here al = 0 and a2 == 1.

Fig. 9. Same as Fig. 8, except that al = 1, a2 = 0.

Figure 10 shows a graph of the quantity 1 1000 bn = 1 E an n 1, ..., 100 1000 __ k=l where ak is the nth member of the kth sample sequence for case (c) al = -1, a2 = 1. The average increases with n, but here n remains small.

Figure 11 is a graph of the function 500 b- >1 f \ak-E(a) I bn = 500 E Jan- E(an) I k=l for n = 1, ..., 1000 with a = 0, a₂ =1. [ak] are 500 sample sequences of case (c). This function appears to be parabolic.

Fig. 10. The averages over 1000 sequences generated by a+ =ai - aj,i,j=1, . . ., n with al =-1, a2 = 1 as a function of n.

Fig. 11. The function bn = 1/500 E^0° = la -E(an)l for n = 1, 2,..., 1000 with an+l = ai - aj and al = 0, a2 = 1. {a }k=i ...,500 is a random sampling of 500 sequences.

IV—
Number Theoretical Games

In this section, we discuss sequences of positive integers generated by a history-dependent random process which might have densities like the primes.

Let di, d₂ , ..., dk be the sequence of differences of the first k primes. Let dk+l = di (i = 1, ..., k) and dk+₂= d + 2 ( = 1,..., k;,e i). In general, let dk+2j-1 = di (i = 1, . .., k + 2j-2) and dk+2j = ' +2 (E= 1,..., k + 2j - 2; £ $ i) for j = 1, 2, ... Let so = 1, si = si-i + di (i = 1, ..., n).

The following results are obtained.

1. For k = 8 and dl,..., d₈ , the first eight real prime differences, i.e., 1, 2, 2, 4, 2, 4, 2, 4, and n = 1000, 200 "games" are played. The average of s1ooo is 7986. The one-thousandth real prime is 7919. The average number of real primes in sl, ..., s1ooo is about 26%.

2. k = 49, d ₁, ..., d ₄₉ are the first 49 real prime differences. Two hundred "games" are played and the average s1ooo is 7643. The average number of real primes in sl, ..., slooo is about 30%. In 20 "games" it is found that the average slo,ooo = 101,356, whereas the ₁₀ ,₀₀₀ th prime is 104,723. The average number of primes in sl, ..., 10oo000 is 20%.

A second history-dependent random process for generating sequences with perhaps prime-like density is similar to the above. Again, let di,d ₂, ..., dk be the sequence of differences of the first k primes with k even. Let dk+i = (i = 1, 2, ..., k/2) and dk+l = / +2 (i = k/2 + 1, ..., k), etc., as before. For n = 1000 with 200 "games," the average value of slooo is 7892 and an average of 26% of the primes up to 7892 were obtained.

P. Stein⁵ has used such sequences to generate minimal binary additive bases for the even integers.

V—
A Problem

Starting with Ulam and von Neumann⁶ and Pitt⁷ there have been various versions of the random ergodic theorem stated and proved. One version is as follows.⁸ Let [X, -, li] be a measure space and [J, A, p] be a probability measure space with J a set of measure-preserving transformations defined on X. Let J* = 11i Ji where i = J for

all i and p* = pli° Pi where pi = p for all i. Then if F(x) e L1(X, t), it follows that {2 I n p* lim -1 f(Tk ... Tx) = f*(x) foralmostall xinX =1 k=1l where (T₁ , T ₂, ...) is a point of J* and f*(x) e L1(X, /).

In the above theorem, the transformations Ti are chosen independently. One can ask, "Suppose the Ti are not chosen independently, but are chosen to form a Markov chain, or are chosen in accordance with the rules governing a P6lya process (see above), does some random ergodic theorem hold?" One could specialize to the case of choosing from two transformations, each of which is a rotation of the circle.

Appendix

Theorem.E(k) = O(1/k) where k(k + 1)Ek+₂ - 2k² Ek+l + (k² - k + 2)Ek =(k = 1, 2,...) . (1)

Proof. Define the vector (k)-k)¹ (k)1 - e(k)1) e2(k) [e(k )J Then Eq. (1) is equivalent to the vector equation 0 1 £(k + 1) = k) (2) kc-l 2 2k k+l -k(k+l) k+l

We first show that any vector solution to Eq. (2) is bounded for all k> 1 in the norm I1£1j = 111 + \121. Write 0 1B(k)_k-l 2kk+l k+l and o o A(k) _ 2 0 k(k+l)

Consider £(k + 1) = B(k)£(k). (3) A fundamental matrix set of solutions for Eq. (3) (Ref. 9) is given by Yk= (4) 1 k for k = 2, 3, .... From work shown in Ref. 9, any solution £(k) to Eq. (2) for k> 2, can be expressed by k-1 £(k) = E Yk-s+1 Y2 ^- A(s + 1) £(s) + YkY2-£(2) . (5) s=2 Since 00 EA(k) < oo, k=l

it follows from Lemma 3.2 of Ref. 9, p. 21, and Eq. (5) that l£(k)l is bounded for k> 2 and hence for k> 1. Now a computation gives 1 11- k -s 4£1(s) Yk-s+l Y₂A(s + 1) (s) = - - k1(s +1)(22) - k sl (s 1)(2 4 2) ' and • - (2) + 22(2) + £1(2) - £2(2) YkY2-'£(2) =k = 1(2) + 252(2) + k21(2) - e2(2)

Thus k-i [1 - 41(s)5(k) =-Zk-s41 (s) s=21( + 1)( + 2) k-s+l + (6) [-£1(2) + 2£2(2) + k21£1(2) - £ E2(s) -£1(2) + 2(2) 1(2 ) + (2 - 1(2)

for k > 2. The summation can be written k-¹ _ k-s 44£1(s) =2 L+11(s + 1)(s + 2) 1 51k-s+l k-I £ l()k-I [ 1_](s) -4 + 4 k⁺⁴ (7) -4 + ( 1)(s + 2) '1 [ 1(s + 1)(s +2) ,9=2 8=2 k-sZl

Since 5¹ (s) is bounded, it is seen that the first term on the right of Eq. (7) can be written k-—1 oc 0i k-4¹ (s) -4¹ (s) 4 (s+1)(s+2) -4 (s+1)(s+2) s-2s=2 +4 ( 1£)(s (8) s=-k (k) where C N= -4 s -i 1)($ (+ 2) is a constant. E (s + 1)(s + 2)

Now consider the second term on the right of Eq. (6). For 2 <s<k/2 1 2 k-s ki11 2 = O(k₎k-s+1 - (k/2 +1) k+2 k and for k/2 < s <k- 1 1 1 4 (k- s)(s + l)(s + 2) (k/2 + 1)(k/2 + 2) (k + 2)(k + 4)' 1 2 (k-s+1)(s+1)(k+2)^- (k+2)(k+4)

Thus, with [] denoting integer part, we have k-1 l(s)[k/2](s)(k-s)(s+ )(s+2) (k - s)(s + l)(s + 2) (9) k-1El¹ (s) +(k 1)(s = 2) -0(1/k) -i kO(1/k² ) = O(1/k) Z=[kE]+ (k-s)(s + 1)(s + 2) s=[k/2]+1

and similarly for k-i 1(s) (k - s + 1)(s+l)(s +2) s=2 Hence, from Eqs. (6), (8), and (9), we have £(k) = A + 0(1/k), where A = [] is a constant vector. Therefore, c(k) = A + 0(1/k) Substituting this into Eq. (1), we obtain [k(k + 1) - 2k²+k² -k +2] A+[k(k + ) - 2k²+k² - k + 2]0(1/k) = 2(A + 0(1/k)) = 0 Thus, A = 0, which completes the proof.

An Approximating Differential Equation

If Eq. (1) is written in the form E(k + 2Ak) -2E(k + Ak) + E(k) 2 E(k + Ak)-E(k)(Ak)²k +1 Ak + k(k E(k) = 0, Ik(k + 1)

with Ak = 1, then the differential equation 2 2 e" + - k(k+ e=O, (10) could be regarded as an approximation to Eq. (1).

In Eq. (10), k is regarded as a continuous variable. It can be shown that the general solution to Eq. (10) is a linear combination of the functions F( ^l+i V^x/1- i.2, k and ( 2 2 2,1 + where F denotes the Gauss hypergeometric function. Hence, for large k, the solution to Eq. (10) is a linear combination of the functions el and e2 where ei= - cos[ i log k + 0 e 7L2 kJ³ /² and e2=- sin[ -log k ]+ ( k/

Regarding Eq. (10) as an approximation to Eq. (1), it can be shown from the theory of difference approximations to differential equations that if ko is fixed and E(ko) = e(ko) as well as E(ko + 1) = e(ko + 1), then for k > ko we have Ie(k)-E(k)\< ,e(~) [e,_] (11) e(k) I 41 e(k)I e where ko< 5 < k. Thus, while it is true that e" (W) -x oo, the exponential factor in Eq. (11) prevents Eq. (10) from being a very good approximation to Eq. (1) for arbitrarily large k. This difference equation can also be discussed, perhaps more successfully, by using generating functions.

References

1. C. Domb, "Self-avoiding Walks on Lattices," to be published in Advan. Chem. Phys.

2. W. Feller, An Introduction to Probability Theory and Its Applications, (John Wiley and Sons, New York 1968), Vol. 1, 3rd ed.

3. W. Feller, An Introduction to Probability Theory and Its Applications, John Wiley and Sons, New York 1966), Vol. 2.

4. L. V. Gallagher and S. Windmer, "Monte Carlo Study of Flexible Branched Macromolecules," J. Chem. Phys. 44, 1139 (1966).

5. P. Stein, private communication. 6. S. M. Ulam and J. von Neumann, "Random Ergodic Theorem," Bull. Am. Math. Soc. 51, 660 (1945).

7. H. R. Pitt, "Some Generalizations of the Ergodic Theorem," Proc. Cambridge Phil. Soc. 38, 325 (1942).

8. P. Revesz, "A Random Ergodic Theorem and its Application in the Theory of Markov Chains," in Ergodic Theory, (Academic Press, New York, 1963) Fred. B. Wright, Ed.

9. K. S. Miller, Linear Difference Equations, (W. A. Benjamin, Inc., New York, 1968).

Everett and Ulam in Madison, Wisconsin, in 1941.