A program in BASIC for facies-by-facies Markov chain analysis

Computers & Geosctemes Vol. 15, No. I. pp. 143-155. 1989 Pnnted in Great Britain. All rights reserved

0098-3004 89 $300 + 0,00 Copyright ~ 1989 Pergamon Press pie

A PROGRAM IN BASIC FOR FACIES-BY-FACIES MARKOV CHAIN ANALYSIS NElL A. WELLS Geology Department, Kent State University, Kent, OH 44242. U.S.A. (Received 22 December 1987; accepted 15 June 1988)

Abstract--This paper presents a BASIC program for Markov analysis following a method proposed in 1982 by D. W. Powers and R. G. Easterling. which is appropriate for facies-by-facies analyses, in which a facies may not succeed itself. It estimates levels of significance by Powers" and Easterling's chi-square method and also by C. W. Harper's binomial probability method. A sample analysis of a series of Eocene fluvial deposits from Pakistan shows the benefits of running several versions of each data set. and of considering whether "apparent randomness" of transitions is expected and irreducible. Randomness is taken usually as the null hypothesis in Markovian analyses, but many sedimentological processes should cause distinctly nonrandom facies successions. Therefore. more effort should be put into considering the meaning and causes of "'randomness." Key Words: Stratigraphy, Facies successions, Markov chain analysis, Randomness, BASIC.

INTRODUCTION The purpose of this paper is to help make the use of Markov chain analysis more convenient and thus more widespread in stratigraphic interpretation. Markov analysis can indicate whether transitions between particular facies might be occurring randomly or whether they sccm to be favored or disfavored. Specifically, it is used to test whether available data are consistent with the null hypothesis that tr.'msitions between particul:tr tacies occurred at frequencies indistinguishable from random succession. Strongly favored transitions then can be linked to determine an ideal series of successions. This paper presents a computer program, in BASIC, for Markov analysis, and offers some comments on the use of the method. Important points are the need to use the PowersEasterling or Carr methods in analyses where transitions from a facies to itself are prohibited, and the need to consider the meaning of random transitions, which generally are ignored. A computer program for Markov analysis following the Powers and Easterling (1982) method seems to be needed for two reasons. First, after too many arduous and error-prone hand-calculations finally goaded me to write a Markov analysis program. I discovered that many colleagues (who, I had assumed, had long since computerized their analyses) also were doing them by hand. Second, experiences at recent geological conferences indicate that too many researchers are analyzing facies-by-facies successions by the type of procedure used by Gingerich 11969), Read (1969), and Miall (1973, 1984), even though Carr (I 982) and Powers and Easterling (1982) showed these methods to be inappropriate for such analyses. In brief, the problem is that early methods do not make corrections for the empty cells ("structural zeros") that are created diagonally across the transition macAozo Is;,-#

143

trix because a facies-by-facies analysis by definition excludes instances of a I:acics succeeding itself. NII.~TIIOI) OF ANAI.YSI.% A Markov analysis begins with a matrix of ohserved transitions ("raw-data matrix" in sample output, Fig. I), which is a compilation of the nnmbcr of times a facies on the left passes up into a facies across the top. The data can be edited into the program alter line 3010, which is most convenient, or they can be fed in via the keyboard by answering N (no) in line 70 and Y (yes) on line 310. Expected values then are calculated for each type of transition ("expected values", Fig. I; lines 1250 and 1570), using the assumption that all tr:msitions (except for self-succession) occur randomly and therefore depend only on the relative abundances of the two facies involved. The matrix of expected values then is subtracted from the matrix of observed values, yielding the difference matrix (lines 1260 and 1580). Negative values indicate disfavored transitions (ones occurring less than randomly), and positive values identify favored transitions. After this stage, a variety of subsequent manipulations now is possible (see suggestions by Powers and Easterling, 1982; Harper, 1984a, 1984b). The first refinement provided by this program is calculation of the chi-square matrix in Figure I (line 1270), using the equation [ o b s e r v e d - expected]"/ expected. The sum of the chi-square scores and the degrees of freedom (calculated in lines 1280-1290 and 1680-1720) should be compared to a table of chisquare values: if the sum exceeds the given value, then the overall matrix is nonrandom. Another refinement is a test of symmetry (the asymmetry matrix: Powers and Easterling, 1982, p. 920-921). This tests whether transitions generally are reversible, which is to say whether transitions from

144

N.A.

WELLS

MARKOV ANALYSIS LOWER KULDANA FORMATION, EOCENE, PAKISTAN, all f a c i e s Number of facies

I0

FACIES: l 2 3 4 5 6 7 8 9 I0

Pebbles Calc. Granules Festoons Tabular Dunes Plane Beds Ripple Troughs Other Ripples Non-red Clay Red Clay Channel Floor

RAN-DATA MATRIX # l 2 3 4 5 6 7 8 9 10

# l 0 0 I 2 0 0 O 0 O 14

# 2 4 0 8 I 1 I 2 2 12 51

# 3 8 30 0 I 6 2 O O 4 27

# 4 2 5 7 0 0 O 3 O 0 4

# 5 0 5 22 4 0 5 0 I 5 5

# 6 0 4 17 3 12 0 O 2 6 8

# 7 1 3 4 2 6 4 0 0 2 2

# 8 0 8 3 0 5 0 2 O 12 0

# 9 2 28 20 5 19 29 14 20 O 0

#10 O O 3 I 3 2 O I 104 0

17 83 85 19 52 43 21 26 145 11!

17

82

78

21

47

52

24

30

137

114

602

ROW & COLUMN FACTORS

Row T o t a l s 1 I 2 2 3 3 4 4 5

r c r c r e r c r

17 1.09 0.26 1.64 0.26 1.64 0.26 1.64 0.26 !.64

83 9.22 1.42 9.04 1.40 9.01 1.40 9.00 1.40 9.00

85 9.44 1.36 9.20 1.34 9.16 1.34 9.16 1.34 9.16

19 2.11 0.32 1.85 0.32 !.84 0.32 1.84 0.32 1.84

52 5.78 0.77 5.29 0.76 5.27 0.76 5.27 0.76 5.27

43 4.78 0.84 4.40 0.82 4.39 0.82 4.39 0.82 4.39

21 2.33 0.37 2.05 0.37 2.05 0.37 2.05 0.37 2.05

26 2.89 0.47 2.57 0.46 2.56 0.46 2.56 0.46 2.56

145 16.11 2.70 18.35 2.79 18.52 2.80 18.54 2.80 18.54

111 12.33 2.09 13.04 2.10 13.02 2.09 13.02 2.09 13.02

66.89 I0.60 67.44 10.62 67.47 10.62 67.47 10.62 67.47

EXPECTED VALUES # 1 2 3 4 5 6 7 8 9 10

1 2.33 2.36 0.48 1.36 1.13 0.53 0.66 4.79 3.36

2 2.30 12.84 2.59 7.39 6.16 2.87 3.59 26.00 18.26

3 2.19 12.04

2.47 7.05 5.87 2.74 3.42 24.80 17.41

4 0.52 2.88 2.93 1.69 1.40 0.66 0.82 5.93 4.17

5 1.24 6.80 6.92 1.39

3.32 1.55 1.93 14.01 9.84

6 1.35 7.42 7.55 1.52 4.35 1.69 2.11 15.28 10.73

7 0.60 3.30 3.36 0.68 1.93 1.61 0.94 6.80 4.78

8 0.76 4.16 4.23 0.85 2.44 2.03 0.95

$ -1.24 -1.80 15.08 2.61 -

6 -1.35 -3.42 9.45 1.48 7.65

7 0.40 -0.30 0.64 1.32 4.07

8 -0.76 3.84 -1.23 -0.85 2.56

8.57 6.02

9 4.59 25.21 25.64 5.16 14.76 12.29 5.73 7.17

36.45

tO 3.43 18.85 19.17 3.86 11.04 9.19 4.29 5.36 38.81 -

DIFFERENCE MATRIX 1 -2.33 -1.36 1.52 -1.36

2 1.70

-4.84 -1.59 -6.39

3 5.81 17.96 -1.47 -1.05

4 1.48 2,12 4.07

-1.69

Figure I

9 10 -2.59 -3.43 2.79 -18.85 -5.64 -16.17 -0.16 -2.86 4.24 -8.04

Continued opposite

145

M a r k o v c h a i n analysis 6 7

8 9 10

-1.13 -5.16 -3.87 -0.53 -0.87 -2.74 -0.66 -1.59 -3.42 -4.79 -14.00 -20.80 10.64 32.74 9.59

-1.40 2.34 -0.82 -5.93 -0.17

1.68 -1.55 -0.93 -9.0l -4.84

-1.69 -0.11 -9.28 -2.73

5 1.24 0.48 32.87 4.87

6 1.35 1.58 11.84 1.44 13.49

2.39 -0.94 -4.80 -2.78

-2.03 1.05

16.71 8.27 12.83 3.43 -6.02 -36.45

-7.19 -4.29 -4.36 65.19

SEQUENCE CHI-SQUARE MATRIX # 1 2 3 4 5 6 7 8 9 10

1 2.33 0.79 4.87 1.36 1.13 0.53 0.66 4.79 33.67

2 1.25 1.83 0.97 5.53 4.32 0.26 0.70 7.54 58.72

3 15.36 26.77 0.87 O.16 2.55 2.74 3.42 17.44 5.28

4 4.14 1.56 5.65 1.69 1.40 8.39 0.82 5.93 O.O1

0.85 1.55 O.45 5.79 2.38

1.69 O.O1 5.64 0.69

7 0.26 0.03 0.12 2.59 8.55 3.55 0.94 3.39 1.61

-

8 0.76 3.54 0.36 0.85 2.70 2.03 1.17 1.37 6.02

9 1.46 0.31 1.24 O.01 1.22 22.72 11.92 22.99

8

9

10 3.43 18.85 13.64 2.12 5.85 5.62 4.29 3.54 109.47

36.45

ASYMMETRY CHI-SQUARE MATRIX # 1 2 3 4 5 6 7

8 9 In

1

2

3

4

5

6

7

4.00 5.44 0.OO O.OO 0.00 1.O0 O.00 2.00 14.O0

12.74 2.67 2.67 1.80 0.20 3.60 6.40 51.OO

4.50 9.14 11.84 4.00 3.00 10.67 19.20

4.00 3.00 0.20 O.OO 5.OO 1.80

2.88 6.00 2.67 8.17 O.50

4.00 2.00 15.11 3.60

2.00 9.OO 2.00

sequence chi-square

71

sum = 6 2 8 . 6 3 3

asymmetry chi-square

= 348.796

41

degrees

2.00 1.O0 1 0 4 . 0 0

degrees of freedom

of f r e e d o m

PRINCIPAL SUCCESSION FOR EACII FACIES binomial probabilities listed at right give probability (confidence level o I - level of significance) o f t h e o b s e r v e d f r e q u e n c y o c c u r r i n g by c h a n c e

Pebbles t o F e s t o o n s (conf. > .99; BP- 0 . 0 0 0 6 ) Calc. Granules to Festoons (conf. > .99; BP- O.OOO0) Festoons to Plane Beds (conf. > .99; BP= 0.OOOO) T a b u l a r Dunes t o P l a n e Beds (conf. .95-.99; EP0.0461) Plane Beds to Ripple Troughs (conf. > .99; BP- O.0010) Ripple Troughs to Red Clay (conf. > .99; BP- O.0000) Other Ripples to Red Clay (conf. > .99; BP- 0.0002) Non-red Clay to Red Clay (conf. > .99; BP- 0.0000) Red Clay to Channel Floor (conf. > .99; BP- O.OOO0) Channel Floor to Calc. Granules (conf. > .99; BP- 0.0OO0) Other Significant

Favored Transitions

( c o n f . > c a . 0 . 9 5 & <20% p r o b a b i l i t y (NB: b i n o m i a l p r o b a b i l i t i e s between

Pebbles t o T a b u l a r BP= 0 . 0 9 5 5

Significantly

Disfavored

Transitions

of selection by c h a n c e ) .1 and .2 d o n ' t mean much)

Dunes

Pebbles to C h a n n e l Floor BP" 0.O216 Calc.

to Non-red Clay BP- 0 . 0 5 6 5

Granules

Calc. Granules t o C h a n n e l F l o o r BP- O.OOOO Festoons t o Tabular Dunes BP- 0 . 0 2 7 6 F i g u r e I.

Continued overleaf

I0

146

N.A. 't~'ELLS Festoons t o R i p p l e Troughs BP= 0.0012 Festoons to Channel F l o o r BP= O.O001 T a b u l a r Dunes t o Pebbles 8P= O . 0 8 l l Plane

Beds t o C a l c . BP= 0 . 0 0 3 3

Granules

Plane Beds to Other Ripples BP = 0.0124 Plane Beds to Channel Floor BP= 0.0022 Ripple Troughs to Calc. Granules BP= O.0107 Ripple Troughs to Other Ripples BP = 0.0765 Ripple Troughs to Channel BP= 0.0026 Other Ripples to Festoons 8P= O.OS31

Floor

Other Ripples to Tabular Dunes BP= 0.0266 O t h e r R i p p l e s t o Channel F l o o r BP= 0 . 0 0 8 3 Non-red Clay to Festoons BP= 0 . 0 2 5 5 Non-red Clay to Channel F l o o r BP- 0 . 0 1 9 2 Red Clay to Pebbles BP= 0.0077 Red C l a y t o C a l c . G r a n u l e s BP- 0 . 0 0 0 8 Red Clay to Festoons BP= O.OOO0 Red C l a y t o T a b u l a r Dunes BP- 0 . 0 0 2 3 Red Clay to P l a n e Beds BP° 0.OO40 Red Clay to Ripple Troughs BP- 0.0045 Red C l a y to O t h e r R i p p l e s BP- O.O315 Channel Floor to Pebbles BP= 0 . 0 0 0 0 Channel Floor to F e s t o o n s BP" 0.0118 Channel Floor to Non-red Clay BP- 0 . 0 0 2 1 C h a n n e l F l o o r t o Red C l a y BP- 0.0000 BINOMIAL PROBABILITY

FOR CELL

2

BINOMIAL PROBABILITY

FOR CELL

3 , 9

.110711

BINOMIAL PROBABILITY

FOR CELL

5 , 9

.125967

,

9

=

.288408

Figure I. Sample output from program, using data from Early to Middle Eocene redbeds in Pakistan facies A to fifties B occur as frequently as transitions from B to A. For each such pair of facies, the program calculates the dilt'crence between the number of transitions from A to B and the number from B to A, v, hich it squares and then standardizes by dividing by the number of transitions involved (A --* B plus B ---, A).Thisisdone in line 1330. Lines 1310and 1320 prevent division by zero when no transitions were recorded for either of a pair, and lines 1350-1360 count the total number of empty and identical pairs for purposes of calculation of degrees of freedom (lines 17101720). The matrix would be symmetrical about its diagonal, so only one-half is produced (lines 1300 and 1600). The values are summed (line 1340) and are

treated as chi-square scores. The higher the score, the more unidirectional (nonreversible) the transition. The program makes several initial selections of "'significant" results. First, the difference matrix is searched once for the principal preferred succession from each facies, and again for any other "significantly" preferred and disfavored transitions. Significance is not defined uniquely and simply in Markov analyms. Harper (1984b) notes that residuals (the square root of the individual chi-square components) have values that are distributed approximately normally and therefore give an approximate indication of levels of significance: by this reasoning, chi-square scores of 2.71 are significant at the 0.05 level, whereas scores of

Markov chain analysis 5.38 have a 0.01 significance. Approximate levels of significance are listed for principal successions. The program then calculates the binomial probability of getting at least the observed n u m b e r of transitions by chance (lines 2..190--2730: see Harper, 1984a). Harper suggests accepting as significant only transitions that have odds < 0.10. Binomial probabilities are listed for all occurrences. Here, transitions other than principal transitions are listed as significant if their chi-square score exceeds 2.71 a n d if their binomial probability exceeds a I in 5 chance of random occurrence. Both limits are generous deliberately so that the list may be over-inclusive (but binomial probabilities are listed and scores can be located in the chi-square matrix). Also note that because the searching routine first uses the difference matrix, "other transitions" for some facies in fact may have more significant chi-square scores and binomial probabilities than the matching "principal transition". Before ending, the program gives the user the option of obtaining the binomial probability for any desired transitions. The heart of the binomial probability calculation is the calculation n C r • p ' . (I - p)q"-'~, in which n = the total n u m b e r of transitions from a facies, r = the observed n u m b e r of those transitions that go to another particular facies, p -- the expected frequency of that transition, and n C • = the number of combinations o f n things taken r at a time (see tlarper. 1984a). For preferred transitions, we are interested in knowing the probability of getting the observed number of transitions or nacre, so the equation has to be calculated and summed from r = • to r = n. (Disfavored transitions, however, require the probability of the observed n u m b e r or h'ss, so calculations must be evaluated from r = r to r .= 0 instead.) For example, if the calculated probability of a single greaterthan-random transition from one facies to another is 0.20 and if 7 of the 10 transitions from the first facies were indeed observed to go to the second facies, binomial probability calculations, summed for r = 7 to r = 10, tell us that the chance of getting at least 7 such transitions by chance is 0.000864. This explanation is needed to explain some programming shortcuts that may not be understandable on first examination. Although n C • could be calculated in each iteration (from r = r to r = n or to • = 0) by n!/[r! * (n - r)!], it is solved by lines 26202660. First, j,CT, or 10!!(7! * 3!), can be simplified to 10 * 9 * 8/3 • 2 * I, a calculation that can be produced clearly by a routine that runs n - r times. In essence, the routine is F O R J = ! TO n - r: N U M E R A T O R = N U M E R A T O R • (n - (J I)): D E N O M I N A T O R = D E N O M I N A T O R • J: N E X T J. The calculation is done in logarithms, as variables otherwise rapidly exceed maximum values. Because of the symmetry of combinations (e.g. ~oC, = roCk), this routine also can be used also for disfavored transitions, if n - r is substituted for r in "'FOR J = 1 TO n - r'" in line 2630 and ifthe entire routine is run from n - r, instead of from r, to n times

147

(line 2600). In the routine, n is represented by Y2 and F(R), which is the r o t total for facies R, and, if the transition is favored, r is represented by A(R. C) (the observed number of transitions from R to C). If it is disfavored. C(R, C) will be < 0 in line 2590 and the equivalent of n - r will be substituted. Because Z in line 2600 can be based either on r or n - r. two different probability calculations are needed (lines 2680-2710). Calculation of probabilities for frequencies such as 50 in 150 takes a long time, but because the routine ensures that each n e t iteration contributes less to the total binomial probability ( = BP in line 2720), it is possible to truncate the iterations once iterations change BP by < I,'1000 of its value (line 2730). The difference between the Gingerich-Read method and the P o t e r s and Easterling method lies in the calculation of the expected values and their appropriateness depends on the type of data set under consideration. The Gingerich-Read method calculates the expected cell value in row i, column j (i.e. the expected n u m b e r of transitions from [attics i to facies j ) by dividing the total number of facies i beds by the total number of non-./' beds. In bed-by-bed or "foot-by-loot'" analyses, where data are collected from catch bed or at specilic intervals, this analysis is appropriate. (The interval must be selected cltrcfully as if it is too large, transitions will be missed and the matrix will be randomized, but if it is too small the matrix will be swamped by transitio,ls from each facies to itself.) l-lowcver, most analyses will bc lacicsby-facies, and successions of one facies by itself wilt not be recognized. As described earlier, this results in structural, or embedded, zeros along the diagonal of the matrix. Powers and Eastcrling (1982) and Carr (1982) show that the Gingerich-Read method, on occasion, can yield important errors when used with such at data set. Powers and Eastcrling propose instead an itcrative procedure for calculating row and column factors that can be multiplied to give expected cell values. The first iteration for each row factor a, begins by dividing each row total by the number of columns less 1 (lines 620-650). Next, each column l:actor bj is calculated as the column total for column j divided by the sum of previous row factors from the previous iteration less the factor for r o w j (lines 660760, especially 721)). Beyond the first iteration, each row factor a, equals the row total for row i divided by the sum of column factors from the previous iteration less the factor for column i (lines 810-870, especially 830). Iterations are continued until each factor differs by < 0.01 from its equivalent in the previous iteration (lines 700, 750-760, 850-860, and 930). An example of the differences resulting from the two methods can be seen by comparing the results of Harper (1984a, 1984b) and this program for analyses of the Battery Point Sandstone data of Cant and Walker (1976). Results can be obtained from this program by converting D A T A statements to REM statements and vice versa in lines 3020-3290. Harper

N. A. ~Id~-LLS

I.aS

(1984a, fig. I) calculated binomial probabilities for the data by using "'Gingerich-Read'" expected values and he accepted seven transitions as significant. This program, using "Po~ers and Easterling'" values in otherwise equivalent calculations, rejects two of these transitions (,4 --. B, and F -+ Ssl at the levels of significance selected b,', Harper. DISCUSSION Much has been v, ritten about Markov analysis by other authors, and apparently authoritative statements are being refined or refuted continually. Therefore it is with considerable trepidation that 1 wish to add several minor points. First. in a talk at the Third Internatiomtl Fluvial Sedimentology Conl~'rence in 1985. Miall expressed the opinion that Markov analysis could no longer be considercd a valuable method of investigation. However. improved methods of analysis have surpassed the method used by Miall (It)84. p. 147){e.g. Powers and Easterling, 19S2: llarper, It)S4a, I q84b). Furthermore, fi~r all of the problems, t, ncertaintics, and diffict, lty of Markov analysis, publication of observed. dilt'crcnce and chi-square matrices seems to be the clearest, shortest, most complete, and least biased way of both presenting data concerning tacies transitions and suggesting ideal successions. In the absence of such a presentation, x~,c are faced with detailed and untligcstcd sections, purportedly "'typical" st, bscctions, or, worst of all, unsupported assertions of "'ideal'" sequences. For cxarnple, in live 1986-1987 issues of tile Jottrnal of N<'dimentary Petroh(gy, seven ol'eight papers (more-or-lcss) that were concerned signiticantly with recognition of typical sequences forwent statistical analysis in favor of such unsatisfactory alterm, lives. In short, more Markov :malyscs should be done, not fewer. Second, Markov analysis usually is performed only on long sections, because mzmy transitions are needed to make the results meaningful, and because use of a single section precludes the possibility of double-counting of spccilic successions, which could create the false appearance of"preferred" transitions. tlowever, it seems reasonable that with appropriate care SOille short equivalent sections could be added together, when they represent dissimilar scqt, ences or separate lithosomes. For example, it would seem legitimate to use Markov analysis to assess typical bedform successions in a series of pro-l'limalay:m rivers by combining one sequence from each river. Another example might be combining one beach protile from each Caribbean beach or island. When interpreting the rest, Its. one would have to kccp in mind the possibility of apparent cyclicity that is event+ rchtted rather than process-related. For example, all the rivers might contain a type of sequence related to one particular basin-wide phenomenon such as the last rnonsoon season. In the situation of the beaches, different island sections might reflect regional changes

in scale, el. and same island sections might be influenced by local uplift or sediment availability. If such possibilities can be dismissed, then each section can be slewed as a separate succession, just as if the', all occurred in a vertical sequence. Again with the proviso of needing to consider event-related effects, it likewise would seem reasonable to add laterally equivalent sections itindividual subunits have changed so much that the? cannot be correlated from one section to the next. Even if the independence of each sequence is in doubt, such a compilation is ~orth~shile because 700 km:, few sections are closer than 3 kin, and the sandstones occur reasonably equally at all levels (Wells, 1984, figs. 4-10), so the sandstones are viewed as separate entities. Despite the typical and sedimentologically expectable fining-upward successions, many transitions were observed from most major facies to rcd mud. even though such transitions to mud turn out to occur "'randomly". In fact, 66 of 114 channels show

Markov chain analysis A.

Preferred transitions between Lower guldana facies, with chi-square scores ) 2.7 or binomial probabilities < 5%, from output #I. 34

5

~.PEBBLES¢- . . . .

4

8

-~TABULAR DUNES(

RED CLAY

OTHER RIPPLES

13 NE BEDS

S

" B.

149

," j . . , R E D CLAY :RIPPLETROUGHS~'~23~23

--

3..'LNON-RED CLAY

Simplified relationships, from output #I (chl-square scores above 12 & B.P. (1%).

RED CLAY

....~PEBBLES~ ~~' .~. r. ~ r cr ~ vn ~CHANNEL FLOOR ~ . .n. .~. c .

) PLANE BEDS

)RIPPLETROUGHS

--GRANULES~

C.

7)RED CLAY NON-RED CLAY

Preferred transitions, with chi-square scores, from output #2 (excluslon of channel floors). Minor facies show unstable results. 25 5 PEBBLES, TABULAR DUNES~ OTHER RIPPLES

12 RED CLAY

...s 3

, GRANULES

) FESTOONS

7

s2

--.

~PLANE BEDS

_ _12. . . . .

;'RIPPLETROUGHS

"" ~'". s

t, ~RED CLAY~i 0 NON-RED CLAY~ II

D. Preferred transitions between simplified (combined) facies, from output #3.

CLAY

PEBBLES & CHANNELJGRANULES"~MEGARIPPLES "~PLANE BEDS"~RIPPLES i FLOOR ~ " ~

~CLAY

Figure 2. Graphic summary of results ot" Markov analysis in Figure I. A --"signilicam'" transitions selected by computer~ B--strongest transitions, with chi-squares scores > 12~ C--rerun of data withoul category "'erosive chann,:l floor"; D-simplilied matrix wherein sever~,l facies are combined (l + 2.3 + 4, 6 + 7, ;,nd 8 + 9).

some sort of abridgement of the ideal succession by transition to clay. However, usual, essentially random, curtailment of fluvial sequences is not expected in normally developing fluvial deposits. The author has suggested that this high frequency of abridgment and the apparent under-development of the sandstones (almost none show evidence of meandering despite an approximate 18:i overbank clay to channel sand ratio) indicate that avulsion and abandonment were such important aspects of this system that about one-half the streams became abandoned before becoming mature and meandering. This would be an example of randomly occurring transitions that occur more frequently than expected because of addition of facies (red mud) by a stochastically interfering process (avulsion), Genuinely random transitions of course also could be occurring less frequently than expected: randomness could result from destruction or deletion of intervening facies. Wells (1984, chap. 5) describes a situation where the end members of a highly symmetrical and repetitive succession of small transgressions and

regressions are juxtaposed more than expected, which he ascribes to erosion of intermediate units under a marine transgression and pedogenic destruction of intermediate units under peak-regression redbeds. In short, we should consider carefully geological logic and ramifications when we use randomness as a null hypothesis. Fourth, most data sets can be organized legitimately in several different fashions, so running several different versions usually is instructive. (Using several different analytical procedures likewise helps to spotlight robust transitions that always show up and to identify additional marginally significant ones.) The sample data set in Figure I was rerun without the category "erosive channel floors", because not all authors would utilize this as a separate facies category and because it contributes substantially to the overall chi-square score. Rerunning showed the variability of results concerning minor categories (Fig. 2C). The data were rerun yet once more by combining several facies: "'pebbles" and "granules" were combined as "'coarse channel lag". tabular and tangential large-

150

N . A . WELLS

scale crossbeds were c o m b i n e d as "large-scale crossbeds", small-scale crossbeds were c o m b i n e d as "'ripples", a n d red a n d white clays were c o m b i n e d as "'clay"( Fig. 2D). C o m b i n i n g facies has the a d v a n t a g e s of decreasing the n u m b e r s o f m i n o r transitions a n d m a k i n g the matrix more robust, but it involves a potential loss o f information. Instead o f worrying a b o u t the trade-off, the data can be run easily b o t h ways. a n d the results c o m p a r e d . CONCLUSION This p r o g r a m is designed to reduce m u c h o f the work o f d o i n g a M a r k o v analysis a n d will save a lot of time ordinarily spent in looking at matrices in o r d e r to construct preferred successions, but the analyst is nevertheless urged to spend time p o n d e r i n g the raw data a n d difference matrices, r e r u n n i n g different versions o f the data, a n d considering the sources a n d m e a n i n g o f a p p a r e n t randomness.

,4,'knmvh'd.t,ments--! would like to thank Dr. Philip Gingerich for starting and funding my work in Pakistan (via grants from the Smithsonian Institution Foreign Currency Program and the National Science Foundation) and for persuading me to redo my Markov analyses following the method of Powers and Easterling. I also would like to thank l)r. tlolly Smith-Gingerich for explaining the method to me. and. especially, l)r. Charles Ilarpcr fi~r it helpful review that considcrahly improved both the paper and the program. Two other anonymous reviewers made helpful comments. I also would like to thank Britt Leath;tm for comments on adapting this progr;tm to BASICA. Please note that the individuals named have not seen the results of their contribulions, so any errors or misrepresentations herein are the fault of the anthor ;done.

REFERENCES Cant. D. J.. and Walker, R. G., 19"76. Development of a braided fluvial facies model for the Devonian Batter.,,' Point Sandstone. Quebec: Can. Jour. Earth SCi., v. 13, no. I, p. 102-119. Carr, T. R.. t982. Log-linear models. Markov chains and cyclic sedimentation: Jour. Sed. Pet., v. 52. no. 3, p. 905-912. Gingerich, P. D., 1969, Markov anal?sis of cyclic alluvial sediments: Jour. SCd. Pet.. v. 39. no. I. p. 330-332. Harper, C. W., 1984a. lmpro,,ed methc,ds of facies sequence analysis, in Walker, R.G.. ed., Facies models: Geoscience Canada Reprint Series no. l {2nd ed,), p. 11-13. Harper, C. W,, 1984b. Facies models re,.isited: an examination of quantitative methods: Geoscient:e Canada. v. I I, no. 4. p. 203-207. Mi/.ll, A. D., 1973. Markov chain anal?sis applied to an ancient alluvial plain succession: Sedimentology. v. 20, no. 3. p. 347-364. Miall. A. D.. 1984. Principles ofsedimemary basin analysis: Springer-Verlag. New York. 490 p. Powers, D. W.. and R. G. Easterling. 1982. Improved methodology for using embedded Markov chains to describe cyclical sediments: Jour. Sed. Pet.. v. 52. no. 3, p. 913924. Read. W. A.. 1969. Analysis and simulation of Namurian rocks oft:cntral Scotland using a Markov-process model: Jour. Math. Geology. v. I, no. 2. p. 199-219. Wells. N. A.. 1983, Transient streams in sand-poor redbeds: |-arly-Middle Eocene Kuldana Formation of northern Pakistan: Inter. Assot:. Sed. Spec. Publ.. no. 6, p. 3934O3. Wells. N. A.. 1984, Marine and continental sedimenhttion in the Early Cenozoic Kohat Basin and adjacent northwestern Indo-Pakistan: unpubl, dot:toral dissertation. Univ. Michig;,n, 465 p. Wells, N. A.. and Wahced, A., 1985, Randomness and disI'avorcd transitions m Markov chains (;.IBM.): Geol. Soc, America, Abst. with Program, v. 17, no. 7, p. 747.

APPENDIX

Markov AnalysL~" Program This program is written in Microsoft GW-BASIC version 3.2 and it runs perfectly on a Zenith Z-150 PC. It also runs in BASICA on an IBM XT, and, unless line I 110 is invoked, in "Sanyo-BASIC" on a Sanyo 555 PC. It therefore should Ix t:o,npatiblc reasonably with most BASICs on any IBM PC clone. Differences with other types of BASIC may involve the following: ( I ) different methods of combining statements in one line (instead of colons). (2) Different ways of routing output to printer or screen (here, LPRINT = printer, PRINT = screen). (3) Here, dimensioning is not required for variables with < 10 subscripts, but combination statements such as DIM T(20). U(20L B(15, 15) and LET C = X = O are not allowed. (4) Some types of BASIC will allow MAT statements (matrix manipulations), which would simplify parts of the program. (5) There may be different formats for PRINT USING statements. At 10cpi. an 80-column-wide printer can accommodate output for 10 facies: exceeding this results in a correct but barely interpretable mess. For easy reanalysis, it is convenient to store several data sets at the end of the program, but note that temporarily unwanted data statements must be inactivated by conversion into rein statements. This program is an update of an earlier, Sanyo-BASIC version handed out during the presentation of Wells and Waheed (1984).

Markov chain analysis

Program Listing i0 R E M *** P o w e r s & E a s t e r l i n g ' s M a r k o v A n a l y s i s , Nell Walls' v e r s i o n 2 20 R E M 30 P R I N T "Type 3 for full p r l n t - o u t of m a t r i c e s a n d f a v o r e d t r a n s i t i o n s ; " 40 P R I N T "2 for s e l e c t i o n of n o n - r a n d o m t r a n s i t i o n s only; or" 50 P R I N T "1 for b i n o m i a l p r o b a b i l i t i e s of s e l e c t e d c e l l s o n l y " 60 I N P U T Z4 70 P R I N T "IS Y O U R D A T A A L R E A D Y IN T H E PROGRAM, Y O R N " : I N P U T AS 80 R E M if p r o g r a m h a s data, t h e n it b r a n c h e s to d a t a - r e a d l n g r o u t i n e #2840 90 IF A S - "N" T H E N 130 i00 IF A S - "n" T H E N 130 Ii0 G O S U B 2860 120 G O T O 250 130 P R I N T "DO YOU W I S H T O I N P U T D A T A ON KEYBOARD, Y OR N " ; : I N P U T B$ 140 IF B$="Y" T H E N 190 150 IF B$-"y" T H E N 190 160 P R I N T " P R O G R A M W I L L N O W STOP. E D I T S A M P L E N A M E , " 170 P R I N T " N U M B E R O F FACIES, N A M E S OF FACIES, A N D D A T A A T P R O G R A M ' S END." 180 S T O P 190 P R I N T "TITLE " ; : I N P U T T$ 200 P R I N T "HOW M A N Y F A C I E S " ; : I N P U T N : G O S U B 2780 210 P R I N T "NAMES OF FACIES " ~ : F O R I=1 TO N : I N F U T F $ ( I ) : N E X T I 220 P R I N T "INPUT DATA, BY ROWS "; 230 F O R R=I TO N : F O R C-I TO N : I N P U T A ( R , C ) : N E X T C : N E X T R 240 REM . . . . . . . . . . . . . . . 250 R E M end input routine, s t a r t a n a l y s i s 260 L P R I N T " M A R K O V A N A L Y S I S " : L P R I N T T $ : L P R I N T 270 R E M c a l c u l a t i n g row t o t a l s F(R), c o l u m n t o t a l s G(C), & g r a n d total T 280 F O R R - 1 TO N 290 FOR C - i TO N 300 LET F(R) - F(R)+A(R,C) 310 LET G(C) - G ( C ) + A ( R , C ) 320 NEXT C 330 LET T=T+F(R) 340 N E X T R 350 REM . . . . . . . . . . . . . . . 360 R E M p r i n t i n g matrix; #420 - facies names, #450 - col. heads, 480 - row #, 370 R E M #500 - o b s e r v e d t r a n s i t i o n s (r to c), #520 - row totals, 380 R E M #560 & 570 - row t o t a l s & g r a n d total. 390 IF Z4<3 T H E N 590 400 L P R I N T " N u m b e r of facies ";N 410 L P R I N T : L P R I N T " F A C I E S : " : L P R I N T 420 F O R I-I TO N : L P R I N T U S I N G " # # " ; I ; : L P R I N T " ";F$(I):NEXT I 430 L P R I N T : L P R I N T : L P R I N T "RAW-DATA MATRIX":LPRINT 440 L P R I N T " #"1 450 FOR C-I TO N:LPRII4T " " ; : L P R I N T " # " ; : L P R I N T U S I N G "##"; C ; : N E X T C 460 L P R I N T 470 FOR R -i TO N 480 L P R I N T U S I N G "##"JR; 490 FOR C - I TO N 500 LPRINT " " ; : L P R I N T U S I N G "###"; A(R,C); 510 NEXT C 520 LPRINT " " ; : L P R I N T U S I N G "###"; F(R) 530 N E X T R 540 L P R I N T 550 L P R I N T " ";: 560 F O R C - 1 TO N : L P R I N T " ";:LPRINT USING "###";G(C);:NEXT C 570 LPRINT " ";:LPRINT USING "####";T 580 REM 590 R E M s t a r t of i t e r a t i o n s for row and c o l u m n f a c t o r s 600 R E M f i r s t i t e r a t i o n only, Iml, row f a c t o r s X(I,R) & t o t a l of r o w factors 610 L E T I - I 620 F O R R - I TO N 630 LET X(I,R) - F(R)/(N-I) 640 L E T T(I) - T(I) + X(I,R) 650 N E X T R 660 R E M c o l u m n factors Y(I,C) and total of c o l u m n f a c t o r s U(I) 670 R E M H - counter: if all d i f f e r e n c e s <0.01 t h e n H-O in line 930 680 R E M H2 - c o u n t e r to p r e v e n t t o o m a n y i t e r a t i o n s 690 L E T H2 - 0 700 LET H - 0 710 F O R C - 1 T O N 720 L E T Y(I,C) - G ( C ) / ( T ( I ) - X ( I , C ) ) 730 L E T U(I) - U(I) + Y(I,C) 740 IF I - 1 T H E N 770 750 IF Y(I,R) - Y(I-1,R) < .01 T H E N 770 CAGZO IS:|-K

151

152

N.A. WELLS

760 LETH-H÷I 770 N E X T C 7 8 O R E M R O W f a c t o r i t e r a t i o n s b e y o n d f i r s t i t e r a t i o n (I>l) 790 R E M R o w f a c t o r s , r o w f a c t o r t o t a l s , d i f f e r e n c e b e t w e e n i t e r a t i o n s , a n d 8 0 0 R E M c h e c k a g a i n s t f a i l u r e to c o n v e r g e (causes abandonment of p r o g r a m ) 810 L E T I - I + I 82O FOR R - 1 TO N 83O L E T X(I,R) - F ( R ) / ( U ( I - I } - Y ( I - I , R ) ) 840 L E T T(I) - T(I) + X ( I , R ) 85O IF X ( I , R ) - X(I-1,R) < .01 T H E N 8 7 0 LET H w H ÷ 1 860 870 NEXT R 88O IF I < (19 + H2) T H E N 930 PRINT H2+20;" ITERATIONS, NO CONVERGENCE YET" 89O P R I N T " H O W M A N Y M O R E I T E R A T I O N S DO Y O U W I S H T O T R Y ? " 90O 910 I N P U T H2 DIM T(H2):DIM U(H2):DIM X(H2,H2):DIM Y(H2,H2) 920 93O IF H > 0 T H E N 700 REM . . . . . . . . . . . . . . . 940 printing row and column factors 95O R E M 9 6 0 IF Z 4 < 3 T H E N 1 1 7 0 9 7 0 R E M 1 0 5 0 - i t e r a t i o n #, r f o r r o w # ; 1 0 7 0 & 1 0 9 0 - r o w f a c t o r s & t o t a l 9 8 0 R E M 1 1 1 0 - i t e r a t i o n #, c o l u m n #; 1 1 3 0 & 1 1 5 0 - c o l u m n f a c t o r & t o t a l 990 LPRIBT I000 LPRINT "ROW & COLUMN FACTORS" i010 LPRINT:LPRINT "Row Totals"; 1020 LPRINT:LPRINT R:LPRINT " " 1030 FOR R-1 TO N:LPRINT USING "######";F(R);:NEXT 1040 FOR G - i TO I 1050 LPRINT " ";G;"r"; 1060 FOR R - i TO N 1070 LPRINT USING "###.##";X(G,R); 1080 NEXT R 1090 LPRINT " " ; : L P R I N T U S I N G " # # # . # # " ; T(G) 1100 IF G - I T H E N 1 1 7 0 1110 LPRINT " ";G;"c" ; 1120 FOR C - 1 TO N 1130 LPRINT USING "###.##";Y(G,C) ; 1140 NEXT C 1150 LPRINT " " ; : L P R I N T U S I N G " # # # . # # " ; U(G) 1160 N E X T G 1170 REM 1180 REM printing expected, difference and chl-square matrices 1190 R E M B - e x p e c t e d , C - d l f f e r e n c e , D-chi-equare, 1200 REM 1280-asslgnment of 0 t o r-c; X2 - c h i - s q u a r e sum, E - a s y m m e t r y 1 2 1 0 R E M X3 - a s y m m e t r y c h l - s q u a r e sum; # 1 3 5 0 - 1 3 6 0 - d e a l i n g w i t h z e r o v a l u e s 1 2 2 0 R E M f o r d e g r e e s of f r e e d o m 1230 FOR R - 1 TO N 1240 FOR C - 1 TO N 1250 LET B(R,C) - X(I,R) * Y(I-I,C) 1260 LET C(R,C) - A(R,C) - B(R,C) 1270 L E T D ( R , C ) - C(R,C)^2/B(R,C) 1280 IF R - C T H E N 1 3 0 0 1290 L E T X2 - X2 + D(R,C) 1300 IF R <- C T H E N 1 3 7 0 1310 IF A ( R , C ) - A ( C , R ) < > 0 THEN GOTO 1330 LET E(R,C) - 0:GOTO 1340 1320 1330 LET E(R,C) - (A(R,C)-A(C,R))^2/(A(R,C)+A(C,R)) 1340 L E T X3 - X3 + E ( R , C ) 1350 IF A ( R , C ) O R A ( C , R ) <> 0 T H E N 1 3 7 0 1360 L E T Z0 ~ Z0 + 1 1370 NEXT C 1380 NEXT R 1390 REM . . . . . . . . . . . . . . . 1 4 0 0 R E M p r i n t i n g of f o u r m a t r i c e s 1410 IF Z4<3 THEN 1670 1420 FOR Z - 1 TO 4 LPRINT:LPRINT:LPRINT: ON Z GOTO 1440,1450,1460,1470 1430 1440 LPRINT "EXPECTED VALUES":GOTO 1480 L P R I N T "DIFFERENCE M A T R I X " : G O T O 1480 1450 LPRINT "SEQUENCE CHI-SQUARE MATRIX":GOTO 1480 1460 LPRINT "ASYMMETRY CHI-SQUAREMATRIX":GOTO 1480 1470 1480 LPRINT " " 1490 LPRINT " # "; FOR R-1 TO N:LPRINT USING "#######";R;:NEXT R:LPRINT 1500 1510 FOR R " 1 TO N 1520 L P R I N T U S I N G "##"; R; 1530 LPRINT " ";

Markov chain analysis FOR C " 1 TO N IF R " C T H E N 1620 ON Z G O T O 1 5 7 0 , 1 5 8 0 , 1 5 9 0 , 1 6 0 0 L P R I N T U S I N G "|###.##"~ B ( R , C ) ; : G O T O 1630 L P R I N T U S I N G "#||#.##"; C ( R , C ) ; : G O T O 1630 L P R I N T U S I N G "###|.##"; D ( R , C ) ; : G O T O 1630 IF R <" C T H E N 1620 L P R I N T U S I N G "####.##"; E ( R , C ) ; : G O T O 1630 LPRINT " - "; NEXT C LPRINT " " NEXT R NEXT Z REM . . . . . . . . . . . . . . . R E M DF" c h i - s q u a r e d e g r e e s of freedom; DOF- a s y m m e t r y d e g r e e s of freedom L E T DF - (N-1)^2 - N L P R I N T : L P R I N T "sequence c h l - s q u a r e sum -";X2, DF;" d e g r e e s of freedom" LET DOF - (N*(N-I)/2)-Z0 L P R I N T : L P R I N T " a s y m m e t r y c h i - s q u a r e - " ; X 3 , D O F ; " d e g r e e s of freedom" REM REM . . . . . . . . . . . . . . . IF Z4<2 T H E N 2340 R E M L i n e s 1 8 9 0 - 1 9 6 0 c o m p a r e c e l l s w i t h i n each row to i d e n t i f y p r i n c i p a l REM s u c c e e d i n g facies & to label l a r g e r v a l u e & c o l u m n #: M A X ( R ) & N3 R E M Lines 1990-2020 identify tied m a x i m u m values, 2 0 4 0 - 2 0 5 0 i d e n t i f y R E M w h e n m a j o r t r a n s i t i o n is not significant. R E M S u b s e q u e n t lines sort for o t h e r favored & d i s f a v o r e d t r a n s i t i o n s : R E M by d i s c a r d i n g s t r u c t u r a l l y e m p t y c e l l s (2180), p r e v i o u s l y s e l e c t e d REM c e l l s (2200), c h l - s q u a r e scores <2.71, & b i n o m i a l p r o b s <0.2 R E M #2270 p u l l s d i s f a v o r e d t r a n s i t i o n s & 2230 p u l l s f a v o r e d ones REM L P R I N T : L P R I N T : L P R I N T " P R I N C I P A L S U C C E S S I O N FOR EACH F A C I E S " L P R I N T "binomial p r o b a b i l i t i e s llsted at right give p r o b a b i l i t y " 1870 L P R I N T " ( c o n f i d e n c e level - 1 - level of significance) 1880 L P R I N T " of the o b s e r v e d f r e q u e n c y o c c u r r i n g by c h a n c e " : L P R I N T 1890 LET N 2 - N - 1 1900 FOR R-I T O N 1910 LET M A X ( R ) - C ( R , 1 ) : L E T N3-1 1920 FOR C-1 T O N2 1930 IF MAX(R) > C(R,C+I) THEN 1960 1940 LET M A X (R)-C(R, C+I) 1950 LET N3-C+I 1960 NEXT C 1970 LET N3 (R)-N3 1980 LET C - N3 1990 FOR J - 1 T O N 2000 IF J - N3 T H E N 2020 IF MAX(R) - C(R,J) THEN L P R I N T F$(R);" E Q U A L L Y T O "; F $ ( J ) ; " &" 2010 2020 NEXT J 2030 L P R I N T F$(R);" to "; F$(N3); IF D(R,N3) >- 5.38 THEN L P R I N T " (conf. > .99;";: G O T O 2070 2040 IF D(R,N3) >- 2.71 THEN L P R I N T " (conf..95-.99;";:GOTO 2070 2050 2060 LPRINT " (conf. < .95;"; 2070 G O S U B 2560 2080 LPRINT " BP~ ";:LPRINT U S I N G "#.####"; B P ; : L P R I N T "}" 2090 N E X T R 2100 L P R I N T : L P R I N T 2110 L P R I N T " O t h e r S i g n i f i c a n t F a v o r e d T r a n s i t i o n s " ; Significantly Disfavored Transltlons":LPRINT 2120 L P R I N T " 2130 L P R I N T "(conf. > ca. 0.95 & <20% p r o b a b i l i t y of s e l e c t i o n by c h a n c e ) " 2140 L P R I N T "(NB: b i n o m i a l p r o b a b i l i t i e s b e t w e e n .1 and .2 d o n ' t m e a n much)" 2150 L P R I N T 2160 FOR R-1 T O N 2170 FOR C-1 TO N 2180 IF R - C T H E N 2300 2190 IF D ( R , C ) < 2 . 7 1 T H E N 2300 IF C-N3(R) T H E N 2300 2200 2210 G O S U B 2560 2220 IF BP > .2 T H E N 2300 IF C(R,C) >O T H E N 2270 2230 LPRINT " ";F$(R);" to ";F$(C) 2240 2250 LPRINT " "; G O T O 2290 2260 2270 IF C ( R , C ) < O T H E N 2300 2280 L P R I N T F $ ( R ) ; " to ";F$(C) LPRINT " BP- " ; : L P R I N T U S I N G "#.####"; BP 2290 2300 NEXT C 2310 N E X T R

1540 1550 1560 1570 1580 1590 1600 1610 1620 1630 1640 1650 1660 1670 1680 1690 1700 1710 1720 1730 1740 1750 1760 1770 1780 1790 1800 1810 1820 1830 1840 1850 1860

153

154 2320 2330 2340 2350 2360 2370 2380 2390 2400 2410 2420 2430 2440 2450 2460 2470 2480 2490 2500 2510 2520 2530 2540 2550 2560 2570 2580 2590 2600 2610 2620 2630 2640 2650 2660 2670 2680 2690 2700 2710 2720 2730 2740 2750 2760 2770 2780 2790 2800 2810 2820 2830 2840 2850 2860 2870 2880 2890 2900 2910 2920 2930 2940 2950 2960 2980 2990 3000 3010 3020 3030 3040 3050 3060 3070 3080 3090 3100

N. A. WELLS REM . . . . . . . . . . . R E M 2 3 4 0 - 2 4 7 0 c a l c u l a t e d e s i r e d b i n o m i a l p r o b a b i l i t i e s f o r any cell P R I N T "Do y o u w a n t the b i n o m i a l p r o b a b i l i t y for a n y c e l l ? " P R I N T "Type y (lower case) If so" I N P U T C$ IF C$ <> "y" T H E N 3320 P R I N T " I n p u t R O W # of c e l I " : I N P U T R P R I N T "Input C O L U M N # of c e l I " : I N P U T C PRINT:LPRINT G O S U B 2490 R E M 2 4 3 0 - 2 4 6 0 send o u t p u t to s c r e e n if Just a n a l y z i n g p r o b a b i l i t i e s IF Z 4 < > 1 T H E N 2460 PRINT '*BINOMIAL P R O B A B I L I T Y F O R C E L L " ; R ; " , " ; C ; " - ";BP GOTO 2470 L P R I N T " B I N O M I A L P R O B A B I L I T Y FOR C E L L " ; R ; " , " ; C ; " - ";BP P R I N T "Do y o u w a n t a n o t h e r one? T y p e y (lower case) if y e s " : G O T O 2360 REM . . . . . . . R E M s t a r t of s u b r o u t i n e to c a l c u l a t e s u m m e d b i n o m i a l p r o b a b i l i t i e s , R E M from n to N for f a v o r e d t r a n s i t i o n s , or n to 0 if d i s f a v o r e d (Crc<0) R E M n-~ o b s e r v e d - A ( r , c ) ; N= max. p o s s i b l e - r o w t o t a l F(r) R E M 2 6 0 0 - 2 7 4 0 s u m b i n o m i a l p r o b a b i l i t i e s for at l e a s t obs. s u c c e s s e s R E M 2 6 2 0 - 2 6 6 0 c a l c u l a t e s f a c t o r i a l s for the s m a l l e r of (N-n)! & n! and REM for NI d i v i d e d by the l a r g e r of n! a n d (N-n) l. R E M 2 6 9 0 - 2 7 3 0 s u m b i n o m i a l p r o b a b i l i t i e s for e a c h i t e r a t i o n LET BP - 0 : L E T B P Z - 0 LET Y2 - F(R) LET Y1 - A(R,C) IF C(R,C) < 0 T H E N Y 1 - F ( R ) - A ( R , C ) FOR Z - Y1 TO Y2 LET A1 - O : L E T A2 - 0 I F Z - Y2 T H E N 2670 FOR J - 1 TO (F(R)-Z) LET A1 - A1 + L O G ( F ( R } - ( J - I ) } LET A2 - A2 + LOG(J) NEXT J P = B(R,C)/F(R) IF C(R,C) > 0 T H E N 2710 BPZ - (F(R)-Z)*LOG(P) + Z'LOG(I-P) + A1 - A2 G O T O 2720 BPZ - Z'LOG(P) + ( F ( R ) - Z ) * L O G ( I - P ) + A1 - A2 LET BP - BP + EXP(BPZ) IF B P > ( E X P ( B P Z ) * 1 0 0 0 ) T H E N 2750 NEXT Z RETURN REM REM . . . . . . . . . . . . . . . REM Dimensioning routine DIM T ( 2 0 ) : D I M U ( 2 0 ) : D I M M A X ( N ) : D I M N 3 ( N ) : D I M F ( N ) : D I M G(N) FOR C-I TO N : F ( C ) - O : G ( C ) - 0 : S ( C ) - O : T ( C ) - 0 : U ( C ) - 0 : M A X ( C ) - 0 : N 3 ( C ) - 0 : N E X T C DIM A ( N , N ) : D I M F$(N) D I M B ( N , N ) : D I M C ( N , N ) : D I M D(N,N) DIM X ( 2 0 , 2 0 ) : D I M Y(20,20) RETURN REM . . . . . . . R E M Data r e a d i n g routine: 2880- t i t l e & n u m b e r of f a c i e s REM 2890-dimensionlng, 2 9 0 0 - f a c l e s names, 2 9 2 0 - 2 9 6 0 - t r a n s i t i o n s READ T$,N G O S U B 2790 FOR F-1 TO N : R E A D F $ ( F ) : N E X T F R E M r e a d i n g o b s e r v a t i o n s for each cell (row r, c o l u m n c) FOR R - 1 TO N FOR C - 1 TO N R E A D A(R,C) NEXT C NEXT R R E M Data set: 1st, t i t l e of a n a l y s i s ; 2nd, n u m b e r of facies; R E M 3rd, facies n a m e s R E M 4th, t r a n s i t i o n f r e q u e n c y for e a c h cell (row r, c o l u m n c), BY R O W R E M several sets c a n be s t o r e d by r e l a b e l l i n g u n w a n t e d s e t s w i t h REM R E M DATA "CANT & W A L K E R ' S B A T T E R Y P O I N T S A N D S T O N E " R E M DATA 8 R E M DATA "SS s c o u r " , " A p e b b l e s " , " B t r o u g h s " , " C p l a n a r - t a b " R E M DATA "D small p l a n a r " , " E a s y m . s m . t r . " , " F r i p p l e s " , " G low a n g l e ss" R E H DATA 0 , 1 2 , 2 , 1 , 0 , 0 , 0 , 0 R E M DATA 2 , 0 , 6 , 3 , 0 , 1 , 0 , 1 R E M DATA 4 , 1 , 0 , 2 , 2 , 1 , 1 , 2 R E M DATA 0 , 2 , 4 , 0 , 1 , 0 , 0 , 0 R E M DATA 1 , 0 , 0 , 0 , 0 , 0 , 2 , 0

Markov chain analysis 3110 3120 3130 3150 3160 3170 3180 3190 3200 3210 3220 3230 3240 3250 3260 3270 3280 3290 3310 3320

DATA 0 , 0 , 0 , 1 , 0 , 0 , 1 , 0 DATA 2 , 0 , 0 , 0 , 0 , 0 , 0 , 1 DATA 5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 DATA "LOWER KULDANA FORMATION, EOCENE, PAKISTAN, all facies" DATA I0 DATA "Pebbles","Calc. Granules","Festoonm","Tabular Dunes" DATA "Plane Bed6","Ripple Troughs","Other Ripples" DATA "Non-red Clay","Red Clay","Channel Floor" DATA 0,4,8,2,0,0,1,0,2,0 DATA 0,0,30,5,5,4,3,8,28,0 DATA 1,8,0,7,22,17,4,3,20,3 DATA 2,1,1,0,4,3,2,0,5,1 DATA 0,1,6,0,0,12,6,5,19,3 DATA 0,1,2,0,5,0,4,0,29,2 DATA 0,2,0,3,0,0,0,2,14,0 DATA 0,2,0,0,1,2,0,0,20,1 DATA 0,12,4,0,5,6,2,12,0,104 DATA 14,51,27,4,5,8,2,0,0,0 RETURN END

155