# Could any equation have predicted the results of this simulation?

Suppose there is a coin that has a 5% chance of landing on HEADS and a 95% chance of landing on TAILS. Based on a computer simulation, I want to find out the following :

• The average number of flips before observing HEADS, TAILS, HEADS (note: not the probability, but the number of flips)

Using the R programming language, I tried to make this simulation (this simulation keeps flipping a coin until HTH and counts the number of flips until this happens - it then repeats this same process 10,000 times):

results <- list()

for (j in 1:10000) {

response_i <- ''
i <- 1

while (response_i != 'HTH') {
response_i <- c("H","T")
response_i <- sample(response_i, 3, replace=TRUE,
prob=c(0.05, 0.95))
response_i <- paste(response_i, collapse = '')

iteration_i = i
if (response_i == 'HTH') {
run_i = data.frame(response_i, iteration_i)
results[[j]] <- run_i
}
i <- i + 1
}
}

data <- do.call('rbind', results)


We can now see a histogram of this data:

hist(data$iteration_i, breaks = 500, main = "Number of Flips Before HTH")  We can also see the summary of this data: summary(data$iteration_i)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
1.0   119.0   288.0   413.7   573.0  3346.0


My Question:

• Could any "mathematical equation" have predicted the results of this simulation in advance? Could any "formula" have shown that the average number of flips to get HTH would have been 413? Can Markov Chains be used to solve this problem?

• Based on the "skewed" shape of this histogram, is the "arithmetical mean" (i.e. mean = sum(x_i)/n) a "faithful" representation of the "true mean"? Looking at the above histogram, we can clearly see that you are are more likely to see HTH before 437 iterations compared to seeing HTH after 437 iterations, e.g. (on 100,000 simulations, the new average is 418):

nrow(data[which(data$iteration_i <418), ]) 63184 nrow(data[which(data$iteration_i > 418), ])

36739

For such distributions, is there a better method to find out the "expectation" of this experiment?

Thanks!

• Could you please use sentence case for the titles of your questions? It makes them uncomfortable to read otherwise. I've updated the title of your question to use sentence case. Nov 23, 2021 at 14:21
• This is Penney's Game with an unfair coin. The solution methods at stats.stackexchange.com/questions/12174 all apply, with suitable changes to the chances.
– whuber
Nov 23, 2021 at 20:32
• At the question @whuber links to, I gave an answer of $2^3+2$ using Conway's algorithm. With the chances changed, this becomes $\frac{20}{1}\times \frac{20}{19}\times \frac{20}{1} +\frac{20}{1} =441+\frac1{19}$ Nov 24, 2021 at 14:38
• Sanity check for those calculating: I crunched 1 billion runs, and got a value of 441.0557, and an absolutely beautiful exponential distribution. Nov 25, 2021 at 13:26
• Unless I'm missing something, your code doesn't match what you seem to be trying to do, and what all the answers assume you want to do. Your code counts the number of triple-flips until one triple-flip is HTH; so, for example, if the sequence of flips is THT HTT TTT TTH TTT HTT HTH, your code will give '7' (= the seventh triple-flip is HTH), but I think what you would really want is '4' (= the second through fourth flips are HTH). The two averages are actually pretty close -- 421 vs. 441 or so -- but that's basically a coincidence. Nov 26, 2021 at 7:03

At any given point in the game, you're $$3$$ or fewer "perfect flips" away from winning.

For example, suppose you've flipped the following sequence so far: $$HTTHHHTTTTTTH$$

You haven't won yet, but you could win in two more flips if those two flips are $$TH$$. In other words, your last flip was $$H$$ so you have made "one flip" worth of progress toward your goal.

Since you mentioned Markov Chains, let's describe the "state" of the game by how much progress you have made toward the desired sequence $$HTH$$. At every point in the game, your progress is either $$0$$, $$1$$, or $$2$$--if it reaches $$3$$, then you have won. So we'll label the states $$0$$, $$1$$, $$2$$. (And if you want, you can say that there's an "absorbing state" called "state $$3$$".)

You start out in state $$0$$, of course.

You want to know the expected number of flips, from the starting point, state $$0$$. Let $$E_i$$ denote the expected number of flips, starting from state $$i$$.

At state $$0$$, what can happen? You can either flip $$H$$, and move to state $$1$$, or you flip $$T$$ and remain in state $$0$$. But either way, your "flip counter" goes up by $$1$$. So: $$E_0 = p (1 + E_1) + (1-p)(1 + E_0),$$ where $$p = P(H)$$, or equivalently $$E_0 = 1 + p E_1 + (1-p) E_0.$$ The "$$1+$$" comes from incrementing your "flip counter".

At state $$1$$, you want $$T$$, not $$H$$. But if you do get an $$H$$, at least you don't go back to the beginning--you still have an $$H$$ that you can build on next time. So: $$E_1 = 1 + p E_1 + (1-p) E_2.$$

At state $$2$$, you either flip $$H$$ and win, or you flip $$T$$ and go all the way back to the beginning. $$E_2 = 1 + (1-p) E_0.$$

Now solve the three linear equations for the three unknowns.

In particular you want $$E_0$$. I get $$E_0 = \left( \frac{1}{p} \right) \left( \frac{1}{p} + \frac{1}{1-p} + 1 \right),$$ which for $$p=1/20$$ gives $$E_0 = 441 + 1/19 \approx 441.0526$$. (So the mean is not $$413$$. In my own simulations I do get results around $$441$$ on average, at least if I do around $$10^5$$ or $$10^6$$ trials.)

In case you are interested, our three linear equations come from the Law of Total Expectation.

This is really the same as the approach in Stephan Kolassa's answer, but it is a little more efficient because we don't need as many states. For example, there is no real difference between $$TTT$$ and $$HTT$$--either way, you're back at the beginning. So we can "collapse" those sequences together, instead of treating them as separate states.

Simulation code (two ways, sorry for using Python instead of R):

# Python 3
import random

def one_trial(p=0.05):
state = 0 # states are 0, 1, 2, representing the progress made so far
flip_count = 0 # number of flips so far
while True:
flip_count += 1
if state == 0: # empty state
state = random.random() < p
# 1 if H, 0 if T
elif state == 1: # 'H'
state += random.random() >= p
# state 1 (H) if flip H, state 2 (HT) if flip T
else: # state 2, 'HT'
if random.random() < p: # HTH, game ends!
return flip_count
else: # HTT, back to empty state
state = 0

def slow_trial(p=0.05):
sequence = ''
while sequence[-3:] != 'HTH':
if random.random() < p:
sequence += 'H'
else:
sequence += 'T'
return len(sequence)

N = 10**5
print(sum(one_trial() for _ in range(N)) / N)
print(sum(slow_trial() for _ in range(N)) / N)

• Nice, this is the best answer so far. Nov 23, 2021 at 20:30
• Agreed, I was also wondering why people just didn't bring out the ol' Markov and chain... Nov 24, 2021 at 20:25
• Very neat state machine solution! Reminds me of the many varieties of best time to buy and sell stocks. Dec 10, 2021 at 6:52

First, you can refactor your R code to be (IMHO) a little more legible, also using pbapply::pbreplicate() to get a nice progress bar:

n_sims <- 1e5
library(pbapply)

results <- pbreplicate(n_sims,{
flips <- NULL
while(length(flips)<3 || !identical(tail(flips,3),c("H","T","H"))){
flips <- c(flips,sample(c("H","T"),size=1,prob=c(.05,.95)))
}
length(flips)
})

hist(results,breaks=seq(-0.5,max(results)+0.5))


Growing the flips vector in each step makes this code very slow. It would be much faster to grow it in large batches - but that would be at the expense of legibility. As it is, this code runs its 100,000 replicates in about ten minutes, which we can use to google around a bit for an abstract solution.

Second, there indeed is an abstract solution. The details are a bit complicated, and I unfortunately don't have the time to write them all up right now, but I'll give the gist.

Specifically, we can model your sequence of unfair coin flips as a Markov chain, where each state corresponds to the last three flips, and where we stop after hitting HTH. Thus, we have a single absorbing state and seven transient states, and if we use $$p$$ to denote the probability of next getting a Head, we get a state transition matrix as follows:

(Apologies for pasting a picture; I don't know how to get the annotations to work in MathJaX. LaTeX code below.)

This is already in canonical form, with the single absorbing state at the bottom and the right. We can partition the matrix into a block matrix $$T$$ for the transient states and the row and column corresponding to the absorbing state: $$\begin{pmatrix} T & T_0 \\ 0^t & 1 \end{pmatrix}.$$

The initial states are all equally probable - this is our initial distribution $$\tau$$.

The distribution of the number of steps necessary to reach the absorbing state is the discrete phase-type distribution, which depends on the initial state $$\tau$$ and the submatrix $$T$$. You can take a look at Dayar (2005) on how to calculate the expectation of this distribution (also see Expected number of steps between states in a Markov Chain, where I got the pointer to this paper), but note that (1) you already start after three steps, so you would need to add $$3$$ to the expectation, and (2) Dayar (2005) assumes that the probability of starting in the absorbing state is zero, and in your case, it isn't, it's already $$p^2(1-p)$$, so you need to correct for that.

LaTeX code for the matrix, kudos to Bernard:

\documentclass{article}
\usepackage{mathtools}
\usepackage{blkarray, bigstrut}

\begin{document}

$\setlength{\bigstrutjot}{4pt} $$\begin{blockarray}{r@{\enspace\vrule}rcccccccc} \phantom{HHH}& & HHT & HHH & THT & THH & TTT & TTH & HTT & HTH \\ \BAhline \begin{block}{r@{\enspace\vrule}r(cccccccc)} HHT & & 0 & 0 & 0 & 0 & 0 & 0 & 1-p & p \bigstrut[t]\\ HHH & & 1-p & p & 0 & 0 & 0 & 0 & 0 & 0 \\ THT & & 0 & 0 & 0 & 0 & 0 & 0 & 1-p & p \\ THH & & 1-p & p & 0 & 0 & 0 & 0 & 0 & 0 \\ TTT & & 0 & 0 & 0 & 0 & 1-p & p & 0 & 0 \\ TTH & & 0 & 0 & 1-p & p & 0 & 0 & 0 & 0 \\ HTT & & 0 & 0 & 0 & 0 & 1-p & p & 0 & 0 \\ HTH & & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\bigstrut[b]\\ \end{block} \end{blockarray}$$$

\end{document}

• Can you include the final formula in your answer as well instead of just linking to Math SO answer? Good find btw! Nov 23, 2021 at 11:28
• @justhalf: I would love to, but as I wrote I don't really have the time to dig through that paper and also wrap my head around the caveats I noted. Also, before putting up a formula someone may use without checking, I would really prefer to double-check it with a couple of simulations, again something that would take some time. I'll see whether I find the time to update this in the next days. Nov 23, 2021 at 11:32
• I think you can simply the matrix to 5x5 by considering last two steps only, plus a terminal/absorbing state i.e. HH, HT, TH, TT, HTH. HT would either transition to TT (1-p) or HTH (p). Nov 23, 2021 at 18:34
• @NeilSlater Yes, and the initial state is really TT (no progress toward goal). But then there's also no real difference between states HH and TH: both amount to "You have an H so far". So actually 3 transient states suffice, plus an absorbing state if you want. (See my answer for more details.) Nov 23, 2021 at 19:47
• I agree, @mathmandan's answer is shorter, simpler and more elementary, very nice (+1). Nov 24, 2021 at 6:38

There is a fun way to answer this problem using martingales, and in particular using https://en.wikipedia.org/wiki/Optional_stopping_theorem. I first saw this trick in the book A First Look at Rigorous Probability Theory by Jeffrey S. Rosenthal, in the martingale chapter. (I don't have the book in front of me at the moment but I'll edit and add a page or equation number when I do.)

Imagine that at each time step $$t$$ a person arrives and places a bet on the outcome of the next coin toss. The payoff to their first bet is:

$$$$\begin{cases} -1 &\mbox{if toss } t \mbox{ is T} \\ +19 &\mbox{if toss } t \mbox{ is H} \\ \end{cases}$$$$

If they are wrong (i.e. lose money) on this bet, they stop betting. Note that their net gain if they exit at this point is $$-1$$. If they are correct, they continue to betting on toss $$t+1$$ and receive the following payoff:

$$$$\begin{cases} +\frac{20}{19} &\mbox{if toss } t+1 \mbox{ is T} \\ -20 &\mbox{if toss } t+1 \mbox{ is H} \\ \end{cases}$$$$

If they are wrong on this second bet, they stop betting. Again, note that their net gain if they exit at this point is (by design) $$19-20 = -1$$. If they were correct for this second bet, they proceed to betting on the outcome of toss $$t+2$$ and receive:

$$$$\begin{cases} -\frac{400}{19} &\mbox{if toss } t+2 \mbox{ is T} \\ +400 &\mbox{if toss } t+2 \mbox{ is H} \\ \end{cases}$$$$

If they are incorrect, they exit with a net gain of $$19 + \frac{20}{19} - \frac{400}{19} = -1$$. If they are correct, everything stops: we have just seen a HTH sequence, and the person who started betting at the beginning of that sequence has just won $$19 + \frac{20}{19} + 400 = 420.0526$$.

Note that two people have started betting after this big winner: one person whose first bet was incorrect (they exit with $$-1$$), and a second person whose first bet was correct (they win $$19$$ but nothing more because the process stops).

Let $$\tau$$ denote the stopping time, and let $$X_t$$ denote the cumulative net amount won by all gamblers up to and including time $$t$$. Let $$X_0 = 0$$, and note that $$X_t$$ is a martingale because all bets are (by construction) fair, with an expected payoff of $$0$$. This will let us use https://en.wikipedia.org/wiki/Optional_stopping_theorem and in particular $$\mathbf{E}[X_\tau] = \mathbf{E}[X_0] = 0$$.

$$X_\tau$$ is the total amount won when we reach the stopping time. At this point $$\tau - 3$$ people will have lost their bets and exited with a net gain of $$-1$$, one person will have won $$420.0526$$, and we also have to account for the last two people who start betting after the winner. We have:

$$$$\mathbf{E}[X_\tau] = 0 = (\mathbf{E}[\tau]-3)*(-1) + 420.0526 + (-1) + 19$$$$

Which leads to $$\mathbf{E}[\tau] = 420.0526 + 3 + 18 = 441.0526$$, which agrees with the answers posted earlier.

Here is a somewhat clumsy brute-force method to obtain the probabilities and order statistics. Getting the mean will take more work.

So first just generate the possible sequences and associated probabilities where "HTH" are the last 3 flips (with that sequence not occurring previously). Then look for patterns. For integer patterns the go to place is http://oeis.org.

Note: To make some of the calculations easier I've used 1 for H and 0 for T. Using Mathematica the probabilities for the flip that ends with "HTH" (or "101") with general $$p$$ is generated as follows:

s[1] = {{1, 0, 1}};
pr[1] = p^2 (1 - p);
Print[{1, pr[1]}]
Do[s[i] = Select[Flatten[{Prepend[#, 1], Prepend[#, 0]} & /@ s[i - 1], 1], ! (#[[1]] == 1 && #[[2]] == 0 && #[[3]] == 1) &];
pr[i] = Total[p^Total[#] (1 - p)^(Length[#] - Total[#]) & /@ s[i]],
{i, 2, 8}]
Table[{i, pr[i]}, {i, 1, 8}] // TableForm


The pattern of the powers of $$p$$ and $$1-p$$ seems obvious but maybe not the associated coefficients. That's where http://oeis.org comes in. We plug in the coefficients for $$n=8$$

1,6,12,13,9,6,1,1


and find sequence A124279. That translates to the formula for the probability of flip $$n$$ ending in "HTH" (and not containing any previous "HTH" sequence):

pr[n_] := p^2 (1 - p) Sum[p^(n - k - 1) (1 - p)^k Sum[Binomial[j, k - j] Binomial[n - k - 1, k - j], {j, 0, n}], {k, 0, n - 1}]


or

$$pr(n) = p^2 (1-p) \sum _{k=0}^{n-1} (1-p)^k p^{-k+n-1} \sum _{j=0}^n \binom{j}{k-j} \binom{-k+n-1}{k-j}$$

The median is between flip 303 and 304 as the associated cumulative probabilities are 0.49891 and 0.500051, respectively, when $$p=0.05$$.

To calculate the probabilities in R you'll either need to use multiple precision arithmetic or reduce round-off errors by using logs.

Disclosure: I wrote the samc R package used in this answer

In it, he showed how to construct a transition matrix representing the problem: Image credit: Stephan Kolassa's answer

Now, as indicated in comments/another answer, this matrix can be simplified a variety of different ways, but I'm going to keep it as-is because there are other things about the experiment you can learn from it that you can't from the others.

Let's call this transition matrix $$P$$, which is how I define it in the samc package: samc overview

This can be broken down as: $$P = \begin{bmatrix} Q & R \\ 0 & I \end{bmatrix}$$

Using this, there a variety of things you can calculate about the model. For example: $$z=(I-Q)^{-1}{\cdot}1=F{\cdot}1$$ Will tell you the expected time to absorption given a starting transient state. In this example, that means that given your last 3 flips, how many more flips do you expect it to take before $$HTH$$ occurs. Some code:

library(samc)

p <- 0.05
q <- 1 - p

p_mat <- matrix(c(0, 0, 0, 0, 0, 0, q, p,
q, p, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, q, p,
q, p, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, q, p, 0, 0,
0, 0, q, p, 0, 0, 0, 0,
0, 0, 0, 0, q, p, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1),
8, byrow = TRUE)

rownames(p_mat) <- c("HHT", "HHH", "THT", "THH", "TTT", "TTH", "HTT", "HTH")
colnames(p_mat) <- c("HHT", "HHH", "THT", "THH", "TTT", "TTH", "HTT", "HTH")

# A samc object is the core of the package
samc_obj <- samc(p_mat)

# Given the last 3 flips, how many more flips until we hit HTH (absorption)?
survival(samc_obj)
[1] 420.0000 421.0526 420.0000 421.0526 441.0526 421.0526 441.0526


So, if we start with $$HHT$$, we expect it would take 420 flips on average to end up with $$HTH$$. $$THT$$ is the same, and they represent the best case scenario where we start only one flip away from $$HTH$$. Conversely, $$TTT$$ and $$HTT$$ are a worst case scenario where we start off needing 3 perfect flips; the result for these is 441.0526, the same as other answers.

As Neil Slater pointed out in a comment, the transition matrix can be reduced to 5 elements if this all you want. But there are other metrics to explore. For example, we can calculate the number of times we expect a sequence to occur before we hit $$HTH$$:

$$F = (I-Q)^{-1}$$

Here's some code (I include some code to show the relationship to the previous metric as well):

# Given a starting point (in this case, a sequence of 3 flips), how many times
# would we expect the different combinations of 3 flips to occur before absorption?
visitation(samc_obj, origin = "HHT")
[1]   1.95   0.05  18.05   0.95 361.00  19.00  19.00

sum(visitation(samc_obj, origin = "HHT")) # Compare to survival() result
[1] 420

# Instead of a start point, we can look at an end point and how often we expect
# it to occur for each of the possible starting points
visitation(samc_obj, dest = "THT")
[1] 18.05 18.05 19.05 18.05 19.00 19.00 19.00

# These results are just rows/cols of a larger matrix. We can get the entire matrix
# of the start/end possibilities but first we have to disable some safety measures
# in place because this package is designed to work with extremely large P matrices
# (millions of rows/cols) where these types of results will consume too much RAM and
# crash R
visitation(samc_obj)
Error: This version of the visitation() method produces a large dense matrix.
See the documentation for details.

samc_obj$override <- TRUE visitation(samc_obj) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.95 0.05000000 18.05 0.95 361 19 19 [2,] 1.95 1.10263158 18.05 0.95 361 19 19 [3,] 0.95 0.05000000 19.05 0.95 361 19 19 [4,] 1.95 0.10263158 18.05 1.95 361 19 19 [5,] 1.00 0.05263158 19.00 1.00 381 20 19 [6,] 1.00 0.05263158 19.00 1.00 361 20 19 [7,] 1.00 0.05263158 19.00 1.00 380 20 20 rowSums(visitation(samc_obj)) # equivalent to survival() above [1] 420.0000 421.0526 420.0000 421.0526 441.0526 421.0526 441.0526  Another one: $$D=(F-I)diag(F)^{-1}$$ would calculate the probability of a sequence of flips occurring before you hit $$HTH$$. And there's a "short-term" version of this: $$\tilde{D}_{jt}=(\sum_{n=0}^{t-1}\tilde{Q}^n)\tilde{q}_j$$ which calculates the same thing, but within a given number of time steps (or coin flips in this context). So let's say you're interested in the probability of flipping $$TTT$$ with only 3 flips: dispersal(samc_obj, dest = "TTT", time = 3) [1] 0.902500 0.857375 0.902500 0.857375 NA 0.857375 0.950000  Depending on your last 3 flips, you basically have 3 options (4 if you count having already flipped $$TTT$$): $$0.95$$, $$0.95^2=0.9025$$, or $$0.95^3= 0.857375$$. Obviously, some of these things can be easier or more intuitive to calculate using other methods, but once set up, absorbing markov chains give you a lot of flexibility in exploring a scenario. A full list of the things you can calculate (at least using my package) are available in the function reference: samc functions. Which calculations are relevant depends on the context. When writing the functions, I used the book “Finite Markov Chains” by Kemeny and Snell as reference, which includes proofs. You can find the pdf of it online for free (legally, I believe) pretty easily. Tying things back to OP's questions: as the other answers have shown, there are a variety of ways to mathematically model your experiment, including Markov chains. One thing Stephen Kolassa alluded to, and what I was hoping to show, is there is a LOT more you can learn about this experiment than you might have realized. Note: the package was originally written for spatial applications, so a lot of the terminology in it is biased towards that. However, it is usable for pretty much any application of absorbing Markov chains Seems very close to Shannon - related theorems. If you posit "HTH" as your "end of message" string, you want to estimate the chance that "HTH" shows up in random data. I suspect a little digging into his work will provide the equations/ formulas of interest. And because I can't resist, "HTH" • The answers for HTH and HHT differ by$20$, though you might think they are equally likely to turn up at a given position of a random string Nov 24, 2021 at 16:21 • @Henry but that's because of the relative probabilities in this case. But you raise a good point. To apply Shannon theory, I probably need to create 19 "T-type" messages and one "H-type" message, and then calculate the probabilities for any$HT_jH\$ string Nov 24, 2021 at 17:54

Final formula appears to be

sum(i = 1 to length(pattern) : if (first i flips of pattern match the last i flips of pattern) then P(series of i flips matches first i flips of pattern)^-1 else 0)

(For a fair coin, this reduces to Conway's algorithm; it can probably be proven by a similar method.)