A serious in-depth problem of probabilities for flipping coins Lets say I'm doing 10,000 flips of a coin. I would like to know the probability of how many flips it takes to get 4 or more consecutive heads in a row. 
The count would work as the following, you would count one successive round of flips being only heads (4 heads or more). When a tails hits and breaks the streak of heads the count would start again from the next flip. This would then repeat for 10,000 flips.
I'd like to know the probability of not just 4 or more heads in a row but, 6 or more, and 10 or more. To clarify if a streak of 9 heads is achieved it would be tallied as 1 streak 4 or more (and/or 6 or more), not 2 separate streaks. For example if the coin came THTHTHTHHHHHH///THAHTHT....the count would be 13 and begin once again on the next tails. 
Let's say the data comes out to be highly skewed to the right; the mean being 40 flips on average it takes to achieve a streak of 4 or more, and the distribution is u = 28.  Obviously skewed. 
I'm doing my best to find a way to make some sense out of descriptive data, except as of now I've found nothing. 
I want to find some way of getting some sensible probability out of it. Like a normal curve where +/- 1 SD is 68% etc..I've looked into log normalizing and that's only really used for a parametric testing which is not my goal. 
I have been told beta distributions but every suggestion I've had has been quite confusing. I've asked this question a year ago and got some insight but unfortunately I still do not have an answer. Thank you to any of you who have ideas. 
 A: I am not sure the beta will be likely to be particularly suitable as a way of dealing with this problem -- "The number of plays until ..." is clearly a count. It's an integer, and there is no upper limit on values where you get positive probability. 
By contrast the beta distribution is continuous, and on a bounded interval, so it would seem to be an unusual choice. If you moment match a scaled beta, the cumulative distribution functions might perhaps approximate not so badly in the central body of the distribution. However, some other choice is likely to substantially better further into either tail. 
If you have either an expression for the probabilities or simulations from the distribution (which you presumably need in order to find an approximating beta), why would you not use those directly?

If your interest is in finding expressions for probabilities or the probability distribution of the number of tosses required, probably the simplest idea is to work with probability generating functions. These are useful for deriving functions from recursive relationships among probabilities, which functions (the pgf) in turn allow us to extract whatever probabilities we require.
Here's a post with a good answer taking the algebraic approach, which explains both the difficulties and makes good use of pgfs and recurrence relations. It has specific expressions for mean and variance in the "two successes in a row" case:
https://math.stackexchange.com/questions/73758/probability-of-n-successes-in-a-row-at-the-k-th-bernoulli-trial-geometric
The four successes case will be substantially more difficult of course. On the other hand, $p = \frac{1}{2}$ does simplify things somewhat.
-- 
If you just want numerical answers, simulation is relatively straightforward. The probability estimates could be used directly, or alternatively it would be reasonable to smooth the simulated probabilities.
If you must use an approximating distribution, you can probably choose something that does pretty well. 
It's possible a mixture of negative binomials (the 'number of trials' version rather than 'the number of successes' version) might be reasonable. Two or three components should be expected to give a good approximation in all but the extreme tail. 
If you want a single continuous distribution for an approximation, there may be better choices than the beta distribution; it would be something to investigate.

Okay, I've since done a little algebra, some playing with recurrence relations, some simulation and even a little thinking.
To a very good approximation, I think you can get away with simply specifying the first four nonzero probabilities (which is easy), computing the next few handfuls of values via recurrence (also easy) and then using a geometric tail once the recurrence relation has smoothed out the initially less smooth progression of probabilities.
It looks like you can use the geometric tail to very high accuracy past k=20, though if you're only worried about say 4 figure accuracy you could bring it in earlier.
This should let you compute the pdf and cdf to good accuracy. 
I'm a little concerned - my calculations give that the mean number of tosses is 30.0, and the standard deviation is 27.1; if I understood what you mean by "x" and "u", you got 40 and 28 in your tossing. The 28 looks okay but the 40 seems quite a way off what I got... which makes me worry I did something wrong.
====
NOTE: Given the complexities between first time and subsequent times that we encountered, I just want to be absolutely certain now that we are counting the same thing.
Here's a short sequence, with the ends of the '4 or more H' sequences marked (pointing to the gap between flips immediately after the last H)
       \/                     \/
TTHHHHHHTTHTTTTTHHTTHTTHHTHHHHHT...
       /\                     /\

Between those two marks I count 23 flips; that is as soon as the previous sequence of (6 in this case) H's ends, we start counting at the immediately following T and then we count right to the end of the sequence of 5 H's (in this case) that ends the next sequence, giving a count of 23 in this case.
Is that how you count them?

Given the above is correct, this is what the probability function of the number of tosses after one run of at least 4 heads is complete until the next run of at least 4 heads is complete looks like:

At first glance it looks like it's flat for the first few values, then has a geometric tail, but that impression is not quite accurate - it takes a while to settle down to an effectively geometric tail.
I am working on coming up with a suitable approximation you can use to answer whatever questions you'd have about probabilities associated with this process to good accuracy that is at the same time as simple as possible. I have a pretty good approximation that should work (that I have already checked against a simulation of a billion coin tosses) but there's some (small but consistent) bias in the probabilities the approximation gives in part of the range and I'd like to see if I can get an extra digit of accuracy out of it. 
It may be that the best way to do it is simply to give you a table of the probability function and cdf out to a point beyond which a geometric distribution can be used.
However, it would help if you can give some idea of the range of things you need to use the approximation for.

I hope to follow through on the pgf approach, but it's possible someone else will be more proficient with them than me and can do not just the 4-case but other cases.
A: You want the geometric distribution. From Wikipedia:

The probability distribution of the number $X$ of Bernoulli trials
  needed to get one success, supported on the set { 1, 2, 3, ...}.

Let heads H be a failure and tails T be a success. The random variable $X$ is then the number of coin flips needed to see the first tails. For example, $X=4$ will be the sequence HHHT. Here is the probability distribution for $X$:
$$
P(X=x)=(1-p)^{x-1}p
$$
However, we only want the number of heads. Let's instead define $Y=X-1$ as the number of heads. Here is it's distribution:
$$
P(Y+1=x) = (1-p)^{x-1}p \\
P(Y=x-1) = (1-p)^{x-1}p \\
P(Y=y) = (1-p)^yp
$$
for $y=0,1,2,3...$. We assume a fair coin, making $p=0.5$. Therefore:
$$
\begin{align}
P(Y=y) &= (0.5)^y(0.5) \\
&= 0.5^{y+1}
\end{align}
$$
This all assumes the number of flips $n$ is sufficiently large (like 10,000). For smaller (finite) $n$, we would need to add a normalization factor to the expression. Put simply, we need to ensure that the total sum is equal to 1. We can do this by dividing by the sum of all probabilities, defined here as $\alpha$:
$$
\alpha = \sum_{i=0}^{n-1}{P(Y=i)}
$$
This means the corrected form of $Y$, denoted $Z$, will be:
$$
\begin{align}
P(Z=z) &= \frac{1}{\alpha}(1-p)^zp \\
&= \frac{1}{\sum_{i=0}^{n-1}{(1-p)^ip}}(1-p)^zp
\end{align}
$$
Again, with $p=0.5$, we can reduce this even further, using a geomoetric series summation:
$$
\begin{align}
P(Z=z) &= \frac{1}{\sum_{i=0}^{n-1}{0.5^{i+1}}} 0.5^{z+1} \\
&= \frac{1}{1-0.5^n} 0.5^{z+1} \\
&= \frac{0.5^{z+1}}{1-0.5^n}
\end{align}
$$
And we can see that, as $n \rightarrow \infty$, our modified version $Z$ approaches $Y$ from earlier.
