62
$\begingroup$

This is an interesting problem I came across. I'm attempting to write a Python program to get a solution to it; however, I'm not sure how to proceed. So far, I know that I would expect the counts of heads to follow a binomial, and length of runs (of tails, heads, or both) to follow a geometric.

Below are two sequences of 300 “coin flips” (H for heads, T for tails). One of these is a true sequence of 300 independent flips of a fair coin. The other was generated by a person typing out H’s and T’s and trying to seem random. Which sequence is truly composed of coin flips?

Sequence 1:

TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHH TTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHH TTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHT THHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHT HTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTT HHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT

Sequence 2:

HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTH THTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHH TTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTT THTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTH HHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHH HTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT

Both sequences have 148 heads, two less than the expected number for a 0.5 probability of heads.

$\endgroup$
17
  • 3
    $\begingroup$ Naive qustion : Does the question means: can I find with that porbablity that such sequence is generated by human is big enough ? Can I compute such probablity ? Every sequence is possible, so this is not a proof but only possibility. Am I right ? $\endgroup$
    – Adam
    May 7 at 8:53
  • 20
    $\begingroup$ (I find it truly amusing that somebody voted to close this thread because it asks for opinions. If that were true, we should close all threads here on CV that use statistical analysis to compare data--and that would leave us with nothing but lists of references!) $\endgroup$
    – whuber
    May 7 at 13:42
  • 3
    $\begingroup$ Compress both strings and see which one comes out shorter. $\endgroup$ May 7 at 13:58
  • 5
    $\begingroup$ @whuber That comment is probably a reference to (approximate) Kolmogorov complexity. See algorithmically random sequence. $\endgroup$
    – user76284
    May 7 at 15:45
  • 3
    $\begingroup$ @whuber With very high probability, a truly random sequence cannot be compressed. Compression methods can leverage many kinds of bias, not only if one symbol appears more overall, but also if some subsequences appear more than others (as here, H is more likely to be followed by T and vice versa, per your answer). $\endgroup$
    – nanoman
    May 7 at 18:20

11 Answers 11

85
$\begingroup$

This is a variant on a standard intro stats demonstration: for homework after the first class I have assigned my students the exercise of flipping a coin 100 times and recording the results, broadly hinting that they don't really have to flip a coin and assuring them it won't be graded. Most will eschew the physical process and just write down 100 H's and T's willy-nilly. After the results are handed in at the beginning of the next class, at a glance I can reliably identify the ones who cheated. Usually there are no runs of heads or tails longer than about 4 or 5, even though in just 100 flips we ought to see a longer run that that.

This case is subtler, but one particular analysis stands out as convincing: tabulate the successive ordered pairs of results. In a series of independent flips, each of the four possible pairs HH, HT, TH, and TT should occur equally often--which would be $(300-1)/4 = 74.75$ times each, on average.

Here are the tabulations for the two series of flips:

   Series 1    Series 2
      H   T        H  T
  H  46 102       71 76
  T 102  49       77 75

The first is obviously far from what we might expect. In that series, an H is more than twice as likely ($102:46$) to be followed by a T than by another H; and a T, in turn, is more than twice as likely ($102:49$) to be followed by an H. In the second series, those likelihoods are nearly $1:1,$ consistent with independent flips.

A chi-squared test works well here, because all the expected counts are far greater than the threshold of 5 often quoted as a minimum. The chi-squared statistics are 38.3 and 0.085, respectively, corresponding to p-values of less than one in a billion and 77%, respectively. In other words, a table of pairs as imbalanced as the second one is to be expected (due to the randomness), but a table as imbalanced as the first happens less than one in every billion such experiments.

(NB: It has been pointed out in comments that the chi-squared test might not be applicable because these transitions are not independent: e.g., an HT can be followed only by a TT or TH. This is a legitimate concern. However, this form of dependence is extremely weak and has little appreciable effect on the null distribution of the chi-squared statistic for sequences as long as $300.$ In fact, the chi-squared distribution is a great approximation to the null sampling distribution even for sequences as short as $21,$ where the counts of the $21-1=20$ transitions that occur are expected to be $20/4=5$ of each type.)


If you know nothing about chi-squared tests, or even if you do but don't want to program the chi-square quantile function to compute a p-value, you can achieve a similar result. First develop a way to quantify the degree of imbalance in a $2\times 2$ table like this. (There are many ways, but all the reasonable ones are equivalent.) Then generate, say, a few hundred such tables randomly (by flipping coins--in the computer, of course!). Compare the imbalances of these two tables to the range of imbalances generated randomly. You will find the first sequence is far outside the range while the second is squarely within it.

Figure

This figure summarizes such a simulation using the chi-squared statistic as the measure of imbalance. Both panels show the same results: one on the original scale and the other on a log scale. The two dashed vertical lines in each panel show the chi-squared statistics for Series 1 (right) and Series 2 (left). The red curve is the $\chi^2(1)$ density. It fits the simulations extremely well at the right (higher values). The discrepancies for low values occur because this statistic has a discrete distribution which cannot be well approximated by any continuous distribution where it takes on small values -- but for our purposes that makes no difference at all.

$\endgroup$
11
  • 4
    $\begingroup$ Why do you say, "This case is subtler" (than just looking for runs)? OP's sequence 1 has length 300 (longer than the 100 you mention) and its longest run is a single HHHH. So it doesn't seem subtle at all, but rather a slam dunk for the "at a glance" method you describe. $\endgroup$
    – nanoman
    May 7 at 18:11
  • 9
    $\begingroup$ It's not proper to apply the chi squared test to the contingency table given since consecutive pairs are overlapping and not independent, so p-values will be exaggerated. I suspect the result will be similar after addressing this though. $\endgroup$
    – Paul
    May 7 at 18:11
  • 5
    $\begingroup$ @Paul That's an excellent point, which I have neglected to discuss (and it's no excuse to claim, as I was hoping to do, that it's implicitly handled correctly in the last paragraph, because I did not advertise that fact). $\endgroup$
    – whuber
    May 7 at 18:13
  • 3
    $\begingroup$ @klm123 But as long as you don't actually try a bunch of tests before finding an improbable result (p-hacking), the result is still significant. The test(s) should be chosen in advance. There are also ways of correcting p-values if you run multiple tests. Yes, there is the potential to be misleading if you only report the significant test and not the others -- this has contributed to the "replication crisis". $\endgroup$
    – nanoman
    May 7 at 23:07
  • 3
    $\begingroup$ My answer talks about why this test in particular is well-motivated. People tend to think that non-repeated values are more random than repeated values. (Which is true in some sense, but people believe it to a greater extent / in more generality than is actually true.) $\endgroup$
    – Paul
    May 7 at 23:25
32
$\begingroup$

There are two very good answers as of writing this, and so let me add a needlessly complex yet interesting approach to this problem.

I think one way to operationalize the human generated vs truly random question is to ask if the flips are autocorrelated. The hypothesis here being that humans will attempt to appear random by not having too many strings of one outcome, hence switching from heads to tails and tails to heads more often than would be observed in a truly random sequence.

Whuber examines this nicely with a 2x2 table, but because I am a Bayesian and a glutton for punishment let's write a simple model in Stan to estimate the lag-1 autocorrelation of the flips. Speaking of Whuber, he has nicely laid out the data generating process in this post. You can read his answer to understand the data generating process.

Let $\rho$ be the lag 1 autocorrelation of the flips, and let $q$ be the proportion of flips which are heads in the sequence. A fair coin should have 0 autocorrelation, so we are looking for our estimate of $\rho$ to be close to 0. From there, we only need to count the number of occurrences of $H,H$, $H, T$, $T, H$ and $T,T$ in the sequence.

The Stan model is shown below

data{

  int y_1_1; //number of concurrent 1s
  int y_0_1; //number of 0,1 occurrences
  int y_1_0; //number of 1,0 occurrences
  int y_0_0; //number of concurrent 0s


}
parameters{
  real<lower=-1, upper=1> rho;
  real<lower=0, upper=1> q;

}
transformed parameters{

  real<lower=0, upper=1> prob_1_1 = q + rho*(1-q);
  real<lower=0, upper=1> prob_0_1 = (1-q)*(1-rho);
  real<lower=0, upper=1> prob_1_0 = q*(1-rho);
  real<lower=0, upper=1> prob_0_0 = 1 - q + rho*q;
}
model{
  q ~ beta(1, 1);
  target += y_1_1 * bernoulli_lpmf(1| prob_1_1);
  target += y_0_1 * bernoulli_lpmf(1| prob_0_1);
  target += y_1_0 * bernoulli_lpmf(1| prob_1_0);
  target += y_0_0 * bernoulli_lpmf(1| prob_0_0);
  
}

Here, I've placed a uniform prior on the autocorrelation

$$ \rho \sim \mbox{Uniform}(-1, 1) $$

and on the probability of a head

$$ q \sim \operatorname{Beta}(1, 1) $$

Our likelihood is Bernoulli, and I have weighted the likelihood by the number of occurrences of each pair of outcomes. The probabilities of each outcome (e.g. probability of observing a heads conditioned on the previous flip being a heads) is provided by Whuber in his linked answer. Let's run our model and compare posterior distributions for the two sequences

enter image description here

The estimated auto correlation for sequence 1 is -0.36, and the estimated autocorrelation for sequence 2 is -0.02 (close enough to 0). If I was a betting man, I'd put my money on sequence 1 being the sequence generated by a human. The negative autocorrelation means that when we see a heads/tails we are more likely to see a tails/heads! This observation lines up nicely with the 2x2 table provided by Whuber.

Code

The plot I present is made in R, but here is some python code to do the same thing since you asked

import matplotlib.pyplot as plt
import cmdstanpy

# You will need to install cmdstanpy prior to running this code

# Write the stan model as a string.  We will then write it to a file
stan_code = '''
data{

  int y_1_1; //number of concurrent 1s
  int y_0_1; //number of 0,1 occurences
  int y_1_0; //number of 1,0 occurences
  int y_0_0; //number of concurrent 0s
  

}
parameters{
  real<lower=-1, upper=1> rho;
  real<lower=0, upper=1> q;
}
transformed parameters{
  real<lower=0, upper=1> prob_1_1 = q + rho*(1-q);
  real<lower=0, upper=1> prob_0_1 = (1-q)*(1-rho);
  real<lower=0, upper=1> prob_1_0 = q*(1-rho);
  real<lower=0, upper=1> prob_0_0 = 1 - q + rho*q;
}
model{
  q ~ beta(1, 1);
  target += y_1_1 * bernoulli_lpmf(1| prob_1_1);
  target += y_0_1 * bernoulli_lpmf(1| prob_0_1);
  target += y_1_0 * bernoulli_lpmf(1| prob_1_0);
  target += y_0_0 * bernoulli_lpmf(1| prob_0_0);
  
}
'''

# Write the model to a temp file
with open('model_file.stan', 'w') as model_file:
    model_file.write(stan_code)
    

# Compile the model
model = cmdstanpy.CmdStanModel(stan_file='model_file.stan', compile=True)

# Co-occuring counts for heads (1) and tails (0) for each sequence
data_1 = dict(y_1_1 = 46, y_0_0 = 49, y_0_1 = 102, y_1_0 = 102)
data_2 = dict(y_1_1 = 71, y_0_0 = 75, y_0_1 = 76, y_1_0 = 77)

# Fit each model
fit_1 = model.sample(data_1, show_progress=False)
rho_1 = fit_1.stan_variable('rho')

fit_2 = model.sample(data_2, show_progress=False)
rho_2 = fit_2.stan_variable('rho')




# Make a pretty plot
fig, ax = plt.subplots(dpi = 240, figsize = (5, 3))

ax.set_xlim(-1, 1)
ax.hist(rho_1, color = 'blue', alpha = 0.5, edgecolor='k', label='Sequence 1')
ax.hist(rho_2, color = 'red', alpha = 0.5, edgecolor='k', label='Sequence 2')
ax.legend()
```
$\endgroup$
3
  • 4
    $\begingroup$ +1 I developed my tabular explanation by going through the exercise you lay out here. I plotted the PACF functions of random sequences and compared them, visually, to the PACF functions of the two sequences in the question. One plot--the first sequence--stood out: it had a very negative lag-one partial autocorrelation coefficient. But motivating and explaining the PACF seemed like overkill for this problem ;-). The more enduring lesson, illustrated here and in the post by @COOLSerdash, is that by finding a suitable way to visualize data, we can discover otherwise hidden things about them. $\endgroup$
    – whuber
    May 7 at 13:28
  • 2
    $\begingroup$ (+1) Interesting resolution that does not rely on insufficient statistics, contrary to the others!, but I would have gone fully Markov and put a prior distribution on both $(p_{11},p_{10})$ and $(p_{01},p_{00})$ rather than introducing a correlation $\rho$, but it is unlikely the result would differ. More generally, the alternative to being an iid Uniform sequence could be anything, so restricting to an order one Markov is favouring the null. $\endgroup$
    – Xi'an
    May 9 at 8:09
  • $\begingroup$ @Xi'an You're right in that the result is not much different. The reason I prefer the $\rho$, $q$ parameterization is because we actually have good priors on these (despite what my answer uses lol). Humans have an intuitive sense for when they are acting too correlated and when they are acting too biased, so we know $\rho$ is going to be close to 0 and $q$ close to 0.5. But, you could rewrite this model using a multinomial likelihood and place priors on the probabilities directly too. That's what I love about Bayes! $\endgroup$ May 9 at 17:48
26
$\begingroup$

This is a class activity I've first read about in the book Teaching Statistics. A Bag of Tricks, 2nd ed. by Andrew Gelman and Deborah Nolan (they recommend 100 flips, though). Their reasoning to detect the fabricated sequence is based on the combination of the longest run and the number of runs. For the following plot, I simulated 5000 fair coin tosses of length 300 and plotted the longest run on the y-axis and the number of runs on the x-axis (I once asked a question about the explicit joint probability). Each dot represents the result of 300 fair flips. For better visibility, the points are jittered. The numbers for the two sequences are plotted in color. The conclusion is obvious.

Longest_Run

For a quick calculation, recall that a rule of thumb for the longest run of either heads or tails in $n$ tosses is$^{[1]}$ $l = \log_{1/p}(n(1-p)) + 1$. For an approximate 95% prediction interval, just add and subtract $3$ from this value. Surprisingly, this number (i.e. $\pm 3$) does not depend on $n$! Applied to a fair coin with $n=300, p=1/2$, we have $l=\log_2(300/2) + 1=8.22$. So we expect the longest run to be round $8$ and reasonably in the range of $8\pm 3$, so between $5$ and $11$. The longest run in sequence 2 is $8$, whereas it is $4$ in sequence 1. As this is outside the approximate prediction interval, we'd conclude that sequence 1 is suspicious under the assumption of $p=1/2$.

$[1]$ Schilling MF (2012): The Surprising Predictability of Long Runs. Math. Mag. 85: 141-149. (link)

$\endgroup$
2
  • 1
    $\begingroup$ Excellent answer and graph. "Surprisingly, this number does not depend on n!" You're talking about +-3, right? I was confused for a while since you define $l$ just before, which obvisouly depends on $n$. $\endgroup$ May 9 at 14:48
  • $\begingroup$ @EricDuminil Yes, sorry. The value that does not depend on $n$ is the $\pm 3$ for the prediction interval. $\endgroup$ May 9 at 15:45
14
$\begingroup$

The runs test (NIST page) is a nonparametric test designed to identify unusual frequencies of runs. If we observe $n_1$ heads and $n_2$ tails, the expected value and variance of the number of runs are:

$$\mu = {2n_1n_2 \over n_1+n_2} + 1$$ $$\sigma^2 = {2n_1n_2(2n_1n_2 - n_1 - n_2) \over (n_1+n_2)^2(n_1+n_2+1)}$$

As a rule of thumb, for $n_1, n_2 \geq 10$ the distribution of the observed number of runs is reasonably well-approximated by a Normal distribution.

Edit: (incorporating Eric Duminil's work below)

For sequence 1, we have 148 heads, 152 tails, 205 runs, and for sequence 2, we have 148 heads, 152 tails and 154 runs. Plugging these numbers into our formulae above give us $z$-scores of 6.5 for the first sequence and 0.58 for the second sequence - extremely strong evidence that the first sequence is fake.

When people fake sequences like this, they tend to greatly underestimate the probability of longer runs, so they don't create as many long(ish) runs as they should. This in turn tends to increase the number of runs beyond that which would be expected. Consequently, when testing for faked data, we might prefer a one-sided test of the alternative hypothesis that there are "too many" runs vs. the null hypothesis that the number of runs is average - at least if we think the sequence was created by a human being.

$\endgroup$
7
  • 4
    $\begingroup$ Yep, this is a classic application of gambler's fallacy. $\endgroup$
    – AdamO
    May 6 at 21:54
  • $\begingroup$ There are 62 runs in the first two rows, AFAICT. For the whole sequence #1 : 148 heads, 152 tails, 205 runs. Vs 148 heads, 152 tails and 154 runs for sequence #2. $\endgroup$ May 9 at 14:56
  • 1
    $\begingroup$ Thanks, @EricDuminil - I've included your efforts in the body of the answer with a citation. $\endgroup$
    – jbowman
    May 9 at 16:46
  • 1
    $\begingroup$ That's good. But it leaves unanswered why you would formulate a question in terms of longest run and then solve it by looking at an indirect proxy (number of runs). Why not just use the longest run as the test statistic? $\endgroup$
    – whuber
    May 9 at 18:13
  • 2
    $\begingroup$ @whuber - it's not just the longest run that's informative, it's the number of longer runs in general, and I don't have a good idea of a cutoff for run length for testing. Having said that, using the longest run as a test statistic with 100,000 randomly generated strings with 148 heads and 152 tails for calculating an approximate p-value gives a p-value of 0.00004 for the first sequence. Maybe I'll expand on that over (my) lunch. $\endgroup$
    – jbowman
    May 9 at 18:28
12
$\begingroup$

Here's an empirical approach, based on compression as a proxy for algorithmic complexity:

import bz2
import random
import statistics

s1 = "TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHHTTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHHTTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHTTHHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHTHTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTTHHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT"
s2 = "HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTHTHTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHHTTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTTTHTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTHHHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHHHTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT"

def compressed_len(s):
    return len(bz2.compress(s.encode()))

trials = []
for x in range(100000):
    sr = "".join(random.choice("HT") for _ in range(300))
    trials.append(compressed_len(sr))

mean = statistics.mean(trials)
stddev = statistics.stdev(trials)

print("Random trials:")
print("Mean:", mean)
print("Stddev:", stddev)

l1 = compressed_len(s1)
l2 = compressed_len(s2)

o1 = l1 - mean
o2 = l2 - mean

d1 = o1 / stddev
d2 = o2 / stddev

print("Selected trials:")
print("Seq", "Len", "Dev", sep="\t")
print("S1", l1, d1, sep="\t")
print("S2", l2, d2, sep="\t")

Roughly speaking:

  1. Compress a bunch (100k) of random coinflips.
  2. Observe the resulting length distribution. (In this case I'm approximating it as a normal distribution of lengths; a more thorough analysis would check and pick an appropriate distribution instead of blithely assuming normality.)
  3. Compress the input sequences.
  4. Compare with observed distribution.

Result (note: not exactly reproducible due to the use of random trials; if you want to be reproducible add a random seed):

Random trials:
Mean: 105.05893
Stddev: 2.6729774976956002
Selected trials:
Seq Len Dev
S1  88  -6.381995364609942
S2  109 1.4744119632124217

Based on this, I'd say that S1 is the non-random one here. 6.38 standard deviations below the mean is rather improbable.


The nice thing about this approach is that it's relatively generic, and takes advantage of the pre-existing work of a bunch of smart people.

Just be aware of its limitations and quirks:

  1. You want a compression algorithm that's designed for space over compression speed. BZ2 works well enough here.
  2. This doesn't work if the compression algorithm simply gives up and writes a raw block to the output.
  3. A null result does not mean that the sequence is random. It means that this compression algorithm is unable to distinguish this input from random.
$\endgroup$
0
5
$\begingroup$

This is probably an overcomplicated way of looking at it, but for me it's fun, so I present to you...

Moran's I

Now, Moran's I was developed to look at spatial autocorrelation (basically autocorrelation with multiple dimensions), but it can be applied to the 1-dimensional case as well. Some of my interpretations might be a little sketchy, but you can consider your coin flips this way.

To summarize Moran's I, it will consider your neighboring values using a pre-defined matrix. How you define the matrix is up to you, but it can actually be used to consider not just the directly neighboring values, but any values beyond that. Moran's I will produce a value ranging from -1 (perfectly dispersed values) to 1 (perfectly clustered values), with 0 being random.

I wrote up some quick R code. First, setup the data (OP's data and a couple generated data sets to test dispersion and clustering):

seq1 = unlist(strsplit("TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHHTTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHHTTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHTTHHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHTHTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTTHHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT",
                     split = ""))

seq2 = unlist(strsplit("HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTHTHTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHHTTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTTTHTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTHHHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHHHTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT",
                     split = ""))


# Alternate T and H. E.g., THTHTHTHT....
# 'perfectly dispersed'
# Moran's I = -1
seq3 = rep(c("T", "H"), times = 50) 


# 50 of T followed by 50 of H
# 'perfectly clustered'
# Moran's I approaches 1 as the sample size increases to infinity
seq4 = rep(c("T", "H"), each = 50) 


# weights must be a vector with an odd length and the middle value set to 0
# weights are relative and do not have to add to 1
moran <- function(x, weights) {
  x = c(`T` = 0, `H` = 1)[x] # convert T/H to 0/1
  
  N = length(x)

  x_mean = mean(x)
  
  den = sum((x - x_mean)^2)
  
  W = 0  
  num = 0
  
  offset = floor(length(weights)/2)
  x_padded = c(rep(NA, 10), x, rep(NA,10)) # padding for sliding windown
  for (i in 1:length(x)) {
    x_slice = x_padded[(i+10-offset):(i+10+offset)]

    W = W + as.numeric(!is.na(x_slice)) %*% weights
    num = num + (x[i] - x_mean) * sum((x_slice - x_mean) * weights, na.rm = TRUE)
  }

  return(unname((N * num)/(as.numeric(W) * den)))
}

Next, I test the generated data sets to illustrate/test that my function is working correctly:

# Test the 'perfect dispersion' scenario (should be -1)
moran(seq3, c(1, 0, 1))
## [1] -1

# Test the 'perfect clustering' scenario (should be ~1)
moran(seq4, c(1, 0, 1))
## [1] 0.979798

Now, let's look at OP's sequences:

# Simple look at seq1. The weights test the idea that the current flip
# is based purely on the last flip (a reasonable model for how a person might react)
moran(seq1, c(1, 0, 0))
## [1] -0.3647031
moran(seq2, c(1, 0, 0))
## [1] -0.02359453

I'm defining my weights matrix such that only the previous flip is considered when testing for autocorrelation. We see that the second sequence is very close to 0 (random), whereas the first sequence seems to lean somewhat toward overdispersion.

But maybe we think someone faking coin flips would consider the last two flips, not just the most recent:

# Maybe the person is looking back at the last two flips
moran(seq1, c(1, 1, 0, 0, 0))
## [1] -0.1726056
moran(seq2, c(1, 1, 0, 0, 0))
## [1] 0.0249505

The second sequence is just as close to 0 as before, but the first sequence had a pretty noticeable shift towards 0. This might be interpretable in a couple of different ways. First, if we know that the first sequence is fake, then maybe it means the person wasn't considering two flips back. A second interpretation is that maybe they were considering the last two flips, and somehow this led them to doing a better job at faking randomization. A third option might just be sheer dumb luck at faking the randomization.

Now, maybe the person considers the last two coin flips but gives the most recent flip more importance.

# Same idea, but maybe the more recent of the two is twice as important
moran(seq1, c(1, 2, 0, 0, 0))
## [1] -0.2367095
moran(seq2, c(1, 2, 0, 0, 0))
## [1] 0.008750762

Here, we see the two sequences react differently. The second sequence (already pretty close to 0), gets noticeably closer to 0, whereas the first sequence shifts noticeably away. I'm not sure I want to try and interpret this, but it's an interesting result, and a similar thing happens if we try to model a scenario where the person is not only considering their previous flips but also thinking ahead to their next flip:

# Maybe the person was thinking ahead to their next flip as well
moran(seq1, c(1, 2, 0, 1, 0))
## [1] -0.2687347
moran(seq2, c(1, 2, 0, 1, 0))
## [1] 0.0006576715

Some of my application/interpretation of Moran's I to the coin flip problem might be a little off, but it's definitely an applicable measure to use.

A related metric is Geary's C, which is more sensitive to local autocorrelation

$\endgroup$
5
$\begingroup$

When people try to generate random sequences, they tend to avoid repeating themselves more than random processes avoid repeating themselves. Thus, if we look at consecutive pairs of flips, we would expect a human-generated sequence to have too many HT and TH and too few HH and TT compared to a typical random sequence.

The code below explores this hypothesis. It splits each sequence of 300 flips into 150 consecutive pairs and plots the frequency of the four possible results (HH, HT, TH, TT).*

library(tidyverse)

a <- "TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHHTTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHHTTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHTTHHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHTHTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTTHHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT"
b <- "HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTHTHTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHHTTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTTTHTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTHHHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHHHTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT"

# split each sequence of 300 into 150 consecutive pairs
# e.g. TTHHTHTT... -> TT, HH, TH, TT, ...
n_pairs <- 150
ap <- tibble(pair = character(n))
bp <- tibble(pair = character(n))
for (i in 1:n) {
  ap$pair[i] <- substring(a, 2*i - 1, 2*i)
  bp$pair[i] <- substring(b, 2*i - 1, 2*i)
}

# get the frequencies of each possible pair and plot
apc <- count(ap, pair)
bpc <- count(bp, pair)
bind_rows(
  `Sequence 1` = apc,
  `Sequence 2` = bpc,
  .id = 'source') %>%
  ggplot(aes(x = pair, y = n, group = source, fill = source)) + 
  geom_col() + 
  facet_grid(vars(source)) + 
  theme_minimal() +
  geom_hline(yintercept = n_pairs/4, linetype = 'dashed') +
  ylab('frequency') + 
  ggtitle('Frequency of consecutive coin flip pairs')

plot

The dotted line at 150/4 = 37.5 is the expected count of each possible pair assuming the coin flips are independent and fair. By the Law of Large Numbers, we expect the bars not to stray too far from the dotted line. Sequence 1 has an above-average number of HT and TH pairs (especially HT), consistent with our hypothesis about human-generated "randomness". The pairs from Sequence 2 are more consistent with average behavior.

To see how unusual this behavior would be under independent, fair flips, we reformat each sequence's pair count data as a 2x2 contingency table (rows = first flip H/T, columns = second flip H/T) and use Fisher's exact test, which checks whether the data is consistent with a null hypothesis in which the first flip is independent of the second:

for (x in list(apc, bpc)) {
  print(x)
  x %>% 
    mutate(f1 = str_sub(pair, 1, 1),
           f2 = str_sub(pair, 2, 2)) %>%
    select(f1, f2, n) %>%
    spread(f2, n) %>%
    select(-f1) %>%
    as.matrix() %>%
    fisher.test() %>%
    print()
}

The contingency table for Sequence 1 has a p value of 0.0002842, while the table for Sequence 2 has 0.5127. This means that pair frequencies skewed to the degree seen in Sequence 1 would only occur by random 1/0.0002842 = 1 out of 3,519 times, while something like Sequence 2 would be seen very commonly. It seems sensible to conclude that Sequence 1 is the human-made sequence, since its pair frequency table is not consistent with random chance but is quite consistent with the behavior we'd expect of humans.

There is a caveat to this analysis: we do not expect random sequences to be perfectly consistent with average behavior. In fact, in some contexts people know that long random sequences should follow the Law of Large Numbers, and they create sequences which follow it too perfectly. A different analysis would be needed to explore whether Sequence 2 looks odd from this opposite perspective.

* Other answers look at all 299 consecutive pairs, which gives you more data points but they become dependent, which prevents us from using standard significance tests. (For example in the sequence TTHHT, you can make pairs like this: (TT)HHT, T(TH)HT, TT(HH)T, .... This gives you more pairs but consecutive pairs are not independent of one another, as the second flip of a pair determines the first flip of the next pair.)

Alternate Analysis

An analysis that uses all 299 pairs of consecutive flips could be more powerful than the one above if the dependency problem can be solved. To do this, @whuber suggests looking at the transitions between consecutive flips, i.e. when H is followed by T or vice versa. If the flips are independent and fair, then after the first flip, each transition can be considered an independent Bernoulli random variable, and there are 299 transitions total. We can use a two-sided test to see whether the number of transitions observed in each sequence is unlikely under fair independent flips. The transitions are counted and the test applied by the code below:

library(tidyverse)

a <- "TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHHTTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHHTTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHTTHHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHTHTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTTHHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT"
b <- "HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTHTHTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHHTTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTTTHTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTHHHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHHHTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT"

# 299 transitions from flip i to i+1 occur in the sequence of 300
# record these transitions in arrays at and bt
n <- 299
at <- logical(n)
bt <- logical(n)
for (i in 1:n) {
  at[i] <- str_sub(a, i + 1, i + 1) != str_sub(a, i, i)
  bt[i] <- str_sub(b, i + 1, i + 1) != str_sub(b, i, i)
}

# two-sided exact binomial test (analogous to z-test)
# gives probability of transition count more extreme than the one observed
pbinom(sum(at) - 1, n, 1/2, lower.tail = F) + pbinom(n - sum(at), n, 1/2)
pbinom(sum(bt) - 1, n, 1/2, lower.tail = F) + pbinom(n - sum(bt), n, 1/2)

Running this code, we find that Sequence 1 has 204 transitions out of a possible 299. The probability of observing a number of transitions at least this imbalanced, on either the left or the right side, is equal to the probability of observing at least 204 transitions, plus the probability of observing at most 299 - 204 = 95 transitions. This probability is 2.696248e-10, on the order of 3 in 10 billion. Sequence 2 has 153 transitions, and the probability of observing at least 153 transitions or at most 299 - 153 = 146 transitions is 0.7286639. The number of transitions in Sequence 1 is extremely improbable, more so than the 150 pair test above suggested.

$\endgroup$
12
  • $\begingroup$ Your code seems to be doing something more complicated than you describe. It is hard to relate your plots, which are unexplained, to the data: those plots clearly do not report on sequences of length 300. $\endgroup$
    – whuber
    May 7 at 14:43
  • $\begingroup$ I'm splitting each length 300 sequence into 150 pairs of adjacent coin tosses and counting how many of the 150 pairs are HH, HT, TH, or TT. That's what the graphs show. I've added some comments to elaborate on this. If that's not enough, could you be more specific as to what you're struggling to understand? $\endgroup$
    – Paul
    May 7 at 16:55
  • $\begingroup$ I am not struggling to understand anything--I just don't care to have to read through code in order to determine what the content of an answer might be, and I believe most readers will feel the same. $\endgroup$
    – whuber
    May 7 at 17:25
  • 1
    $\begingroup$ Thank you: what you are doing is now more apparent. But why? Could you explain what is accomplished with this split? $\endgroup$
    – whuber
    May 7 at 18:10
  • 1
    $\begingroup$ I've added a transition-based analysis and it is indeed a more powerful test. Great idea. $\endgroup$
    – Paul
    May 10 at 0:45
5
$\begingroup$

This answer is inspired by @user1717828's answer which transforms the sequence of coin flips into a random walk. I don't show the two given sequences as random walks here; see @user1717828's answer for that plot.

The random walk approach is interesting because it examines long-run features of the sequence rather than short-range ones (such as the overabundance of H-T and T-H flips). In a way the two approaches are complementary as they analyze the sequences at different scales. At the "local" scale the fake sequence oscillates which creates autocorrelation; at the "global" scale it keeps within a limited range of positions for long stretches of time. Both of these aspects are non-random and as one occurs, so does the other: staying in place means not moving around and vice versa.

A (one-dimensional) simple random walk $S_n$, defined as the sum of n random variables that take value {-1, 1} with probability ½, has many interesting properties, including:

  • The probability that a simple random walk returns to the origin is 1. In fact, it visits every integer infinitely often.
  • The mean of a simple random walk is $\text{E}(S_n) = 0$ and its variance is $\text{Var}(S_n) = n$.
  • A random walk with nonzero mean (which corresponds to flipping a biased coin) is transient: it makes finitely many visits to the origin before diverging, to +∞ if the mean is positive and to -∞ in the mean is negative.

These properties imply that a simple random walk spends a lot of time far away from 0 before eventually returning to it.

We can use this intuition founded on theory (as well as any other property of random walks) to propose characteristics for comparing the two sequences. Then we generate many simple random walks and use the observed distribution of the features to estimate how "extreme" the two given sequences are.

To capture the variance of a simple random walk I use the maximum distance from 0. To capture its tendency to visit every integer I use the number of times the walk crosses the integers from -8 to +8. To make a symmetric grid of scatterplots and to avoid suspicions of HARKing I don't show the crossings of 0. It's important to look at crossing on both sides of 0 (the starting point), so that we don't make assumptions about how the non-randomness manifests in the data. For example, we don't know whether the coin is fair (p = ½), biased towards heads (p > ½) or biased towards tails (p < ½).

Note: HARKing is the practice of performing many analyses to find an interesting hypothesis. I compute many statistics, grounded in the theory of random walks, and report all of them, thus avoiding HARKing.

Here are the results from the simulation of 500 simple random walks. Sequence #1 (in blue) appears as an outlier compared to sequence #2 (in red) in many panels. Due to its overabundance of H-T and T-H pairs, sequence #1 doesn't explore enough and spends too much time between -5 and -8 (shown in the top row). Sequence #1 doesn't advance beyond the band [-9, 9] and it is rare for a true random walk so stay so close to 0 for 300 steps. Visually, sequence #1 is the overall outlier even if it doesn't appear extremely "unusual" is some panels.

enter image description here

$\endgroup$
6
  • 2
    $\begingroup$ +1. This explains why I did not examine this statistic in my answer: I had done so, exactly as you have, and came to the same conclusion that the max deviation statistic was suggestive but only weak evidence of a departure from randomness. But in that investigation I noted there were other visual quirks in the plot of the first series and that led to a more effective statistic for assessing the randomness of the series. It's worth noting that extensive computation is unnecessary: you could draw all your conclusions more quickly by simulating as few as 20 random walks.. $\endgroup$
    – whuber
    May 9 at 13:44
  • $\begingroup$ @whuber After some (over)thinking, I realized that it works quite well to compute two statistics: one to represent central tendency and the other -- extreme values. Similar to COOLSerdash's solution which looks at number of runs (a kind of average) and length of longest run. With the simple random walk, zero crossings vs maximum distance from zero works very well (visually). But then I don't know how to get a p-value. $\endgroup$
    – dipetkov
    May 9 at 21:35
  • $\begingroup$ I think the zero crossings observation comes down to HARKing. (For instance, these could well be realizations of Markov processes in which a transition away from a side of the coin is made with probability $p$; $p\approx 2/3$ in one sequence and $p\approx 1/2$ in the other. In general I wouldn't expect zero-crossing counts to distinguish among these.) A better way to approach this would be to do all your exploration on the first half of each sequence; develop a suitable statistic from that; and apply it to the second halves of the sequences. $\endgroup$
    – whuber
    May 9 at 21:39
  • $\begingroup$ @whuber I don't agree it's HARKing. At least not any more than number of runs. A simple random walk returns to the start infinitely often, so given this fact it makes sense to look at zero crossings. Since the fake sequence doesn't move enough any number will do in theory; here because sequences are only 300 long every number between about -10 to 10 will work. $\endgroup$
    – dipetkov
    May 9 at 22:03
  • $\begingroup$ Ask yourself this: would you have focused on counting zero crossings by exploring just the first halves of the sequences? You might have if both of two things occurred: (1) you observed no zero crossings in one and many zero crossings in the other and (2) in simulated iid sequences, you observed a consistent tendency to have positive numbers of zero crossings. In such a case, you likely have a useful statistic. Unfortunately, you would discover that about 6.5% of the time, a Binomial random walk has no zero crossings. That's why I view this as HARKing: it's accidental. $\endgroup$
    – whuber
    May 10 at 14:32
1
$\begingroup$

In addition to statistical approaches, one visual approach is to plot the sequences as a "drunkards walk". Treat H as a step forwards and T as a step back and plot the sequences. One way in Python is:

import altair as alt
import pandas as pd

seq1 = "HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTHTHTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHHTTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTTTHTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTHHHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHHHTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT"

seq2 = "TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHHTTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHHTTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHTTHHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHTHTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTTHHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT"


def encode(stream):
    return [1 if c == "H" else -1 for c in stream]


df = pd.DataFrame({"seq1": encode(seq1), "seq2": encode(seq2)})
df_cumsum = pd.melt(
    df.cumsum(), var_name="sequence", value_name="cumsum", ignore_index=False
).reset_index()

chart = alt.Chart(df_cumsum).mark_line().encode(x="index", y="cumsum", color="sequence")
chart.show()

enter image description here

When comparing a truly random sequence to a human generated one, it is often easy to tell them apart based on inspection of the walk.

$\endgroup$
6
  • 4
    $\begingroup$ This is a good approach. It was the first thing I tried. However, I created a context for evaluating these random walks: I also generated 18 more sequences using iid flips of a fair coin and plotted all 20. Generally, whenever I repeated this procedure, there was a truly random walk that looked like the orange one here, at least in terms of the relatively small range. Exactly which characteristic(s) suggest to you that the orange curve is not random but the blue is? $\endgroup$
    – whuber
    May 8 at 12:53
  • $\begingroup$ The orange one is what I think is the truly random one based on this visual. I would have guessed based on just the number of zero-crossings alone that blue is human-generated. But you're right, there's always going to be a chance that any sequence is randomly-generated, and the analytical methods are better for quantifying the odds of that. $\endgroup$ May 8 at 20:17
  • 1
    $\begingroup$ Unfortunately for your theory, the orange one is the non-random one. It has an unusually low range; but more to the point, its thickened appearance is due to a huge overabundance of H-T-H oscillations. $\endgroup$
    – whuber
    May 8 at 20:32
  • $\begingroup$ Ahh! In that case, I'm going to rethink my approach of visualizing these things before looking at the statistics next time! $\endgroup$ May 8 at 22:29
  • 1
    $\begingroup$ The thickened appearance is the feature analyzed in the answers by whuber, Dmetri Pananos and Paul. The distance from 0 in a simple random walk (which I think is what this answer is getting at) would be a different feature altogether. A (true) simple random walk will spend a lot of time away from zero until eventually returning (with probability 1). The sequences in this exercise though seem too short for this kind of argument. $\endgroup$
    – dipetkov
    May 8 at 23:55
1
$\begingroup$

Looking at the HH, HT, TH, TT frequencies is probably the most straightforward way to approach the two series presented, given people's tendency to apply HT and TH more frequently when trying to appear random. More generally, however, that approach will fail to detect non-randomness even in sequences with obvious patterns. For instance, a repeating sequence of HHHTHTTT will produce balanced counts of pairs as well as triples (HHH, HHT, etc.).

Here's an idea for a more general solution that was initially inspired by answers discussing random walks. Start with the observation that for a random sequence, the number of heads in any given subsequence is binomially distributed. We can count the number of heads, $n_{i,j}$, between flip number $i$ and flip number $j>i$ for all values of $i,j$. Then compare these to the expected number of counts if all the $n_{i,j}$ were independent. Although they are obviously not independent, the comparison gives rise to a useful statistic: the maximum absolute difference between the observed counts and the expected counts.

Applying this approach to sequence 1 gives us 98.26 as our test statistic: there are 257 subsequences of length 44. If they were all composed of independent Bernoulli trials, the expected number of the 257 that contained exactly 22 heads is ~30.74, whereas sequence 1 contains 129 subsequences with exactly 22 heads (very underdispersed). 129 - 30.74 = 98.26, which is the maximum of these differences for sequence 1.

Performing the same calculations on sequence 2 gives a test statistic of 48.30: there are 197 subsequences of length 104. The expected number containing exactly 54 heads would be ~13.70. Sequence 2 contains 62, so the test statistic is 62 - 13.70 = 48.30.

The test statistics can be compared to those from a large number of random sequences of the same size. In this case, no samples are greater than the test statistic from sequence 1, and about 14% of samples are greater than the statistic from sequence 2.

Here it is all together with R code:

library(parallel)
set.seed(16)

seq1 <- utf8ToInt("TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHHTTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHHTTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHTTHHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHTHTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTTHHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT") == utf8ToInt("H")
seq2 <- utf8ToInt("HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTHTHTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHHTTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTTTHTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTHHHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHHHTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT") == utf8ToInt("H")

# function to calculate test statistic S 
fS <- function(m, p = 0.5) {
  if (!is.matrix(m)) m <- matrix(m, ncol = 1)
  n <- nrow(m)
  steps <- rep.int(1:(n - 1L), 2:n)
  nEx <- dbinom(sequence(2:n, from = 0L), steps, p)*(n - steps)
  idx <- sequence((n - 1L):1)
  idx <- idx*(idx + 1L)/2L
  S <- numeric(ncol(m))
  for (i in 1:ncol(m)) S[i] <- max(abs(nEx - tabulate(idx + dist(cumsum(m[,i])), length(nEx))))
  S
}

(S1 <- fS(seq1))
#> [1] 98.26173
(S2 <- fS(seq2))
#> [1] 48.30014
# calculate S from 1e6 random sequences of length n (probably overkill)
n <- length(seq1)
cl <- makeCluster(detectCores() - 1L)
clusterExport(cl, list("fS", "n"))
system.time(simS <- unlist(parLapply(cl, 1:100, function(i) fS(matrix(sample(0:1, 1e4L*n, TRUE), n)))))
#>    user  system elapsed 
#>    0.00    0.03  339.42
stopCluster(cl)
nsim <- 1e6L
# calculate approximate p-values
sum(simS > S1)/nsim
#> [1] 0
sum(simS > S2)/nsim
#> [1] 0.139855
max(simS)
#> [1] 91.01078
$\endgroup$
0
$\begingroup$

This problem specifies that we have the following information:

  • We observed two coin flip sequences: sequence $S_1$ and sequence $S_2$.
  • Each of these sequences could have been generated by either mechanism $R$ corresponding to independent flips of a fair coin or some other mechanism $\bar{R}$.
  • Exactly one of these two sequences was generated by the mechanism $R$.

Since both sequences $S_1$ and $S_2$ could possibly have been generated by $R$ or $\bar{R}$, we cannot be sure which of these two sequences was generated by $R$. However, probability theory provides us with the tools to quantify our belief $P(R_i|S_i)$ that sequence $S_i$ was generated by mechanism $R$. Then, given that we have been provided with the information $I$ that exactly one of these two sequences has been generated by mechanism $R$, then the probability that sequence 1 was generated by mechanism $\bar{R}$ and sequence 2 was generated by mechanism $R$ is:

$$ P(\bar{R_1}R_2|S_1 S_2 I) = \frac{P(R_2|S_2)}{P(R_1|S_1) + P(R_2|S_2)} $$

This problem explicitly specifies what it means for a sequence to have been generated by mechanism $R$: the sequence $S_i$ is described by independent flips $y_{ij}$ of a fair coin, such that the likelihood is:

$$ P(S_i|R_i) = \prod_j \mathrm{Bernoulli}(y_{ij} | 0.5) $$

However, in order to determine the probability that a given sequence was generated by mechanism $R$, $P(R_i|S_i)$, we also need to specify alternative mechanisms $\bar{R}$ by which the sequences may have been generated. At this point we have to use our own creativity and external information to develop models that could reasonably describe how these sequences might have been generated, if they were not generated by mechanism $R$.

The other answers to this question do a great job at describing various mechanisms $\bar{R}$ that describe how sequences could be generated:

  • Sequences are generated by sampling pairs of coin flips from a distribution where the probability of different pairs of coin flips deviates from uniform: 1 2 3
  • Sequences are generated in a way such that they do not have long runs of only heads or tails: 1
  • Sequences are generated in such a way so that a given compression algorithm results in a small compressed size: 1
  • Sequences are generated by a random walk such that a given coin flip outcome is affected by recent coin flip outcomes: 1 2 3

Note: we can always develop a hypothesis that appears to result in our observed sequence having been "inevitable", according to our likelihood, and if our prior probability for that mechanism of sequence generation is high enough, then using this hypothesis can always result in a low probability that the sequence had been generated by mechanism $R$. This practice is also referred to as HARKing, and Jaynes refers to this as a "sure-thing hypothesis". In general, we may have multiple hypotheses that can describe how a sequence may be generated, but the more contrived hypotheses such as "sure-thing hypotheses" should generally have sufficiently low prior probabilities such that we generally would not infer that those hypotheses have the highest probability of having generated a given sequence.

Here's an example using this procedure with one possible model for $\bar{R}$ sequence generation to quantify the probability that sequence 2 was the sequence that was generated by mechanism $R$. This stan model defines model_r_lpmf as the log likelihood for sequence generation mechanism $R$, and model_not_r_lpmf as the log likelihood for sequence generation mechanism $\bar{R}$. In this case, we are testing a hypothesis $\bar{R}$ in which coin flips can be correlated, and there is a fairly well defined correlation length. Then, we are modeling the probability that the sequence was generated by model $R$ as p_model_r, such that our prior probability considers mechanisms $R$ and $\bar{R}$ sequence generation equally probable.

// model for hypothesis testing whether a sequence of coin flips is generated by either:
// - model R: independent flips of a fair coin, or
// - model not R: unfair coin with correlations

functions {
  // independent flips of a fair coin
  real model_r_lpmf(int[] sequence, int N) {
    return bernoulli_lpmf(sequence | 0.5);
  }

  // unfair coin with correlations
  real model_not_r_lpmf(int[] sequence, int N, int n, vector beta) {
    real tmp = 0;

    real offset = 0;
    for (i in 1:N) {
      offset = beta[1];
      for (j in 1:min(n - 1, i - 1)) {
          offset += beta[j + 1] * (2 * sequence[i - j] - 1);
      }
      tmp += bernoulli_logit_lpmf(sequence[i] | offset);
    }
    return tmp;
  }
}
data {
  int N;                                 // the length of the observed sequence
  int<lower=0, upper=1> sequence[N];     // the observed sequence
}
transformed data {
  int n = N / 2 + 1;                     // number of correlation parameters to model
}
parameters {
  real<lower=0, upper=1> p_model_r;       // probability that the sequence was generated by model R
  // not R model parameters
  real log_scale;                        // log of the scale of the correlation coefficients `beta`
  real log_decay_rate;                   // log of the decay rate for the correlation coefficients `beta`
  vector[n] alpha;                       // pretransformed correlation coefficients
}
transformed parameters {
  // transformed not R model parameters
  real scale = exp(log_scale);           // typical scale of the correlation coefficients
  real decay_rate = exp(log_decay_rate); // typical decay rate of correlation coefficients
  vector[n] beta;                        // correlation coefficients
  for (j in 1:n) {
    beta[j] = alpha[j] * scale * exp(-(j - 1) * decay_rate);
  }
}
model {
  // uninformative haldane prior for the probability that the sequence was generated by model R
  target += - log(p_model_r) - log1m(p_model_r);

  // priors for not R model parameters
  log_scale ~ normal(-1, 1);
  log_decay_rate ~ normal(-1, log(n + 1) / 2.0);
  alpha ~ normal(0, 1);

  // the sequence is was generated by model R with probability `p_model_r`, otherwise it was generated by model not R
  target += log_mix(
    p_model_r,
    model_r_lpmf(sequence | N),
    model_not_r_lpmf(sequence | N, n, beta)
  );
}

Using pystan with this model, we can evaluate the probability p_model_r that each of the two provided sequences were generated by mechanism $R$.

model = pystan.StanModel(model_code=MODEL_CODE)

sequence_1 = "TTHHTHTTHTTTHTTTHTTTHTTHTHHTHHTHTHHTTTHHTHTHTTHTHHTTHTHHTHTTTHHTTHHTTHHHTHHTHTTHTHTTHHTHHHTTHTHTTTHHTTHTHTHTHTHTTHTHTHHHTTHTHTHHTHHHTHTHTTHTTHHTHTHTHTTHHTTHTHTTHHHTHTHTHTTHTTHHTTHTHHTHHHTTHHTHTTHTHTHTHTHTHTHHHTHTHTHTHHTHHTHTHTTHTTTHHTHTTTHTHHTHHHHTTTHHTHTHTHTHHHTTHHTHTTTHTHHTHTHTHHTHTTHTTHTHHTHTHTTT"
samples_1 = model.sampling(
    data=dict(N=len(sequence_1), sequence=[int(c == "H") for c in sequence_1]),
    chains=16,
)
p_model_r_1 = samples_1["p_model_r"].mean()
# p_model_r_1 = 0.027

sequence_2 = "HTHHHTHTTHHTTTTTTTTHHHTTTHHTTTTHHTTHHHTTHTHTTTTTTHTHTTTTHHHHTHTHTTHTTTHTTHTTTTHTHHTHHHHTTTTTHHHHTHHHTTTTHTHTTHHHHTHHHHHHHHTTHHTHHTHHHHHHHTTHTHTTTHHTTTTHTHHTTHTTHTHTHTTHHHHHTTHTTTHTHTHHTTTTHTTTTTHHTHTHHHHTTTTHTHHHTHHTHTHTHTHHHTHTTHHHTHHHHHHTHHHTHTTTHHHTTTHHTHTTHHTHHHTHTTHTTHTTTHHTHTHTTTTHTHTHTTHTHTHT"
samples_2 = model.sampling(
    data=dict(N=len(sequence_2), sequence=[int(c == "H") for c in sequence_2]),
    chains=16,
)
p_model_r_2 = samples_2["p_model_r"].mean()
# p_model_r_2 = 0.790

Here we have quantified the probability that each sequence was generated according to the sequence generating model $R$ as the following:

$$ \begin{aligned} P(R_1|S_1) &= 0.027 \\ P(R_2|S_2) &= 0.790 \\ P(\bar{R_1}R_2|S_1 S_2 I) &= 0.967 \\ \end{aligned} $$

So, we are fairly certain that sequence 2 was generated by independent coin flips and sequence 1 was not, relative to the model of sequence generation $\bar{R}$ that we defined, and we quantify this belief with a probability of 0.967. This is a high probability, but of course we could still be wrong -- both sequences could have been generated by either model. Additionally, we could have selected a model of sequence generation for $\bar{R}$ that can not describe the procedure that was actually used to generate the sequences. This process for quantifying our belief that one of these two sequences was generated by a fair coin with independent flips is the piece of reasoning that currently appears to be missing from the other answers that are currently available.

$\endgroup$
2
  • $\begingroup$ I have to reject your initial dichotomy, because it (a) is too vague to be actionable (what can one do, in any objective and quantifiable way, to characterize any "non-random process"?) and (b) ignores a huge set of other possible mechanisms: namely, random mechanisms that are not iid.. In fact, it is almost surely the case that whoever generated these sequences used a random mechanism for both--they are just different random mechanisms. $\endgroup$
    – whuber
    May 31 at 12:00
  • $\begingroup$ @whuber I framed this problem in terms of hypothesis testing between a model $R$ by which sequences are generated by independent flips of a fair coin, and models of other possible ways that we might believe that sequences could be generated $\bar{R}$. I used labels "random", and "not random" to describe these models, which appears to have been confusing, so I have edited the post to remove that terminology. $\endgroup$
    – b-r-oleary
    Jun 8 at 3:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.