{'fraction_non_alphanumeric': 0.04478453134483656, 'fraction_numerical': 0.02850515910864476, 'mean_word_length': 3.953613228819929, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling Correlated Bernoulli Data Part I: Theory and Run Lengths

Binary data are very common in many applications, and are typically simulated independently via a Bernoulli distribution with a single probability of success. However, this is not always the physical truth, and the probability of a success can be dependent on the outcome successes of past events. Presented here is a novel approach for simulating binary data where, for a chain of events, successes (1) and failures (0) cluster together according to a distance correlation. The structure is derived from de Bruijn Graphs - a directed graph, where given a set of symbols, V, and a 'word' length, m, the nodes of the graph consist of all possible sequences of V of length m. De Bruijn Graphs are a generalisation of Markov chains, where the 'word' length controls the number of states that each individual state is dependent on. This increases correlation over a wider area. To quantify how clustered a sequence generated from a de Bruijn process is, the run lengths of letters are observed along with run length properties.


Introduction
Although binary data are often considered a simple concept, they are typically difficult to model in statistics. Traditionally, binary random variables are modelled independently via a Bernoulli distribution with a single probability of success (π). However, in many underlying processes, true independence of outcomes is not guaranteed. The correlation is often dependent on the distance between outcomes; given a series of binary trials in a one dimensional space it is often the case that the successes in nearby trials are more highly correlated than trials far away. For example, in ecology [Pielou, 1984[Pielou, , 1969 we may be modelling the absence or presence of individuals in a population such as trees or bacteria (denoted 0 or 1 respectively). Clearly, there is a much higher chance of observing trees/bacteria close to spatial locations that are known to already contain trees/bacteria.
The aim of this work is to model binary data from a correlated process such that there is a spatial (or temporal) varying probability of success; the probability of obtaining a 1 is higher having observed other 1's (and similarly for 0's). It is then expected that the successes (1's) will cluster (or anti-cluster in the case of negative correlation). Methods such as logistic regression classification or generalised linear models [Diggle et al., 1998, Chang et al., 2016 are typically used in these scenarios. They allow the marginal probability (π) of a 1 to vary smoothly in space, and therefore at any point x we have a marginal Bernoulli distribution with probability of success equal to π(x). However, again only having information regarding the marginal probability of a success at each trial causes sequence samples to neglect any form of spatial correlation.
Correlated multivariate Bernoulli distributions have been discussed in literature, where it is possible to draw a sample of length n binary random variables jointly. Teugels [1990] introduces a multivariate Bernoulli distribution, such that given a sequence of n Bernoulli trials, the joint probability of each of the possible 2 n sequences of 0's and 1's is specified. The distribution is then parametrised in terms of the central moments so that each of these sequences is further dependent on 2 n−1 parameters. Although the work by Teugels [1990] successfully defines a multivariate Bernoulli distribution, a large number of parameters are required for even short sequences of Bernoulli trials. Secondly, there is no metric to control the spread of correlation across the binary variables. Similar approaches have been seen in both Tallis [1962] and Tallis [1964].
Graphical models have also proven to be useful in simulating correlated variables. Similar to the structures used in Teugels [1990], Dai et al. [2013] use a multivariate Bernoulli distribution to estimate the structure of graphs with binary nodes. The binary nodes are described as random variables, where pairwise correlation is defined in terms of the edges of the graph. Variables are conditionally independent if the associated nodes are not linked by an edge. The authors further extend the model to take account of higher order interactions (not just pairwise interactions) to allow a higher level of structure to be incorporated. Another alternative is to use a Markov chain [Billingsley, 1961, Ching and Ng, 2006, Isaacson and Madsen, 1976, so that each variable is dependent on the previous variable in a graph structure. However, we find that the Markov property by itself is fairly limiting in the possible range of correlation structures. We also point readers to the image processing literature where a similar problem arises when restoring fuzzy images [Besag, 1986, 1974, Abend et al., 1965, Bishop, 2006, Li et al., 2016. This also includes literature relating to cellular automata [Wolfram, 2002, Agapie et al., 2014, 2004.
In this paper, we propose a novel framework for simulating binary trials from a correlated multivariate Bernoulli process. We define this process, the de Bruijn process (DBP), so named since we incorporate a distance correlation between variables using the structure of de Bruijn graphs [De Bruijn, 1946, Good, 1946, Golomb, 1967, Fredricksen, 1992. These are directed graphs, which have nodes consisting of m symbols. We denote the symbols the 'letters' of the de Bruijn graph, and the sequence of these letters at each node a de Bruijn 'word' of length m. We let the set of letters be {0, 1} and incorporate a generalised Markov property [Billingsley, 1961, Ching and Ng, 2006, Isaacson and Madsen, 1976 within the graph. The length m sequences of binary letters then become the states (or words) of a Markov chain where the spread of correlation is controlled by varying m. Figure 1 shows two de Bruijn graphs with the set of letters {0, 1}, and word lengths m = 2 (top) and m = 3 (bottom). The nodes consist of all the possible length two or three sequences of 0's and 1's where there are two edges coming in and out of each node to give a total of 2m nodes and 2m + 1 edges. By defining a distinct probability for each possible transition that can occur in the process (each of the edges), we can control the amount of clustering of each letter when simulating from the generalised Markov chain. Other work involving de Bruijn graphs can be seen in Hauge and Mykkeltveit [1996], Hauge and Helleseth [1996] and Hunt [2002]. With a DBP, we are able to generate a variety of sequences that can be very clustered together for 0's or 1's, or alternatively, can be very structured to continuously be alternating (anti-clustered). We hence analyse the run lengths of letters to help quantify how clustered a sequence is likely to be. In this paper, we define a run length, R, to be the number of consecutive 1's (or 0's) in a row bounded by a 0 (or a 1) at both ends, and use this to calculate a run length distribution. The distribution gives the probability of a run of length n for any n ∈ N + in terms of the word length m and the transition probabilities (as given by the edges in Figure 1). From this we can then calculate expected run lengths, variance of run lengths and generating functions. As previously mentioned, in ecology there is a large amount of research in the measurement of aggregation (e.g. level of clustering of trees in a one dimensional direction). Currently, this often relies on only testing data for randomness instead of levels of clustering. Hence, by analysing expected run lengths, it will be possible to make better predictions on levels of aggregation. In this paper we apply our run length distribution for a DBP to precipitation data collected daily from a station in Eskdalemuir, UK. We translate the data to be binary such that 1 is recorded for when precipitation is present and 0 otherwise. To remove any seasonal effects, we separate the data into seasons and analyse separately. In each season, we successfully show that the data favours a DBP structure rather than an independent Bernoulli equivalent.
The remainder of this paper is organised as follows. De Bruijn graphs are explained further in Section 2. This includes details and relationships that come from analysing the stationary distribution. By introducing correlation between binary variables in a sequence, we can control how clustered the 0's and 1's are. Hence, Section 3 gives details of a run length distribution of consecutive numbers of 1's in the sequence. This includes definitions for the expectation, variance, generating functions. Examples including an application to precipitation data is given in Section 4. Finally, Section 5 concludes with details of future work.

De Bruijn Graphs
De Bruijn Graphs [De Bruijn, 1946, Good, 1946, Golomb, 1967, Fredricksen, 1992 are directed graphs consisting of overlapping sequences of symbols. Given a set of s symbols, V = {v 1 , ..., v s }, the vertices or nodes of the graph consist of all the possible sequences of V . Each graph has s m vertices, where m is the length of each possible sequence given the set of symbols, V . The possible nodes are as follows: Edges in de Bruijn graphs are drawn between node pairs in such a way that the connected nodes have overlaps of m − 1 nodes. An edge is created by removing the first symbol from the node sequence, and adding a new symbol to the end of the sequence from V . Thus, from each vertex, (v 1 , ..., v m ) ∈ V m , there is an edge to vertex (v 2 , ..., v m , v) ∈ V m for every v ∈ V . There are exactly s directed edges going into each node and s directed edges going out from each node. We denote the symbols, v, the 'letters' of the de Bruijn graph, and the sequence of these letters at each node a de Bruijn 'word' of length m.
By travelling along a path through the de Bruijn graph, chains of letters are generated such that each letter is dependent on the m letters that come before [Fredricksen, 1992, Ayyer and Strehl, 2011, Rhodes et al., 2017. Hence, by altering m, we are able to change the dependence structure of the de Bruijn graph. De Bruijn graphs are used in a number of applications including genome sequencing [Compeau et al., 2017], and the mathematics of juggling [Ayyer et al., 2015].
Let X = {X 1 , X 2 , . . . , X n } be a length n sequence of binary random variables such that X i is contained within the set of s = 2 letters, V = {0, 1}. The sequence X can be written in terms of words of length m specified by the de Bruijn structure above; X = {W 1 , W 2 , . . . , W n−m+1 }, where each W i is an overlapping sequence of m letters. X is defined to be a sequence of order m. We thus formally define a de Bruijn process in Definition 2.1.
Definition 2.1 (de Bruijn Process (DBP)). Let random variable words W 1 , W 2 , . . . consist of length m sequences of 'letters' from the set V = {0, 1}. For any positive integer, t, and possible states (words), i 1 , i 2 , . . . , i t , X is described by a stochastic process in the form of a Markov chain with, is the probability of transitioning from the word i t−1 to word i t . For the remainder of this paper p j i is used to denote the probability of transitioning from the word (state) i to the word (state) j.
Since X has a Markov property [Billingsley, 1961, Ching and Ng, 2006, Isaacson and Madsen, 1976 on the word and not the letter, far more structure can be incorporated into the sequence. The length of each word (or order) defines how spread the correlation is over nearby letters (the distance over which letters are correlated). For a length m = 2 de Bruijn graph, each letter is dependent on two letters back and, for a length m = 3 de Bruijn graph, each letter is dependent on three letters back. This structure can be observed in the graph in Figure 1. Each word can only transition to other words in which the initial m − 1 letters are equivalent to the last m − 1 letters of the previous word. For example, if m = 3, then the word 001 can only transition to either the word 010 or 011. If m = 1, then the model collapses down to be classically Markov (Markov on the letter), and if m = 0, this is equivalent to independent Bernoulli trials.
The main aim for developing this de Bruijn process is to be able to model sequences of 0's and 1's so that the amount of correlation between variables can be controlled. This is possible by altering the word length and associated transition probabilities, to control the amount of clustering between like letters. For example, let the probabilities be: {p 01 00 = 0.1, p 11 01 = 0.9, p 01 10 = 0.1, p 11 11 = 0.9} in a length m = 2 DBP (see top graph in Figure 1), this will ensure that there is a high level of correlation forcing the letters to cluster together for both 0's and 1's, and avoiding changes between values. Equivalently, the model can be made very anti-clustered by choosing the transition probabilities that retain the current letter to be small. Examples are given in Section 2.2 to emphasize the full range of flexibility.

Stationary Distribution
Let X = {X 1 , . . . , X n } be a sequence of binary random variables such that X can be specified in terms of its de Bruijn words of length m, X = {W 1 , . . . , W n−m+1 }. The de Bruijn structure for letters implies a Markov chain for words. Hence, the stationary distribution [Isaacson andMadsen, 1976, Jones andSmith, 2001] of a word or state j ∈ V m from an irreducible, persistent, aperiodic Markov chain [Norris, 1997, Ross, 2014 with word length m de Bruijn structure is defined as: where p j i (t) = P [W k+t = j|W k = i] are the transition probabilities at the t th time step. The stationary probabilities, π m , must also satisfy the following: k for k ∈ V m and by letting t → ∞. The stationary distribution is hence found when π m(t) remains unchanged as t increases and gives the marginal or long run time probabilities of the de Bruijn words. Given a long enough sequence X with sufficient burn-in period, π m states the proportion of each word occurring in the sequence.
The stationary distribution for words of length m ≥ 2 are given in Theorem 2.2. To simplify the notation, and to apply to a general word length, the words have been written in terms of the decimal representation of their binary values. This form is repeatedly used for the remainder of this paper and is formally written as, m i=1 k i 2 i−1 , where k i ∈ {0, 1} is each letter in the word; e.g. for the word 010, m i=1 k i 2 i−1 = k 1 2 0 + k 2 2 1 + k 3 2 2 = 0 · 1 + 1 · 2 + 0 · 4 = 2.
Theorem 2.2 (Stationary Distribution m ≥ 2). For a general case m, the system of equations for the marginal probabilities, π m , can be expressed as follows: for i = 0 : (2 m − 1). Solving these equations gives the following relationships: There are also relationships between the stationary distribution and the transition probabilities which we can express using the law of total probability [Feller, 1950]. The law of total probability states that π 1 (j) = n P (j|i n )π m (i n ), so that the probability of a letter (0 or 1) is the weighted average of all the possible words that could generate that letter. In other words, we average over all of the possible starting words to get the probability of obtaining a single 0 or 1.

Examples
We now present several examples of how the transition probabilities (p j i ) effect the distribution of letters in sequences produced from de Bruijn processes. It is assumed that binary sequences are unconditional on the starting letters. We hence simulate for a sufficient amount of time to be in steady state, taking a large enough burn-in to ensure that this does not affect any results.
Each panel in Figure 2 gives one of four examples of running a word length m = 2 DBP for n = 200 time steps. The marginal probabilities are kept the same (π 1 (0) = π 1 (1) = 0.5) so that we can make a better comparison of the spread of letters. We would expect that like letters cluster in larger blocks when the probabilities for remaining at the same letter are kept close to one. The transition probabilities ({p 01 00 , p 11 01 , p 01 10 , p 11 11 })) for the four examples (from top to bottom) take the values {0.9, 0.1, 0.9, 0.1}, {0.5, 0.5, 0.5, 0.5}, {0.25, 0.75, 0.25, 0.75} and {0.1, 0.9, 0.1, 0.9} respectively. The light blue areas represent simulated 0's and the dark blue areas represent simulated 1's.
As the sequences progress from top to bottom, the letters are becoming far more clustered together in blocks. Each sequence has equal numbers of 0's and 1's (50% of each), but vary in their distribution across the sequence. The top sequence is designed to be constantly swapping between letters (anti-clustered), the second sequence is equivalent to independent Bernoulli trials, the third sequence is slightly more clustered to like letters, and the final sequence is very clustered into blocks of 0's and 1's. Figure 3 shows that it is not necessary to keep the marginal probabilities of the letters equal to 0.5. In this plot we have two sequences from m = 2 de Bruijn processes where the top sequence has equivalent marginal probabilities, but the bottom sequence is set to have π 1 (0) = 0.2 and π 1 (1) = 0.8. The corresponding transition probabilities for these are {0.1, 0.8, 0.2, 0.9} and {0.775, 0.8, 0.8, 0.9} respectively. The transition probabilities in this example are defined such that we have p 11 01 = 0.8 and p 11 11 = 0.9 for both, but the frequency of 1's differs. Hence the lengths of continuous 1's in a row is equivalent, but they occur far more often in the latter sequence.  The word length of the de Bruijn process also has an effect on the distribution of the letters across the sequence, but it is not as easy to make comparisons due to the different numbers of transition probabilities required for each word. Selecting a longer word length means that more structure can be incorporated through the more intricate transition probabilities, and so it is a mixture of both the word lengths and the transition probabilities that causes the largest impact on the clustering of the sequence.

Run Length Distribution
A run length, R, is defined as the number of consecutive 1's (or 0's) in a row bounded by a 0 (or 1) at both ends. To quantify how clustered a sequence generated from a de Bruijn process is, we can observe the run lengths of letters. Without loss of generality, we will only quantify run lengths of 1's, but all of the following results can be extended to be in terms of run lengths of 0's. We also note that all results are conditional on a run existing (i.e. at least a run of R = 1).
The distribution of run lengths R for a word length m = 2 de Bruijn process with transition probabilities {p 10 01 , p 11 01 , p 10 11 , p 11 11 } is given in Lemma 3.1. This distribution gives the probability of a run of 1's of length R = n, for any n ∈ N + . The probability of a single 1 is simply the transition probability p 10 01 , since this gives the probability of transitioning from the word 01 to the word 10 (giving the sequence 010). For higher run lengths we again start with the word 01, but now transition to the word 11, with probability p 11 01 . For a run of length R = 2, we finish with the probability p 10 11 , whilst for a run length of R ≥ 3, we would further transition to the same word 11 with transition probability p 11 11 until a run length of R = n is reached. The transition probability is raised to the power n − 2 since a run length of two is already created with the other two transition probabilities in the equation. Every run must finish with the probability p 10 11 since a zero signifies the end of the run.
Lemma 3.1 (Run Length Distribution, m = 2). Let the random variable R denote the run length of a sequence of 1's from a sequence X. The probability density function of R in terms of length m = 2 de Bruijn words is given as: for n = 1 p 11 01 (p 11 11 ) n−2 p 10 11 for n ≥ 2,

Proof. See Appendix
In Theorem 3.2 the run length distribution is extended for de Bruijn words of length m ≥ 3, where the decimal representation of the binary notation is used. The general form of the distribution is the same as that seen in Lemma 3.1. We retain the term that represents consecutively drawing 1's p 2 m −1 2 m −1 , however there is now a longer burn-in period until this point is reached. We refer the reader to Figure 4. Each line of this figure shows how a run of length n 1's evolves through subsequent transitions. To start a sequence of 1's off, the first word must take the form * 01, where * represents any possible sequence of length m − 2. Since there are 2 m−2 possibilities for the letters represented by * , we must average over all possibilities using the law of total probability [Feller, 1950] to get the full run length distribution as follows: where π m−2 ( * i ) is the marginal probability for the ith initial sequence of letters. Since, * is of length m − 2, we must then represent each π m−2 ( * i ) in terms of the length m de Bruijn words.We apply the law of total probability again, such that: For longer word lengths m, there are not only more possible starting sequences (given by * ), but the burn-in periods are also longer (as shown in Figure 4). If the run does not reach the state containing all 1's before a 0 is drawn, we finish by transitioning to a word of the form * 10. Alternatively, all runs of length R = n ≥ m will end up reaching the all 1 state with probability p 2 m −1 2 m −1 . The run then ends with the first 0, resulting in a word of the form 11..10 with transition probability p 2 m −2 2 m −1 . The state with all 1's must also be raised to the power n − m since a length of m 1's is obtained in the run up period.
Theorem 3.2 (Run Length Distribution, m ≥ 3). Let the random variable R denote the run length of a sequence of 1's from a sequence X. The probability density function of R in terms of length m de Bruijn words is given as: Proof. See Appendix Figure 4: Diagram representing the burn-in period for sequences of 0's and 1's which generate a run of length n 1's. Letters represented by * can take either a 0 or a 1 for the run to be valid.
The discrete geometric distribution, given by P (X = k) = (1 − p) k p for k = 0, 1, 2, ... [Johnson et al., 2005, Feller, 1950, gives the probability that the first success (1) from a set of independent Bernoulli trials, X, with probability of success, p, happens on the (k + 1) th trial. We can consider the run length distribution as a generalised geometric distribution with the additional burn-in period of a run. Although the geometric distribution considers independent trials, we still retain correlation through the Markov property on the words, which ensures the transitions have a unique ordering to match the words to the sequence. We will be using this link to the geometric distribution in following calculations.

Expectation and Variance of Run Length
From the run length distribution, we can now calculate subsequent measures including the expectation [Feller, 1950]. The expected run length for a length m = 2 de Bruijn process is given in Lemma 3.3 and the expected run length for a de Bruijn process with word length m ≥ 3 is given in Theorem 3.4. The expectation of a geometric random variable is utilised in both expressions.

Proof. See Appendix
Theorem 3.4 (Expected Run Length, m ≥ 3). Given the run length distribution in Theorem 3.2 for random variable R, the expected run length is given as: p 2 j+2 (i mod 2 m−j−2 )+(2 j+1 −1)−1 j=n 2 j+1 (i mod 2 m−j−1 )+(2 j −1) Proof. See Appendix By using the well known result, ) 2 (as well as results from the geometric distribution), we can calculate the variance for the run length distribution. The squared expectation of the run length distribution for the m = 2 word length case is given in Lemma 3.5 and the squared expectation of the run length distribution for word lengths m ≥ 3 is shown in Theorem 3.6. Proof. See Appendix Theorem 3.6 (Squared Expectation of Run Length, m ≥ 3). Given the run length distribution in Theorem 3.2 for random variable R, the squared expectation of the run length is given as: where, π m (j)

Generating Functions of Run Length
Generating functions [Wilf, 1994, Johnson et al., 2005] are a powerful tool which are often used as an alternative to a distribution when it is too complex to work with. The following generating functions are formed from recurrence relationships which give the general form for the sequence of run lengths. Again, since there is a relationship between the de Bruijn run length distribution and a geometric distribution, both the probability and moment generating functions for the geometric distribution can be utilised.
The probability generating function for the run length distribution takes the form G R (t) = E[t n ] = ∞ n=1 t n P (R = n). The result for m = 2 is given in Lemma 3.7 and the result for m ≥ 3 is given in Theorem 3.8. The latter has the same structure as both the geometric and m = 2 cases, with the slight difference that there is now an additional polynomial which comes from the initial burn-in period before the all-1 word is reached. Once we have reached this word, and continue drawing 1's to the chain, then this becomes equivalent to the geometric distribution.
Lemma 3.7 (Run Length Probability Generating Function, m = 2). Given the m = 2 run length distribution in Lemma 3.1 for random variable R, the probability generating function of the run length is given as follows: G R (t) = (p 11 01 p 10 11 − p 10 01 p 11 11 ) t 2 + p 10 01 t 1 − p 11 11 t Proof. See Appendix Theorem 3.8 (Run Length Probability Generating Function, m ≥ 3). Given the run length distribution in Theorem 3.2 for random variable R, the probability generating function of the run length is given as follows: . . .

Proof. See Appendix
We can also calculate the moment generating function of the run length distribution, which take the form M R (t) = E[e nt ] = ∞ n=1 e nt P (R = n) for n = 0, 1, 2, .... The moment generating functions when m = 2 and m ≥ 3 are given in Lemma 3.9 and Theorem 3.10 respectively. As for the case in Theorem 3.8, we can see the similarities to the geometric case, but with the addition of a polynomial term representing the burn-in period to the de Bruijn word consisting of all 1's. Lemma 3.9 (Run Length Moment Generating Function, m = 2). Given the m = 2 run length distribution in Lemma 3.1 for random variable R, the moment generating function of the run length is given as follows: M R (t) = (p 11 01 p 10 11 − p 10 01 p 11 11 ) e 2t + p 10 01 e t 1 − p 11 11 e t

Proof. See Appendix
Theorem 3.10 (Run Length Moment Generating Function, m ≥ 3). Given the run length distribution in Theorem 3.2 for random variable R, the moment generating function of the run length is given as follows:

Examples
Consider the four m = 2 examples in Figure 2 from Section 2.2. We can compare the theoretical results obtained from the de Bruijn run length against the actual run lengths generated from these simulations. Table 1 shows the probabilities of getting a run of length n for n = 1, . . . , 10 from the run length distribution in Lemma 3.1. They are presented in the same order as above so that DBP 1 refers to the de Bruijn process with transition probabilities {0.9, 0.1, 0.9, 0.1} and DBP 4 refers to probabilities {0.1, 0.9, 0.1, 0.9}. The distributions of run lengths from each sequence are also depicted in the histograms in Figure 5 along with the theoretical run length distributions as seen in Table 1. DBP 1 has the highest chance of a short run length, where the probability for a run length R ≥ 2 quickly becomes very small. There is a 90% chance of a run length of R = 1 since this process is designed to constantly be alternating between the two letters. DBP 2 has a 50% chance of a run of length R = 1. The chance then continues to drop for higher processes until DBP 4 only has a 10% chance of a run of length R = 1. We also notice that all of the given probabilities for DBG 3 and DBG 4 stay fairly similar for all run lengths as compared to DBP 1 and DBP 2 which tend to reduce at a fast rate.  Table 2. Here we have given both the analytical results generated from Theorems 3.4 and 3.6 and the simulated results from the length n = 200 sequences in Figure 2. For all cases, both values are shown to agree with on average a difference of 2.18%. Hence, this gives confidence that both the simulated and analytic run lengths are calculated correctly.
The expected run length of the structured example (DBP 1) is 1.11, which is to be expected since the de Bruijn process is designed to be constantly alternating letters. The variance for this process is 0.12, which confirms that the process is designed to be very structured and does not alter much from having an average run length close to R = 1. These both match exactly to the theoretical results. The random Bernoulli sequence has both an expected run length of 2 and variance of 2. For DBP 3 and DBP 4, since the de Bruijn processes are designed to generate more clustered sequences, the expected run lengths are also shown to increase. This is to be expected since the transition probabilities have a high probability of remaining at the same letter. The variance for these de Bruijn processes also increase significantly; being more likely to generate several very long sequences as well as a few short ones by chance. We observe this effect in the histograms in Figure 5, where the spread of run lengths gets much larger for DBP 3 and DBP 4. Occasionally we observe a very long run length, but we also occasionally draw very short lengths. There is also likely to be more variance in the simulations, especially if the de Bruijn Markov chains are not run for a sufficient amount of time. Table 2 also shows values for two standard deviations of the sample expected run lengths and two standard deviations of the sample variance of run lengths. These are calculated using the expressions, σ 2 n and µ 4 n − σ 4 (n−3) n(n−1) respectively, where n = 200 and µ 4 is the kurtosis. This is the fourth central moment and is found from the fourth differential of the cumulant generating function of the run length distribution evaluated at zero. All the differences between the analytical and simulated expectations and variances are within ± two standard deviations.

Application: Precipitation
As an example application we finally investigate whether the amount of precipitation observed can be successfully described using a de Bruijn process structure. In order to do this, we have obtained daily weather recordings from a human-facilitated observation station used in the Global Historical Climatology Network-Daily (GHCN) database. There are 455 daily measurements of precipitation recorded at a station in Eskdalemuir, UK from January 2021 to March 2022. The data gives the total amount of precipitation (rain, melted snow, etc.) recorded in inches in a 24 hour period ending at the observation time.
We further translate this to binary data, recording if either precipitation was observed (1) or no precipitation was observed (0). The full data is given in Figure 6. Since the amount of precipitation can depend on the current season, the data is first split up into corresponding seasons (spring, summer, autumn and winter). The corresponding distributions of run lengths are presented in the histograms in Figure 7 in the above order of seasons. There are shown to be different distributions of run lengths for each season such that for the colder seasons, there are both many single days of rain as well as long periods of rain lasting up to 27 days.
A de Bruijn process is then fitted separately to each season, estimating both the word length and associated transition probabilities. See the companion paper (PAPER 2) for more information on model fitting and inference. Both spring and summer were estimated to have a word length of m = 2, whilst autumn and winter were instead estimated to have a word length of m = 3. The black dots seen in Figure 7 correspond to the theoretical run lengths generated by applying the estimated transition probabilities to Lemma 3.1 and Theorem 3.2. In general, the theoretical run lengths match well with the run lengths from the data implying that the data can indeed be represented by a de Bruijn structure.
We have further compared the theoretical run lengths from the de Bruijn processes with that of a geometric distribution. The black crosses in Figure 7 give the probabilities of a run of length n generated from a geometric distribution for each given season. It is evident that in each season, the run lengths produced match those from the de Bruijn process far better than those produced from the geometric distribution. We can conclude that the precipitation recorded has correlation structures that are far more complex than that of Bernoulli trials or a Markov property.

Discussion
We have introduced a novel method for modelling correlated binary random variables, where the strength of the correlation is related to distance. It is often the case that variables close in distance are more likely to be highly correlated then those further away. Hence, by altering the correlation we are able to generate both clustered and anti-clustered sequences of 0's and 1's.
The process uses structures from de Bruijn graphs. These are directed graphs where the nodes of the graph are length m sequences of 0's and 1's. If a Markov property is placed on these length m sequences rather than the variables themselves, then we are able to control the amount of correlation included in the neighbourhood. Examples are presented throughout as well as details of the stationary distribution for the length m sequences.
We have also presented a run length distribution in order to determine how clustered a sequence is likely to be from a given de Bruijn process. The run length distribution gives the probability of a run of 1's of length n bounded by a 0 at each end. From this, we are then able to calculate the expected run length, variance of run length and generating functions.
There are many areas for future work on the topic of de Bruijn graphs and correlated Bernoulli processes. Firstly, we need to establish a method of inference. Given a sequence of binary variables known to be generated from a de Bruijn process, it is important that we have a method to estimate both the de Bruijn word length (length m sequence of each node) and the associated transition probabilities. We will then be able to test the proposed de Bruijn process on real-life applications. This is presented in a companion paper (PAPER 2).
Secondly, we must also consider de Bruijn processes in two (or more) dimensions, as well as non-directional de Bruijn processes. Many applications (particularly in ecology) are defined on a 2-dimensional spatial grid. However, when considering the form of a de Bruijn word in 2d, it becomes very difficult to work in a specific direction (something that is very important for de Bruijn graphs). Hence, we believe further research should focus on trying to remove this natural direction in simulation.
De Bruijn graphs can quickly become very complicated with many possible transition probabilities for large word lengths. Therefore, we believe that it would be useful to try and limit the number of transition probabilities, which could be particularly important for inference. We also intend to develop non-stationary de Bruijn processes where we would be able to create chains that are clustered in some places and anti-clustered in others.

Appendix Theorem 2.1 (Stationary Distribution m ≥ 2)
For a general case m, the system of equations for the marginal probabilities, π m , can be expressed as follows: for i = 0, . . . , (2 m − 1). Solving these equations gives the following relationships: Proof. Expanding out the equation π m T = π m for marginal word probabilities π m , and transition matrix T , gives a total of 2 m equations. Due to the form of the transition matrix, the i th expression consists of the sum of two transition probabilities of the form: where j = 0, 0, 1, 1, 2, 2, . . .. Both transition probabilities have the same starting word but their end words differ by exactly 2 m−1 . This is satisfied by letting j = 1 2 (i − i mod 2), giving the following general equation: If i = 0 this gives the following, Alternatively when i = 1, (1) When i = 2 m − 1, Finally, given i = 2 m − 2, Since there are two terms in each row of the transition matrix, there will be pairs of equations with the same marginal probabilities on the left hand side. For the i th and (i + 1) th equations (counting from zero), we have the following: p i 1 2 (i−i mod 2) π m 1 2 (i − i mod 2) + p i 1 2 (i−i mod 2)+2 m−1 π m 1 2 (i − i mod 2) + 2 m−1 = π m (i) and p i+1 1 2 (i−i mod 2) π m 1 2 (i − i mod 2) + p i+1 j+2 m−1 π m 1 2 (i − i mod 2) + 2 m−1 = π m (i + 1).

Lemma 3.1 (Run Length Distribution, m = 2)
Let the random variable R denote the run length of a sequence of 1's from a sequence X. The probability density function of R in terms of length m = 2 de Bruijn words is given as: for n = 1 p 11 01 (p 11 11 ) n−2 p 10 11 for n ≥ 2, Proof. P (R = n) = p 10 01 for n = 1. For n ≥ 2, the run starts with the word 01 and transitions to 11 with probability p 11 01 creating the sequence 011. The run continues by transitioning to the word, 11, with probability p 11 11 , n − 2 times to create a run of length n. The run ends with the first zero with probability p 10 11 , hence: P (R = n) = p 11 01 (p 11 11 ) n−2 p 10 11 , for n ≥ 2.

Theorem 3.2 (Run Length Distribution, m ≥ 3)
Let the random variable R denote the run length of a sequence of 1's from a sequence X.
The probability density function of R in terms of length m de Bruijn words is given as: for n ≥ m.
Proof. Given a sequence X containing the letters V = {0, 1} of length N , the sequence can be written in terms of its de Bruijn words of length m. Each of these words can be expressed by the numerical representation of their binary form. Transitioning from the word i to the word j is denoted by p j i . For R = n, a run starts with a word of the form * 01 irrespective of the word length m. Since, * denotes a sequence of length m − 2, this gives a total of 2 m−2 words of the form * 01. The law of total probability states: where π m−2 ( * i ) is the marginal probability for the i th starting sequence of length m − 2. Runs of length n = 1 are generated by transitioning from a word of the form * 01 to a word of the form * 010, where * is the same sequence as denoted by * but with the first symbol removed. Hence: P (R = 1| * ) = p * 010 * 01 . The m−2 sequences given by * are the binary representations of the numbers 0, 1, . . . , 2 m−2 − 1. Since the sequence * starts in the third position (from right to left) representing 2 2 , all starting words of the form * 01 are given by 4i + 1 ∀ i ∈ (0 : 2 m−2 − 1).
The m−3 sequences given by * are the binary representations of the numbers 0, 1, . . . , 2 m−3 − 1. When * = 0, the word is of the form 00 . . . 010 and translates to 2 in the numerical representation. For values of * higher than 0, the full word increases by 2 3 each time since * starts in the fourth position (from right to left). Hence all words of the form * 010 are given by (2 3 )i + 2 ∀ i ∈ (0 : 2 m−3 − 1).
Since there are only 2 m−3 words to transition to, half of the initial 2 m−2 words must transition to the same following word. This forms a repeated cycle of recurring end words that correspond to the starting words in logical order (1, . . . , m). The cycle repeats every 2 m−3 words, and so depends on the word length m. Combing this with the law of total probability gives: For P (R = n) where n = 2, . . . , m−1, a run must start with a word of the form * 01 and transition to a word of the form * 011 with probability p * 011 * 01 . The n − 1 transitions that follow each reduce * by one letter but include an additional 1 to the end of the sequence.
Each transition causes the sequence * to reduce by one each time and the remaining letters are shifted to the left. This then causes the numerical binary number to increase by a power of 2 for each transition. Since * must start in the third position (from right to left), we have: 2 j+1 ∀ j ∈ (1 : n). This is the increase in value for each i. As for the n = 1 case, a repeated cycle of recurring end words occurs. To begin with, the cycle repeats every 2 m−2 , however as the sequence * decreases with increasing transitions, this decreases the total number of possible words and the length of the cycle to give: 2 j+1 (i mod 2 m−j−1 ) ∀ i ∈ (1 : 2 m−2 − 1), j ∈ (0 : n − 1).
The numerical representation of a sequence of n consecutive 1's is given by 2 n − 1. Again, as the sequence moves to the left with each transition, the numerical binary number increases by a power of 2 for each transition. This gives: 2 j − 1 ∀ j ∈ (1 : n).
Finally, the run of 1's is stopped with a word of the form * 10, where * is a length m − 2 sequence containing n − 1 consecutive 1's. An indicator is used to remove one from the numerical representation of the last transition when a 0 occurs.
Since the end word from the previous transition is the starting word of the next this gives: for n = 2, ..., m − 1.
For run lengths n ≥ m, the initial m − 1 transitions will be identical to the n < m case with the indicator removed. The remaining n − m transition probabilities are of the form p 11...1 11...1 = p 2 m−1 2 m−1 with the last word of the form 11...10 = 2 m − 2, hence: , for n ≥ m. For π m−2 (i) ∀ i ∈ (0 : 2 m−2 − 1), the law of total probability states: where the numerical representation of the binary words are used. Since the run must start with a word of the form * 01, π m−2 (i) represents the marginal probability of all possible sequences given by * of length m − 2. π m is defined in terms of length m words and associated transition probabilities. Finding the conditional probabilities involves m − 2 transition probabilities for each j ∈ (0 : 2 m − 1) since there are 2 m words in an m length de Bruijn graph. First consider the case when i = 0 for any word length m (in the numerical representation of the binary word, 0...0). The first word for each j, is always j itself. For each following transition, for each j, the letters shift to the left, and a 0 is added to the end of the word. Hence, the word becomes an extra factor of 2 larger at each transition i.e. 2 k × j for each k ∈ (0 : m − 3). However, since the word transitions to a word of the form * 0...0 each time, where * reduces in length by one, some words can transition to the same word, creating a cycle on j. Initially all 2 m words are unique, but the cycle reduces with each transition (as * becomes smaller), so the words become: 2 k (j mod 2 m−k ).
If i ∈ (0 : 2 m−2 − 1), the initial pattern is the same as for i = 0, with different contributions to each word for each transition. 1's must be added to the base case (i = 0) when they turn up in the transitions for any i. For all 2 m−2 possible binary sequences of length m − 2, the 1's appear in cycles depending on which column they lie in. From the left entry, the first 2 m−3 words begin with a 0 and the last 2 m−3 words begin with a 1. The next column repeats 0's and 1's every 2 m−4 words. This continues until the last column alternates between 0's and 1's for each i. This reduces to ( 1 2 m−s−2 (i − (i mod 2 m−s−2 )))mod 2, where s ∈ (1 : k) for each transition. Since the current 1's increase by a power of two for each transition, a summation is formed to make sure each is included. Hence the words become: k s=1 2 k−s [( 1 2 m−s−2 (i − (i mod 2 m−s−2 )))mod 2]. Since the end word of the current transition is the start word of the next transition, this gives the following: π m (j)

Lemma 3.3 (Expected Run Length, m = 2)
Given the m = 2 run length distribution in Lemma 3.1 for random variable R, the expected run length is given as: since E[X] = x xP (X = x) = x x(1 − p) x p = 1−p p for a geometric distribution with P (X = x) = (1 − p) x p and x = 0, 1, 2, 3, ....

Theorem 3.4 (Expected Run Length, m ≥ 3)
Given the run length distribution in Theorem 3.2 for random variable R, the expected run length is given as: Proof. π m (j) since E[X] = x xP (X = x) = x x(1 − p) x p = 1−p p for a geometric distribution with P (X = x) = (1 − p) x p and x = 0, 1, 2, 3, ....

Lemma 3.5 (Squared Expectation of Run length, m = 2)
Given the m = 2 run length distribution in Lemma 3.1 for random variable R, the squared expectation of the run length is given as: for a geometric distribution with P (X = x) = (1 − p) x p and x = 0, 1, 2, 3, .... since E[X 2 ] = x x 2 P (X = x) = x x 2 (1 − p) x p = (1−p)(2−p) Theorem 3.8 (Run Length Probability Generating Function, m ≥ 3) Given the run length distribution in Theorem 3.2 for random variable R, the probability generating function of the run length is given as follows: Given the m = 2 run length distribution in Lemma 3.1 for random variable R, the moment generating function of the run length is given as follows: G(x) = (p 11 01 p 10 11 − p 10 01 p 11 11 ) e 2x + p 10 01 e x 1 − p 11 11 e x Proof.
where, a 1 = 2 m−2 −1 i=0 π m−2 (i)p  Figure 2. Black dots give the theoretical run length distribution derived from Lemma (3.1). The run lengths increase as the de Bruijn process gets increasingly clustered. Figure 6: Daily precipitation data recorded at a station in Eskdalemuir, UK. Dark blue (or 1) represents that there was precipitation recorded in the 24 hour period, whilst light blue (or 0) corresponds to no precipitation recorded that day. Figure 7: Histograms showing the seasonal distributions of run lengths of recorded precipitation in a 24 hour period at a station in Eskdalemuir. From top to bottom, the data corresponds to seasons spring, summer, autumn and winter from January 2021 to March 2022. Black dots give the theoretical run length distributions derived from estimated de Bruijn process. Black crosses give a geometric density fitted to the data.