<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Adventures in coin flipping
AKA Introduction to the Binomial Distribution
Step1:  The binomial distribution is closely related to the normal distribution (aka Gaussian distribution)
The probability density for the Gaussian distribution is
Step2:  Could we figure this out analytically?
General Form
$$
p(\text{k heads in n flips}) = (\text{prob. of this many heads}) * \
                               (\text{prob. of this many tails}) * \
                               (\text{how many possible orderings?})
$$
Specifics
$$
\begin{array}{c}
p(\text{k heads in n flips}) \leftarrow {n \choose k} * p^{k} * (1-p)^{(n-k)} \\
p \leftarrow \text{probability of heads in a single flip} \\
p^k \leftarrow \text{total probability of k heads} \\
1-p \leftarrow \text{probabilty of one tails} \\
n-k \leftarrow \text{number of tails} \\
(1-p)^{(n-k)} \leftarrow \text{ probability of all the tails} \\
{n \choose k} \leftarrow \text{ all the possible orderings of k heads in n flips} \\
\end{array}
$$
Reminder
Step3:  How can we use the mean and standard deviation to estimate the probability?
FACT
<ASSISTANT_TASK:>
Python Code:
import random
results = []
for trial in xrange(10000):
    heads = 0
    for i in xrange(100):
        flip = random.randint(0,1)
        if (flip == 0):
            heads += 1
    results.append(heads)
print results[1:10]
import matplotlib.pyplot as plt
plt.figure()
plt.hist(results)
plt.show()
## Plot the histogram using integer values by creating more bins
plt.figure()
plt.hist(results, bins=range(100))
plt.title("Using integer values")
plt.show()
## Plot the density function, notice bars sum to exactly 1 
## Also make the plot bigger
plt.figure(figsize=(15,6))
plt.hist(results, bins=range(100), normed=True)
plt.title("coin flip densities")
plt.show()
Explanation: Adventures in coin flipping
AKA Introduction to the Binomial Distribution
End of explanation
flips_mean = float(sum(results)) / len(results)
print flips_mean
## the numpy package has lots of useful routines: http://www.numpy.org/
import numpy as np
mean = np.mean(results)
print mean
## we could code standard deviation by hand, but numpy makes it easier
stdev=np.std(results)
print stdev
## Overlay a normal distribution on top of the coin flip data
plt.figure(figsize=(15,6))
count, bins, patches = plt.hist(results, bins=range(100), normed=True, label="coin flip histogram")
plt.plot(bins, 1/(stdev * np.sqrt(2 * np.pi)) *
                np.exp( - (bins - mean)**2 / (2 * stdev**2) ),
          linewidth=3, color='red', label="normal distribution")
plt.title("Coin flip densities with normal distribution overlay")
plt.legend()
plt.show()
Explanation: The binomial distribution is closely related to the normal distribution (aka Gaussian distribution)
The probability density for the Gaussian distribution is
End of explanation
prob_heads = .5
num_flips = 100
num_heads = 25
prob_flips = np.math.factorial(num_flips) / \
            (np.math.factorial(num_heads) * np.math.factorial(num_flips-num_heads)) * \
            (prob_heads**num_heads) * ((1-prob_heads)**(num_flips-num_heads))
print "The probability of seeing %d heads in %d flips is %.015f" % (num_heads, num_flips, prob_flips)
## Another super useful package is scipy
import scipy.stats
sp_prob = scipy.stats.binom.pmf(num_heads, num_flips, prob_heads)
print "scipy computed it as %0.15f" % sp_prob
## normal approximatation
print scipy.stats.norm(50, 5).pdf(25)
## Overlay a normal distribution on top of the coin flip data
plt.figure(figsize=(15,6))
count, bins, patches = plt.hist(results, bins=range(100), normed=True, label="coin flip histogram")
plt.plot(bins, scipy.stats.binom.pmf(bins, num_flips, prob_heads),linewidth=3, color='red', label="binomial distribution")
plt.plot(bins, scipy.stats.norm(50,5).pdf(bins),linewidth=3, color='green', linestyle='--', label="normal distribution")
plt.title("Coin flip densities with normal distribution overlay")
plt.legend()
plt.show()
Explanation: Could we figure this out analytically?
General Form
$$
p(\text{k heads in n flips}) = (\text{prob. of this many heads}) * \
                               (\text{prob. of this many tails}) * \
                               (\text{how many possible orderings?})
$$
Specifics
$$
\begin{array}{c}
p(\text{k heads in n flips}) \leftarrow {n \choose k} * p^{k} * (1-p)^{(n-k)} \\
p \leftarrow \text{probability of heads in a single flip} \\
p^k \leftarrow \text{total probability of k heads} \\
1-p \leftarrow \text{probabilty of one tails} \\
n-k \leftarrow \text{number of tails} \\
(1-p)^{(n-k)} \leftarrow \text{ probability of all the tails} \\
{n \choose k} \leftarrow \text{ all the possible orderings of k heads in n flips} \\
\end{array}
$$
Reminder: 
$$
{n \choose k} = \frac{n!}{k!(n-k)!}
$$
End of explanation
expected_mean  = num_flips * prob_heads
expected_stdev = np.math.sqrt(num_flips * prob_heads * (1 - prob_heads)) 
print "In %d flips, with a probability %.02f" % (num_flips, prob_heads)
print "The expected frequency is %.02f +/- %.02f" % (expected_mean, expected_stdev)
print "The observed frequency was %0.2f +/- %0.2f" % (mean, stdev)
Explanation: How can we use the mean and standard deviation to estimate the probability?
FACT: The mean of the binomial distribution is
$$
mean = n * p
$$
FACT: The standard deviation of the binomial distribution is
$$
stdev = \sqrt{np(1-p)}
$$
End of explanation
<END_TASK>