Some remarks on the use of daily rainfall models

Journal of Hydrology, 36 (1978) 295--308

295

© Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands [4]

SOME REMARKS ON THE USE OF DAILY RAINFALL MODELS

T.A. BUISHAND

Department of Land and Water Use and Department of Mathematics, Agricultural University Wageningen, Wageningen (The Netherlands) (Received May 20, 1977; accepted for publication August 11, 1977)

ABSTRACT Buishand, T.A.,1978. S o m e r e m a r k s o n t h e use of daily rainfall models. J. Hydrol.,36: 295--308. Features of daily rainfall processes are described using data from different sites of the world. The process of rainfall occurrence is modelled by an alternating renewal process or by a Markov chain. It is shown that rainfall amounts within a wet spell are often neither independent nor identically distributed. There are, however, features of the rainfall process which are hardly sensitive to the choice of the process for the occurrence of wet days or to some assumptions about the behaviour of the rainfall amounts, for example the distribution of 30-day totals and the distribution of monthly extremes.

INTRODUCTION

In the last two decades techniques have been developed for generating daily rainfall sequences. Synthetic rainfall data can be used as the input to a hydrological system to obtain a better insight into the system performance. Because of the complex nature of the rainfall process a rather sophisticated model is needed to describe the many different aspects of this process well. For solving practical problems one may ask, however, whether one should use a complex model or only a simplified version of it. To answer this question one needs to know which features of the model affect the solution of the problem at hand. For instance, it may be possible that the solution is hardly sensitive to the choice of the marginal distribution for the rainfall amounts. Then it is not necessary to spend much time in modelling this distribution well. The first part of this paper describes features of daily rainfall sequences using data from different climatic areas. The importance of the various aspects of the daily rainfall process is examined in the second part.

296 FEATURES OF DAILY RAINFALL SEQUENCES Most techniques for generating daily rainfall sequences use a separate process for the rainfall occurrence (denoted shortly as wet--dry process) and another process for the rainfall amounts on wet days. To use such techniques one needs a suitable definition of a wet day. When all days with a positive rainfall a m o u n t are called wet there may be some problems since small values can also be caused by fog or dew. On the other hand, small rainfall amounts are often registered as zero. It is possible, therefore, that a change of observer also results in a change of the wet--dry series (Buishand, 1977). To circumvent these problems days with a rainfall a m o u n t smaller than some threshold will also be considered as dry days. T h e p r o c e s s f o r the rainfall o c c u r r e n c e

For modelling wet--dry sequences use is frequently made of one of the following two models (Lowry and Guthrie, 1968; Cole and Sherriff, 1972). (1) A l t e r n a t i n g r e n e w a l processes. This model assumes that lengths of successive wet and dry spells are independent. A wet spell is defined here as a sequence of wet days, on each side bounded by a dry day. A dry spell can be defined analogously. A suitable distribution to fit the lengths of wet and dry spells is the truncated negative binomial distribution (TNBD). The TNBD with parameters p and r can be defined by its probability function: P(X=k) =

k

I 1 _pr .

.

w i t h 0 < p ~ < 1, r ~ > - l a n d q =

k = 1, 2, . .

1-p.

For r = 1 one gets the geometric distribution (GD): p(X=k ) = pqk-1,

k = 1, 2 , . . .

and for r -~ 0 one gets the logarithmic series distribution (LSD): P ( X = k ) = a q k /k,

k = 1, 2 , . . .

with a = - 1 / l n p. The last two distributions have widely been used to fit the distribution of weather spells. Since the wet--dry process shows seasonal variation its parameters were estimated for each m o n t h separately. For Winterswijk (The Netherlands), 1908--1973, Table I shows some features of distributions fitted to lengths of spells, beginning in January. The threshold defining a wet day is 0.8 mm in this table. When parameters for each month are estimated separately one gets 12 X 2 estimates for each type of spell if a TNBD is used. To reduce this number it

297

TABLE I

Features of distributions fitted to lengths of weather spells for Winterswijk (The Netherlands) 1908--1973, January Type of spell

Dry

Wet

Distribution

TNBD LSD GD TNBD LSD GD

Mean (days)

3.55 3.57 3.55 2.54 2.55 2.54

Standard deviation (days)

3.75 4.29 3.01 2.07 2.60 1.98

Coefficient of skewness

2.84 3.42 2.03 2.22 3.21 2.06

Fraction of spells with length (days): 1

>6

>10

0.363 0.409 0.282 0.408 0.490 0.393

0.145 0.145 0.137 0.054 0.065 0.037

0.056 0.072 0.050 0.009 0.021 0.007

was proposed to keep the parameter r constant throughout the year. Estimates of the 12 + 1 parameters were obtained by an approximate maximum likelihood method. The estimate of r was 0.246 for dry spells, and 0.656 for wet spells. From Table I it is seen that the LSD has the longest tail and the largest fraction of 1-day spells; the GD has the shortest tail and the smallest fraction of 1-day spells. The general TNBD lies between these two cases. The fit of the TNBD is good both for wet and dry spells. Special cases of this distribution provide a less satisfactorily fit. For dry spells the ×Z-test of goodness of fit leads in nearly all months to significant values at the 5% level when a GD is assumed; for wet spells the LSD is often rejected at the 5% level (Buishand, 1977).

(2) Markov chains. Here the assumption is made that the probability of a wet or dry day depends on the situation of a certain number (the order of the chain) of previous days. Estimates of the parameters of a Markov chain can be obtained more easily than for an alternating renewal process. Besides this, the generation of synthetic sequences is simpler for a Markov chain. A disadvantage of this model may be that sometimes the order of the chain (and therefore the number of parameters) must be quite large to obtain a reasonable fit to lengths of weather spells. A comparison between an alternating renewal process and a Markov chain is made in Table II. This table gives the distribution of the lengths of dry spells for the historic series of Pasar Minggu (Indonesia) and for some synthetic sequences. The generation models considered are: an alternating renewal process with a LSD for lengths of wet spells and a TNBD for lengths of dry spells ( L S D - T N B D process), a first-order Markov chain (Markov-1) and a second-order Markov-chain (Markov-2).

298

T A B L E II T h e distribution of the lengths of dry spells for Pasar Minggu (Indonesia)

Length

Historic

(days)

series

LSD--TNBD

Markov-1

Markov-2

S~

S2

S~

S2

S~

S~

1 2 3 5 10 20 30

0.328 0.513 0.639 0.783 0.916 0.973 0.989

0.346 0.537 0.644 0.781 0.915 0.973 0.990

0.323 0.516 0.626 0.810 0.906 0.973 0.988

0.269 0.451 0.583 0.745 0.906 0.979 0.995

0.263 0.442 0.574 0.745 0.913 0.982 0.996

0.324 0.476 0.593 0.741 0.899 0.976 0.991

0.316 0.471 0.581 0.739 0.898 0.976 0.994

3rd longest spell 2nd longest spell longest spell

87 115 156

70 77 83

72 88 91

49 49 49

47 53 56

51 52 52

54 59 59

Cumulative relative frequencies and the lengths of the three longest spells are given for the historic series and for two synthetic sequences (S 1 and S 2 ) of different types of wet--dry processes.

The length of the historic series and of the synthetic sequences is 69 years. The threshold defining a wet day is 3.0 mm. This rather large threshold is necessary to obtain a homogeneous wet--dry series. Seasonal variation of the wet--dry process was incorporated by changing monthly its parameters. The first-order Markov chain has the shortest tail and its fraction of 1-day spells is too small. Since for this process lengths of weather spells are distributed geometrically (GD--GD process) this distribution has the same shortcomings here as for weather spells of Winterswijk. When a second-order Markov chain is used the fit is improved markedly. Yet, the right tail of this model is shorter than that of the LSD--TNBD process. The last model is therefore preferable. For rainfall stations with a long dry season it is often difficult to fit distributions to lengths of dry spells during the dry season and during the transition periods. An example is given in Table III, which shows the frequency distribution of the lengths of dry spells beginning in October for the rainfall series of New Delhi (1901--1970). The length of each synthetic series is 70 years. The threshold defining a wet day is 0.2 mm in Table III. The frequency distribution shows two maximums; namely there are many relatively short spells, but there are also some very long spells since it may happen that the dry-monsoon begins in October. To fit a theoretical distribution one needs a large number of parameters in this case. But a first-order Markov chain provides a reasonable fit here.

III

13

9.8

8 10 10 9 12

20

8

7.4

4 9 7 11 6

30

2

6.2

6 2 5 10 8

40

The synthetic data are based on a seasonal first-order Markov chain.

14

H i s t o r i c series

12 10 12 22 15

14.2

S1 S2 S3 S4 Ss

A v e r a g e o f s y n t h e t i c series

S y n t h e t i c series

10

Class w i t h u p p e r b o u n d ( d a y s )

0

7.6

6 7 11 8 6

50

5

4.8

4 4 7 2 7

60

3

6.8

5 12 7 5 5

70

F r e q u e n c y d i s t r i b u t i o n o f t h e l e n g t h s o f d r y spells f o r N e w D e l h i ( I n d i a ) 1 9 0 1 - - 1 9 7 0 , O c t o b e r

TABLE

4

6.0

6 6 7 7 4

80

6

4.2

4 3 5 3 6

90

6

3

3.0

3 5 1 4 2

2 4 2 2 2 2.4

oo

100

64

72.4

60 72 74 83 73

Number of spells

~D

t~

300

The behaviour o f rainfall amounts on wet days (1) Serial correlation. Often a small but significant dependence exists between rainfall amounts within a wet spell. For instance, for rainfall series in The Netherlands the first-order autocorrelation coefficient ranges from 0.10 to 0.15 during winter. A bout the same values were found in the wet season for rainfall series from India and Indonesia. But there is no evidence for a serial correlation for rainfall series during summer in The Netherlands. (2) Marginal distribution. The frequency distribution of rainfall amounts on wet days is usually J-shaped. Since a wet day was defined as a day with a rainfall a m o u n t above some threshold 6 the lower limit of the distribution for the rainfall amounts is 5, or more precisely 5 minus half the unit of measurement. The last value was subtracted from the original rainfall amounts and then the gamma distribution was fitted. A justification of the use of the gamma distribution is given in Table IV. For most rainfall stations the ratio of the coTABLE IV Ratio of the coefficient of skewness to the coefficient of variation for rainfall amounts Rainfall station

Length (years)

~ (mm)

Ratio

Winterswijk (The Netherlands) Paramaribo (Suriname*) Pasar Minggu (Indonesia) Calcutta (India) Bangalore (India) New Delhi (India)

66 70 69 70 68 70

0.8 0.3 3.0 0.2 0.2 0.2

2.23 2.12 1.96 2.09 2.05 2.22

The threshold defining a wet day is denoted by 6. Rainfall amounts are shifted to the interval (0,~). *Formerly Dutch Guiana.

efficient o f skewness to the coefficient of variation is close to its theoretical value 2. Because of seasonal variation the ratios in Table IV were obtained by averaging twelve m o n t h l y estimates.

(3) The mean rainfall a m o u n t on a wet day related to its position in a wet spell. When working with rainfall observations at discrete time points one has the problem that the first and last day of a wet spell contain a part of the preceding or following dry spell. Therefore one may ask w het her the mean of such day differs from~ that of the ot he r days. To test this assumption three different types of wet days were distinguished, namely (a) solitary wet days;

301

(b) wet days at one side bounded by a wet day and at the other side bounded by a dry day; (c) and wet days at each side bounded by a wet day. The rainfall amounts on these different wet days are denoted by type 0, 1, and 2, respectively; so a t y p e - / a m o u n t denotes a rainfall amount on a wet day with i adjacent wet days. Estimates of the mean are given in Table V. From this table TABLE V Means of rainfall amounts of different types Rainfall station

Winterswijk Paramaribo

Pasar Minggu Calcutta Bangalore New Delhi

Period

Dec.--Feb. Jun.--Aug. Dec.--Apr. May--Jun. Sep.--Nov. Jam--Feb. Jun.--Aug. Jun.--Aug. Jun.--Aug.

6 (ram)

0.8 0.8 0.3 0.3 0.3 3.0 0.2 0.2 0.2

Mean (mm)

Critical level of the F-test

type 0

type 1

type 2

3.34 5.56 4.43 6.48 7.72 19.53 11.19 6.60 8.66

4.48 6.30 7.65 11.77 8.07 19.63 11.64 6.64 14.54

5.72 6.43 12.34 12.19 8.83 21.16 16.02 7.75 19.40

0.000 0.056 0.000 0.002 0.133 0.230 0.000 0.029 0.000

A type i amount is a rainfall amount on a day with i adjacent wet days. The last column gives the critical level of the F-test for equality of means of type 0, 1 and 2 amounts (a small value supports the reality of differences in mean).

it is seen that the mean increases with the number of adjacent wet days. The largest differences are found for Winterswijk (Dec.--Feb.), Paramaribo ( D e c . Apr.) and New Delhi. There is also an obvious seasonal variation in the differences. The F-test rejects nearly always the assumption of equal means at the 5% level. There may be some criticism about the applicability of the F-test since the data are highly skewed. However, nearly the same results were obtained with the cube root (normalizing transformation) of the data (Buishand, 1977). THE IMPORTANCE OF SOME FEATURES OF DAILY RAINFALL SEQUENCES

A good description of all features of the rainfall process needs a rather complex model. Such a model, however, should incorporate a large number of parameters. Moreover the generation of synthetic data can be very cumbersome. For instance, generating synthetic sequences with serially correlated non-identically distributed rainfall amounts within a wet spell is rather tedious. As indicated in the introduction, one should verify how the various aspects of the model affect those features of the rainfall process which are relevant in a particular planning situation. First the choice of the wet---dry process is examined; and then the model for the rainfall amounts.

302

Re[ freq 1.0

-

-

i~-~.

-~

//"

.//"

,//

"Historic data - TNBD-TNBD process ..... GD-GD process ....... GD-LSD process

,~/ ,~.'

---

~-~

I

I

-

I

10

r

~

--

i

l

I

i

r

-

20

days Fig.1. Cumulative frequencies of the number of wet days in a 30-day period for the historic series of Winterswijk (Dec.--Feb.)andfor some theoretical wet--dry processes. Smooth curves are drawn through the theoretical values at the integers.

The choice of the process for the occurrence of wet days (1) The distribution o f the number of wet days in a 30-day period. Fig.1 gives the cumulative distribution function (CDF) of the number of wet days (h = 0.8 mm) in a 30-day period for Winterswijk (Dec.--Feb.). Though there may be some differences in the distributions of lengths of wet and dry spells for these processes (see Table I) the differences between the CDF's in Fig.1 are only very small. Besides this, there is a nice correspondence with the empirical distribution function. For the probability distribution of the number of wet days in a n-day period Gabriel (1959) developed a formula for a first-order Markov chain. Elliott (1965) gave this probability distribution for a renewal process (this is an alternating renewal process in which at least one type of spell is geometrically distributed). For an extension of Elliott's method to the general alternating renewal process the reader is referred to Buishand (1977). (2) A first passage time problem. For rainfall series in Indonesia Schmidt and Van der Vecht (1952) considered the cumulative rainfall a m o u n t after August 31 to describe the transition from the dry to the wet season. The date on which this cumulative rainfall amount exceeded the 350-mm level was taken

303

as the beginning of the wet season. The value of 350 mm was chosen because then the soil is, in general, sufficiently wet for preparing the seed beds for the rice crop in the wet season. It may be worthwhile to investigate the sensitivity of the so defined beginning date to the type of wet--dry process. Fig.2 compares the CDF of the beginning date of the historic series of Pasar Minggu with those of synthetic sequences. In one case only rainfall amounts were generated, using the historic wet--dry series. In the other cases the synthetic wet--dry sequences mentioned in Table II were taken. Rainfall amounts within a wet spell were assumed to be independent and identically distributed (i.i.d.). JiShnk's (1964) algorithm was used to generate gamma-distributed rainfall amounts. The CDF's in Fig.2 do not always reach the value 1, since it may happen that the 350-mm level is reached after December 31. Though the best fit is obtained if the historic wet--dry sequence is used the fit of the other cases is only slightly worse. Obviously the choice of the type of wet--dry process is relatively unimportant here. But there are features of the rainfall process which are sensitive to the type of wet--dry processes. For instance, the standard deviations of annual totals were 220 and 282 mm for the synthetic data based on the Markov-1 process; for the synthetic data based on the historic wet--dry series these values were 450 and 476 mm. The last two values correspond quite well to the estimated standard deviation (411 mm) of the historic series.

The choice of the model for the rainfall amounts (1) The distribution of 30-day totals. In Table V significant differences were found between the means of different types of rainfall amounts. One may ask whether it is important to include this phenomenon in the model. For the rainfall station of Winterswijk a comparison was made between a model with i.i.d, rainfall amounts within a wet spell, and a model with independent but non-identically distributed rainfall amounts on the basis of the CDF of 30-day totals. To obtain an insight into the shape of the CDF of 30-day totals for a model with different types of rainfall amounts five synthetic sequences were generated; for a model with i.i.d, rainfall amounts it is not difficult to obtain the CDF by numerical methods. The length of each synthetic sequence was the same as that of the historic series (66 years). The CDF of n-day totals for a model with i.i.d, rainfall amounts can be calculated as follows (Todorovic, 1970). Let Pn (k) be the probability of k wet days in a n-day period, and let H k (x) be the CDF of the sum of k i.i.d. rainfall amounts, then the following expression holds for the CDF F n (x) of n-day totals: t~

F n (x) =pn (0)+ ~

k=l

Pn (k) H k (x)

304 The calculation of the probabilities Pn (k) was discussed earlier; for gammadistributed rainfall amounts the sum of k i.i.d, values is again gamma-distributed. That is the calculation of H k (x) involves the calculation of an incomplete gamma function. Fig.3 shows the CDF's of the different models together with the empirical distribution function. From that figure it is seen that the CDF of a simple Rel freq.

1.0~ Historic

w e t - d r y series

m¢.

L

+ "o

o$÷

o+

0.5-

o,÷ o~÷ 0

•

O44 Oe •

÷

O

September I

I November

October

Oecember

Rel freq. 1.0-

+ o

1

LSD-TNBD process

°+°

+g +o o

+ +

0.5.

+

o °

•

+

oo

**o• O

•

O o . 1•

+

O

Septemberl÷ October

I November

December

]

305 Rel freq.

1.0~-- Markov- 1 process

o<-

/

0

+ +

6

•

o+• + o+ 40.e o

+

0.5 •o4+

•

o4-

4*o

o

O+

September i

October +~[ November 7

December

Rel. freq. 1.0--

i

2 4. °

Markov- 2 process

•

%

0.5•

+o

+o

•

o

• o ,o+ O• •

O

~

September I

• Historic data +$1 oS2

+

_

October

I November

I December

I

Fig.2. Cumulative frequencies of the beginning date of the wet season for the historic series of Pasar Minggu and for synthetic sequences (S~ and S 2) based on various types of wet--dry processes. The first day of the wet season is defined as the day for which the cumulative rainfall amount after August 31 exceeds the 350-ram level.

model with i.i.d, rainfall amounts fits poorly. The fit is hardly improved, however, when a model with different types of rainfall amounts is taken. The CDF's in Fig.3 are based on TNBD's for lengths of wet and dry spells. Since the type of wet--dry process hardly influences the probabilities Pn (k), see Fig.l, the results are approximately the same when a GD--LSD or a GD-GD process is used.

306 Re[ freq. 10

•

:

•

• J • • S• • • ~//@ •

35 / /

w •

•

•

I/"

/

;

°/

Historic data

•

"//

Synthetic data .....

/

Theoretical

/

b

5O

~

T

100

150 mm

Fig.3. Cumulative frequencies of 30-day totals for the historic series of Winterswijk (Dec.Feb.) and for five synthetic sequences. For the synthetic sequences the smallest and the largest value are given at a certain rainfall depth. Points of the synthetic sequences are omitted when they coincide with values of the historic series. The theoretical distribution function is based on a simplified model (i.i.d. rainfall amounts within a wet spell). A s h o r t c o m i n g o f the models is t h a t t h e y u n d e r e s t i m a t e the variance. T h e t h e o r e t i c a l variance o f a m o d e l with i.i.d, rainfall a m o u n t s is 639 mm2; for a m o d e l with d i f f e r e n t t y p e s o f rainfall a m o u n t s this value is 759 m m 2. These values are m u c h smaller t h a n the variance e s t i m a t e d f r o m the historic series (1,167 mm~). F o r m o d e l s with serial c o r r e l a t i o n b e t w e e n rainfall a m o u n t s within a wet spell the variance is o n l y slightly larger than t h a t o f models witho u t serial c o r r e l a t i o n (Buishand, 1977).

(2) Extreme values. Consider daily rainfall observations in a particular m o n t h , for e x a m p l e January. F o r each y e a r the largest rainfall a m o u n t in t h a t m o n t h is selected. It is a c o m m o n practice t o fit some distribution t o these e x t r e m e values, f o r e x a m p l e the G u m b e l distribution. But it is also possible t o obtain the d i s t r i b u t i o n o f the e x t r e m e values by n u m e r i c a l m e t h o d s for rainfall m o d e l s with i.i.d, rainfall a m o u n t s within a w e t spell. F o r the calculation o f the CDF Gn (x) o f the m a x i m u m in a n-day period the following expression can be used ( T o d o r o v i c , 1970): n

V n(x)=pn (0)+ k~=lPn (k) {H1 ( x ) } k

307

TABLE VI Monthly extremes (ram) for different return periods Rainfall station

Month

Return period (years) 10

2

G

M

G

100 M

G

M

Winterswijk

Dec. Jan. Feb. Jul.

14 12 11 19

15 13 12 19

24 23 21 33

24 21 21 32

37 35 33 51

37 32 32 49

Pasar Minggu

Jan.

60

60

Feb. Jul. Aug.

54 32 23

54 35 29

100 89 76 55

99 90 75 59

150 132 132 95

146 135 124 98

G is based on a fitted Gumbel distribution and M is the theoretical value for a model with i.i.d, rainfall amounts within a wet spell.

where H1 (x) stands for the C D F o f the rainfall a m o u n t s on w e t days. Table VI c o m p a r e s t h e o r e t i c a l values o f the simple rainfall m o d e l with t h o s e o b t a i n e d b y fitting a G u m b e l d i s t r i b u t i o n for the rainfall series o f Winterswijk and Pasar Minggu. Estimates o f the p a r a m e t e r s of the G u m b e l d i s t r i b u t i o n were o b t a i n e d b y the m e t h o d o f m a x i m u m likelihood. There is a nice corres p o n d e n c e b e t w e e n the t w o values, indicating t h a t it is n o t necessary t o take a m o r e a d v a n c e d m o d e l for the b e h a v i o u r o f rainfall a m o u n t s here. The t h e o r e t i c a l values in Table VI are based on a L S D - - T N B D process for Pasar Minggu and a T N B D - - T N B D process for Winterswijk. F o r the last station nearly t h e same results were o b t a i n e d with a G D - - G D process (firsto r d e r Markov chain). This was t o be e x p e c t e d since there is h a r d l y a n y difference b e t w e e n the probabilities Pn (k) for these processes, see Fig.1. CONCLUDING REMARK In the f o r e g o i n g it was d e m o n s t r a t e d t h a t s o m e features o f the rainfall process are h a r d l y sensitive t o the choice o f the rainfall model. When these features are i m p o r t a n t in a particular planning situation it does n o t pay to use a c o m p l i c a t e d model. However, o n l y a few features were s t u d i e d in this paper. More research is n e e d e d t o test the p e r f o r m a n c e o f existing rainfall models.

308

ACKNOWLEDGEMENTS This p a p e r d e s c r i b e s t h e m a i n ( a n d s o m e a d d i t i o n a l ) r e s u l t s o f t h e r e s e a r c h d o n e f o r t h e Ph.D. t h e s i s o f t h e a u t h o r . This thesis was p r e p a r e d u n d e r supervision o f t h e P r o f e s s o r s W.H. van d e r Molen a n d P. van d e r Laan. T h r o u g h o u t t h e s t u d y Dr. M.A.J. van M o n t f o r t m a d e m a n y h e l p f u l s u g g e s t i o n s a n d comments. The r a i n f a l l d a t a in this r e s e a r c h were p r o v i d e d b y t h e M e t e o r o l o g i c a l I n s t i t u t e s o f G r e a t Britain, I n d i a , I n d o n e s i a a n d T h e N e t h e r l a n d s .

REFERENCES Buishand, T.A., 1977. Stochastic modelling of daily rainfall sequences. Commun. Agric. Univ. Wageningen, Wageningen, 77-3. Cole, J.A. and Sherriff, J.D.F., 1972. Some single- and multi-site models of rainfall within discrete time increments. J. Hydrol., 17: 97--113. Elliott, E.O., 1965. A model of the switched telephone network for data communications Bell Syst. Tech. J., 44: 89--109. Gabriel, K.R., 1959. The distribution of the number of successes in a sequence of dependent trials. Biometrika, 46: 454--460. J6hnk, M.D., 1964. Erzeugung von betaverteilten und gammaverteilten Zufallszahlen. Metrika, 8: 5--15. Lowry, W.P. and Guthrie, D., 1968. Markov chains of order greater than one. Mort. Weather Rev., 96: 798--801. Schmidt, F.H. and Van der Vecht, J., 1952. East monsoon fluctuations in Java and Madura during the period 1880--1940. Djawatan Meteorol. Geofis., Verh., No. 43. Todorovic, P., 1970. On some problems involving random number of random variables. Ann. Math. Stat., 41: 1059--1063.