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!

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} 

