How can I determine which of two sequences of coin flips is real and which is fake?

Question

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.

Answer 1

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.

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.

Answer 2

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

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()
```

Answer 3

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.

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)

Answer 4

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.

Answer 5

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.

Answer 6

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

Answer 7

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')

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.

Answer 8

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.

Answer 9

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()

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.

Answer 10

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

Answer 11

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.