<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A Practical Application of Counting & Likelihood Estimation
by me, Fiona Pigott
Problem
Step1:  Getting a feel for it
Step2:  How likely is it that we get exactly X% of $tau_1$?
$$ P(t_1 = s\tau_1)  = {\tau_1\choose s\tau_1} s^{t_1}(1-s)^{(\tau_1 - s\tau_1)}$$
Step3:  How likely is it that we get $t_1$ samples, with $t_1$ on some interval? Say, $s\tau_1 \pm$ 1%?
$$ P(.99s\tau_1 < t_1 < 1.01s\tau_1)   = \sum_{i = .99s\tau_1}^{1.01s\tau_1}{\tau_1\choose i} s^{i}(1-s)^{(\tau_1 - i)} $$
Step4:  Standard deviation of the size of the sample
It's relatively simple to work out the standard deviation of the Binomial distribution, if we remember a few things about expected values (averages), standard deviation and variance
Step5:  Using the standard deviation to estimate the most likely interval for $t_1$ to fall in
Step6:  Now let's look at $t_1$ and $t_2$ together
Step7:  Finally, the joint PMF
Step8:  That was interesting, but it wasn't exactly the problem that we wanted to answer
What we just did assumed that we knew $\tau_1$ and $\tau_2$, and we figured out probabilities of getting $t_1$ and $t_2$ given those specific $\tau$ values. Actually, $\tau$ is the value that we don't know (the value we're trying to estimate).
We want to know
Step9:  I've heard people say the words "log-likelihood"
Me too! Notice how the values of the likelihood (y-axis) are very, very small? Taking the log of both sides of the likelihood equation lets us maximize something a little computationally friendlier than 10^-20 or whatever, and since a logarithm is a monotonic function, it doesn't change the position of the max. Remember, the log of a product is a sum of the logs.
$$ L(\theta) = \prod_{i = 0}^n P(X_i= x_i; s = \theta)$$
$$ \log(L(\theta)) = \sum_{i = 0}^n \log(P(X_i= x_i; s = \theta))$$
Step10:  And finally, there's a way to do the above analytically (rather than numerically)
Recall
Step11:  3. Estimating the value of $\tau$
Now, remember the Binomial distribution again for a second
Step12:  Finally
<ASSISTANT_TASK:>
Python Code:
# math
import numpy as np
from scipy.misc import comb
from scipy.stats import binom, norm
from itertools import accumulate
from math import sqrt, ceil, log
# python tools
from collections import Counter
from itertools import zip_longest
from multiprocessing import Process, Pool, Queue
from random import random
import os
# plotting
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: A Practical Application of Counting & Likelihood Estimation
by me, Fiona Pigott
Problem: Find the most Retweeted Tweet
Given: set of Retweets of various original Tweets
We want to use a collected sample of all Retweets to estimate which Tweet was Retweeted the most. We can count up all of the Retweets in our sample to estimate the ordering of the Tweets, but we can't count every Retweet to get the exact number. How big of a sample do we need to be relatively certain that we've correctly identified the most Retweeted Tweet?
Tools that we have:
PowerTrack with the sampling operator
PowerTrack sampling works by flipping a biased coin for each Tweet to decide whether or not to include it in the sample: if the sample is 10%, each Tweet has an independent 10% chance of being included in the sample.
We can select only Retweets using a PowerTrack operator
More information on PowerTrack: http://support.gnip.com/apis/powertrack2.0/rules.html
What we need to know:
What's the probability that, for some sampling factor (using PowerTrack sampling) I get this ordering wrong (i.e., I mistakenly identify some Tweet as the most Retweeted Tweet).
Problem statement:
Let the Tweet most Retweeted in the set have a total of $\tau_1$ Retweets
Let the 2nd most Retweeted in the set have a total of $\tau_2$ Retweets
Take a sample of the whole set of Retweets, including each individual Retweet with a probability $s$ and excluding it with a probability $(1-s)$. Note that the inclusion of any given Retweet is independent of whether or not any other Retweet is included.
Call the number of occurrences of Retweets of the most Retweeted and 2nd most Retweeted Tweets in the sample $T_1$ and $T_2$ respectively, and the 
Some probability vocabulary:
A random variable: "A variable quantity whose value depends on possible outcomes." In this case, $T_1$ and $T_2$ are "random variables"
An event: A measurable outcome of the experiment (i.e., "$T_1 = t_1$," where $t_1$ is some specific number of Retweets in the sample)
Independent events: Two (or more) events whose outcomes do not affect eachother. In this example, all of our events are pretty much independent of eachother.
A PMF or probabilty mass function: For a discrete probabilty distribution (in this case, we're dealing with a discrete distribution, because we cannot collect a sample of 10.5 Tweets), the probability mass function is a function that gives the probability of each possible event, and the sum over all possible events is 1. 
A joint PMF: A function that gives the probabilty of two events occurring together. For two independent events, we can just multiply the two individual PMFs together. Easy.
And now! Counting.
Counting probabilities of $t_1$ and $t_2$ occurences:
The probability mass function of $T_1$ is: 
The probabilty of selecting $t_1$ Tweets, so $s^{t_1}$ times the probabilty of not selecting $\tau_1 - t_1$ Tweets, times the number of different ways you could select $t_1$ Tweets, ${\tau_1\choose t_1}$:
$$ P(T_1 = t_1) = {\tau_1\choose t_1} s^{t_1}(1-s)^{(\tau_1 - t_1)}$$
Similarly, the probability mass function of $T_2$ is:
$$ P(T_2 = t_2) = {\tau_2\choose t_2} s^{t_2}(1-s)^{(\tau_2 - t_2)}$$
This distribution might look familiar: it's called the binomial distribution
"In probability theory and statistics, the binomial distribution with parameters n and p is the discrete probability distribution of the number of successes in a sequence of n independent experiments, each asking a yes–no question, and each with its own boolean-valued outcome" 
from Wikipedia
 $$ B(n,p) =    \textstyle {n \choose k}\,p^{k}(1-p)^{n-k} $$
The probability of exactly $t_1$ and exactly $t_2$ occurrences (joint PMF):
$$ {\tau_1\choose t_1} s^{t_1}(1-s)^{(\tau_1 - t_1)} {\tau_2\choose t_2} s^{t_2}(1-s)^{(\tau_2 - t_2)} = $$
$$ {\tau_1\choose t_1}{\tau_2\choose t_2} s^{t_1 + t_2} (1-s)^{\tau_1 + \tau_2 - t_1 - t_2} $$
End of explanation
# Set up the problem. We have to fix tau_1, tau_2, and s
tau1 = 100
tau2 = 90
s = .1
# We can use the binomial PMF from scipy, which is better than coding it up myself 
# because all of those factorials make for hard computations
# we're gonna use these a lot
binomial_tau1 = binom(tau1,s)
binomial_tau2 = binom(tau2,s)
# Let's show how likely it is to get exactly t_1 of T_1. Just for practice
# do the probability calculations
probabilities_t1 = []
for i in range(0,tau1):
    #probabilities_t1.append((i,comb(tau1,i)*s**(i)*(1-s)**(tau1-i)))
    probabilities_t1.append((i,binomial_tau1.pmf(i)))
probabilities_t2 = []
for i in range(0,tau2):
    #probabilities_t2.append((i,comb(tau2,i)*s**(i)*(1-s)**(tau2-i)))
    probabilities_t2.append((i,binomial_tau2.pmf(i)))
plt.plot([x[0] for x in probabilities_t1],[x[1] for x in probabilities_t1])
plt.title("Probability Mass Function for the Binomial distribution")
plt.xlabel("Number of successes in tau_1 trials (t_1)")
plt.ylabel("Probability of exactly t_1 successes")
Explanation: Getting a feel for it:  What's the probability that $t_2 \geq t_1$ given that $\tau_1 > \tau_2$?
$$ \sum_{j = 0}^{\tau_2}\sum_{i = j}^{\tau_2} {\tau_1\choose j}{\tau_2\choose i} s^{j + i} (1-s)^{\tau_1 + \tau_2 - j - i}$$
Note: $t_1 \rightarrow j$, $t_2 \rightarrow i$
Before we start with the joint probabilty, let's just look at how $t_1$ is distributed:
End of explanation
print("The probability of getting exactly s*tau_1 samples is: {:f}".format(binomial_tau1.pmf(s*tau1)))
Explanation: How likely is it that we get exactly X% of $tau_1$?
$$ P(t_1 = s\tau_1)  = {\tau_1\choose s\tau_1} s^{t_1}(1-s)^{(\tau_1 - s\tau_1)}$$
End of explanation
plus_or_minus_percent = .01
lower_bound_t1 = int((1 - plus_or_minus_percent)*s*tau1)
upper_bound_t1 = int((1 + plus_or_minus_percent)*s*tau1)
probability_t1_interval = 0
for i in range(lower_bound_t1,upper_bound_t1):
    probability_t1_interval += binomial_tau1.pmf(i)
print("P({} < t_1 < {}) = {:f}".format(lower_bound_t1,upper_bound_t1,probability_t1_interval))
Explanation: How likely is it that we get $t_1$ samples, with $t_1$ on some interval? Say, $s\tau_1 \pm$ 1%?
$$ P(.99s\tau_1 < t_1 < 1.01s\tau_1)   = \sum_{i = .99s\tau_1}^{1.01s\tau_1}{\tau_1\choose i} s^{i}(1-s)^{(\tau_1 - i)} $$
End of explanation
# standard deviation varies as the size of tau1 varies
std_dev_sizes = []
max_tau = 10**4
for i in range(0,max_tau,5):
    std_dev_sizes.append((i,sqrt(i*s*(1-s))))
plt.plot([x[0] for x in std_dev_sizes],[x[1] for x in std_dev_sizes])
_ = plt.xlim([0,max_tau])
_ = plt.ylim([0,max([x[1] for x in std_dev_sizes])])
_ = plt.title("Standard deviation of the sample size, varying tau")
Explanation: Standard deviation of the size of the sample
It's relatively simple to work out the standard deviation of the Binomial distribution, if we remember a few things about expected values (averages), standard deviation and variance:
1. Expected value is a linear operator (i.e., $E[A + 2B] = E[A] + 2E[B]$)
2. We'll call $E[X] = \mu$ 
3. Variance: $Var[X] = E[(X - \mu^2)] = E[X^2] - E[X]^2$
4. If two random variables are uncorrelated, Variance is a linear operator: $Var[\sum_i X_i] = \sum_i Var[X_i]$
5. Standard deviation: $\sqrt{Var[X]}$  
And remember one thing about the Binomial distribution:
"In probability theory and statistics, the binomial distribution with parameters n and p is the discrete probability distribution of the number of successes in a sequence of n independent experiments (Bernoulli trials), each asking a yes–no question, and each with its own boolean-valued outcome" 
from Wikipedia
I'm not going to work this out in a lot of detail here, I'll simply say that the expected value ($E[X]$) of a single Bernoulli trial (i.e., for one Tweet, decide whether it is in the sample or outside of the sample) is equal to $s$, and the variance of that single trial is equal to $Var[X] = s(1-s)$.  
Then, if the total number of Tweets in the sample is simply the sum of these independent Bernoulli trials, the expected value of $t_1$ is $E[t_1] = s\tau_1$ and the variance is $Var[t_1] = \tau_1s(1-s)$
The standard deviation is then: $\sqrt{\tau_1s(1-s)} $
End of explanation
# Let's show how likely it is to get exactly t_1 of T_1. Just for practice
# we cen grab the 
normal_distribution_tau1 = norm(tau1*s,sqrt(tau1*s*(1-s)))
normal_prob_t1 = []
probabilities_t1 = []
for i in range(0,tau1):
    probabilities_t1.append((i,binomial_tau1.pmf(i)))
    normal_prob_t1.append((i,normal_distribution_tau1.pdf(i)))
plt.plot([x[0] for x in probabilities_t1],[x[1] for x in probabilities_t1],linewidth = 5)
plt.plot([x[0] for x in normal_prob_t1],[x[1] for x in normal_prob_t1],"--",linewidth = 3)
plt.legend(["binomial","normal"])
plt.ylabel("Probabilty mass a x")
plt.xlabel("t_1")
# the total probability of falling within 3 standard deviations of the mean (s*tau_1):
plus_or_minus_stdv = .01
lower_bound_t1 = max(0,int(s*tau1 - 3 * sqrt(tau1*s*(1-s))))
upper_bound_t1 = int(s*tau1 + 3 * sqrt(tau1*s*(1-s)))
probability_t1_interval = 0
for i in range(lower_bound_t1,upper_bound_t1):
    probability_t1_interval += binomial_tau1.pmf(i)
print("P({} < t_1 < {}) = {:f} (very nearly 99.7%)".format(lower_bound_t1,upper_bound_t1,probability_t1_interval))
Explanation: Using the standard deviation to estimate the most likely interval for $t_1$ to fall in:
Now, I'm not going to work this out explicitly because there are other fish I'd like to fry, but know that as the number of trials increases, the Binomial distribution tends towards the normal distribution.
That means we can use the same rule of thumb to estime how much of the distribution falls within, say, 3 standard deviations of the mean that we use for normal distributions (that is, 99.7% of all of the probabilty mass falls within 3 standard deviations of the mean).
I'm going to use a plot to try to get you to beleive me:
End of explanation
# t1 and t2 most likely fall within some range
# t1 and t2 are selected form binomial distributions, they are symmetric 
# and their standard deviations are easily calculated
fig, ax = plt.subplots(1, 1)
ax.plot([x[0] for x in probabilities_t1],[x[1] for x in probabilities_t1])
ax.plot([x[0] for x in probabilities_t2],[x[1] for x in probabilities_t2])
ax.legend(["t_1","t_2"])
ax.set_ylabel("Probabilty mass at x")
Explanation: Now let's look at $t_1$ and $t_2$ together:
End of explanation
# Now let's code this up so we can simulate it for various values of s, tau_1, and tau_2
probabilities_t2_and_t1 = np.zeros((tau1,tau2))
# i can seriously speed up the amount of time this takes by not calculating the very, very low probabilities
# outside of 3 stdev from the mean
bool_truncate_window = True
if bool_truncate_window:
    stdv1 = ceil(sqrt(tau1*s*(1-s)))
    stdv2 = ceil(sqrt(tau2*s*(1-s)))
    num_stdv = 3
    window_t1 = [max(ceil(s*tau1)-num_stdv*stdv1,0), ceil(s*tau1)+num_stdv*stdv1]
    window_t2 = [max(ceil(s*tau2)-num_stdv*stdv2,0), ceil(s*tau2)+num_stdv*stdv2]
else:
    window_t1 = [0,tau1]
    window_t2 = [0,tau2]
for j in range(window_t1[0],window_t1[1]):
    for i in range(window_t2[0],window_t2[1]):
        probabilities_t2_and_t1[j,i] = binomial_tau1.pmf(j)*binomial_tau2.pmf(i)
figure, ax = plt.subplots(1,1,figsize = (6,12))
ax.imshow(probabilities_t2_and_t1, cmap='Blues', interpolation='nearest', aspect = 'equal', origin = 'lower')
ax.set_xlim([0,tau2])
ax.set_ylim([0,tau1])
if bool_truncate_window:
    ax.fill_between(window_t2,0,tau1,alpha = .2,color = 'navajowhite')
    ax.fill_between([0,tau2],window_t1[0],window_t1[1],alpha = .2,color = 'lightsteelblue')
ax.plot([0,tau2],[0,tau2],'--', color = 'gray')
# The total probabilties should sum to 1
# because I truncated, they sum to something slightly less than one--you can see how small the actual 
# values are in those other areas 
print("Total probabilty represented: {:f}".format(probabilities_t2_and_t1.sum()))
# Now, sum only the lower triangle--the piece where t_2 >= t_1
print("The probabilty that t_2 >= t_1 given tau1 = {} and tau2 = {}: {}".format(
    tau1, tau2, np.triu(probabilities_t2_and_t1).sum()))
Explanation: Finally, the joint PMF:
$$ {\tau_1\choose t_1} s^{t_1}(1-s)^{(\tau_1 - t_1)} {\tau_2\choose t_2} s^{t_2}(1-s)^{(\tau_2 - t_2)}$$
What's the probability that $t_2 \geq t_1$ given that $\tau_1 > \tau_2$?
$$ \sum_{j = 0}^{\tau_2}\sum_{i = j}^{\tau_2} {\tau_1\choose j}{\tau_2\choose i} s^{j + i} (1-s)^{\tau_1 + \tau_2 - j - i}$$
Note: $t_1 \rightarrow j$, $t_2 \rightarrow i$
End of explanation
def bernoulli_experiment(s):
    if random() < s:
        return True
    else:
        return False
def bernoulli_pmf(x, s):
    if x:
        return s
    else:
        return 1-s
def eval_bernoulli_likelihood_for_all_x_i(x_i, s):
    likelihood = 1
    for x in x_i:
        likelihood = likelihood * bernoulli_pmf(x,s)
    return likelihood
# set up the experiment
# number of trials
n = 100
# probabilty, s. let's assign it randomly and check later
s_sample = random()
s_range = list(np.linspace(0,1,100))
# results, x_i of n experiments
x_i = [bernoulli_experiment(s_sample) for x in range(0,n)]
# get the likelihood results
likelihood = [eval_bernoulli_likelihood_for_all_x_i(x_i, s_) for s_ in s_range]
estimated_s = max(list(zip(s_range,likelihood)), key = lambda x:x[1])
# plot
plt.plot(s_range,likelihood)
plt.title("Max Likelihood as a function of s")
plt.legend(["actual value of s: s = {:f}".format(s_sample) + "\n" + "estimated value of s: {:f}".format(estimated_s[0])])
plt.xlabel("s")
plt.ylabel("likelihood")
Explanation: That was interesting, but it wasn't exactly the problem that we wanted to answer
What we just did assumed that we knew $\tau_1$ and $\tau_2$, and we figured out probabilities of getting $t_1$ and $t_2$ given those specific $\tau$ values. Actually, $\tau$ is the value that we don't know (the value we're trying to estimate).
We want to know: What's the probability that $\tau_1$ > $\tau_2$ given that $t_1 > t_2$? What's the probability that we get our ordering (of the most popular, second most popular, etc Retweets) correct?
Basically, we know that the number of Retweets of any given Tweet in our sample is going to follow this Binomial distribution, where $\tau$ is the number of Retweets in total (the number of Retweets that we would have collected with a 100% sample) and $t$ is:
$$ P(T = t) = {\tau\choose t} s^{t}(1-s)^{(\tau - t)}$$
We know that, because I've already typed it several times in this notebook. Now, when we actually run this experiment, we know the probability of including any individual Retweet in the sample, that probabilty is $s$, and we set it. We also know the outcome of each experiment--that outcome is the number of Retweets that we actually count in our sample. The parameter of the distribution that we don't know is $\tau$--the total number of Retweets.
Estimating the value of a parameter in a probabilty distribution
One often-used and often-referenced way to find the most likely value of a parameter is called "Maximum Likelihood Estimation," and I think that right now would be a good time to talk about Maximum Likelihood Estimators (because it's a generally useful and interesting concept).
Now, it turns out that finding a nice, tidy Maximum Likelihood Estimator for $\tau$ in this distribution is very difficult (you'll see why soon enough), so I'm going to do two things:
1. Explain, generally, what MLEs are
2. Pretend that I know $\tau$ and don't know $s$, and demonstrate a simple, clean example of using MLEs to estimate $s$ in this Binomial distribution
3. Use numerical Maximum Likelihood Estimation to estimate $\tau$ (closer to the actual problem I'm trying to solve)
1. What is Maximum Likelihood estimation?
Maximum Likelihood Estimation uses the outcomes of many iterations of an experiment to guess (estimate) the paramater that maximizes the probabilty (likelihood) of those outcomes. That's a pretty intuitive idea.
For example: If you flip a single coin 100 times, and that coin comes up heads 30 times, you might guess that the probabilty of that coin coming up heads in general is 0.3. That would be a pretty good guess, because the probabilty of getting 30/100 "heads" results is maximized if the probabilty of the coins coming up "heads" is 0.3.
Let's work that out:
For any set of independent, identically distributed experiments ($X_i$) and outcomes ($X_i = x_i$), the probabiliy of all of the outcomes occurring at once would simply be the joint PMF of all of the experiments (remember, because they're independent, we can just multiply PMFs together). 
Call each PMF $P(X_i = x_i; \theta)$ (where $"; \theta"$ denotes the value of the parameter that we're trying to estimate, which the PMF is dependent on) then the probability (likelihood) of a set of $n$ results is dependent on a independent variable $\theta$:
$$ L(\theta) = P(X_1 = x_1, X_2 = x_2, X_3 = x_3, \dots, X_n = x_n; \theta) = \prod_{i = 0}^n P(X_i= x_i; \theta)$$
Now, if we find the value of the parameter $\theta$ that maximizes the likelihood of our results $x_i$ coming from the distribution, we have found the MAximum Likelihood Estimator for $\theta$. Super.
2. A simple example of MLE estimation
Right now, we're going to do an example of a Bernoulli distribution (basically, a Binomial distribution with only one trial). Think of the coin toss example again.
End of explanation
def eval_bernoulli_loglikelihood_for_all_x_i(x_i, p):
    likelihood = 0
    for x in x_i:
        likelihood += np.log(bernoulli_pmf(x,p))
    return likelihood
# get the likelihood results
loglikelihood = [eval_bernoulli_loglikelihood_for_all_x_i(x_i, s_) for s_ in s_range[1:-1]]
log_estimated_s = max(list(zip(s_range,likelihood)), key = lambda x:x[1])
# plot
# notice that the max is in the same place
plt.plot(s_range[1:-1],loglikelihood)
plt.title("Max Log Likelihood as a function of s")
plt.legend(["actual value of s: s = {:f}".format(s_sample) + "\n" + "estimated value of s: {:f}".format(
    log_estimated_s[0])])
plt.xlabel("s")
plt.ylabel("likelihood")
Explanation: I've heard people say the words "log-likelihood"
Me too! Notice how the values of the likelihood (y-axis) are very, very small? Taking the log of both sides of the likelihood equation lets us maximize something a little computationally friendlier than 10^-20 or whatever, and since a logarithm is a monotonic function, it doesn't change the position of the max. Remember, the log of a product is a sum of the logs.
$$ L(\theta) = \prod_{i = 0}^n P(X_i= x_i; s = \theta)$$
$$ \log(L(\theta)) = \sum_{i = 0}^n \log(P(X_i= x_i; s = \theta))$$
End of explanation
# let's see what we would estimate s to be analytically
print("Estimate of s is: {} \n (The correct value is {})".format(sum(x_i)/n,s_sample))
Explanation: And finally, there's a way to do the above analytically (rather than numerically)
Recall: the maximum of a function occurs when the derivative is equal to zero. 
Also: It tends to be easier to differentiate the log-likelihood (rather than the likelihood). So let's do that.
Write the Bernoulli PMF as: $ s^{x_i}(1−s)^{1−x_i} $, where $x_i$ has one of two values: 0 or 1.
$$ \log(L(\theta)) = \sum_{i = 0}^n \log(P(X_i= x_i; s = \theta))$$ 
$$ \log(L(\theta)) = \sum_{i = 0}^n \log(\theta^{x_i}(1−\theta)^{1−x_i}) $$
$$ \log(L(\theta)) = (\sum_{i = 0}^n x_i)\log(\theta)+(n−\sum_{i = 0}^n x_i)\log(1−\theta) $$ 
$$ \frac{d\log(L(\theta)}{d\theta} = \frac{d}{d\theta} (\sum_{i = 0}^n x_i)\log(\theta)+\frac{d}{d\theta}(n−\sum_{i = 0}^n x_i)\log(1−\theta) $$
$$ \frac{d\log(L(\theta)}{d\theta} = \frac{\sum_{i = 0}^n x_i}{\theta}+\frac{n−\sum_{i = 0}^n x_i}{1−\theta} $$ 
Find where the derivative is equal to 0:
$$ \frac{d\log(L(\theta)}{d\theta} = \frac{\sum_{i = 0}^n x_i}{\theta}+\frac{n−\sum_{i = 0}^n x_i}{1−\theta} = 0 \text { for } \theta = \frac{\sum_{i = 0}^n x_i}{n}$$
End of explanation
# set up the experiment
# number of trials
n = 1000
# probabilty, s. i'll use the same "s" that we're using for the problem setup
print("s = {}".format(s))
# now, fix a "real" tau
# tau could really be anywhere from 0 to infinity, but I'll fix it between 0 and 100
tau_sample = ceil(random()*100)
# possible tau values
possible_tau_range = range(1,101)
print("The real tau is {}".format(tau_sample))
# the result of n experiments
binomial_tau_sample = binom(tau_sample, s)
t_i = [binomial_tau_sample.rvs() for x in range(0,n)]
#print("The results of the esperiments are: {}".format(t_i))
# Now, evaluate the binomial distribution likelihood for this collection of results
def binomial_loglikelihood(t_i, s, tau):
    binomial_distribution = binom(tau,s)
    loglikelihood = 0
    for t in t_i:
        loglikelihood += np.log(binomial_distribution.pmf(t))
    return loglikelihood
        
# there's no sense in testing tau that are lower than the max of t_i
likelihood_tau = [(x, binomial_loglikelihood(t_i, s, x)) for x in possible_tau_range if x > max(t_i)] 
# estimate tau
estimated_tau = max(likelihood_tau, key = lambda x:x[1])
plt.plot([x[0] for x in likelihood_tau], [x[1] for x in likelihood_tau])
plt.title("Max Likelihood as a function of tau")
plt.legend(["actual value of tau: tau = {:f}".format(tau_sample) 
            + "\n" + "estimated value of tau: {:f}".format(estimated_tau[0])])
plt.xlabel("tau")
plt.ylabel("likelihood")
Explanation: 3. Estimating the value of $\tau$
Now, remember the Binomial distribution again for a second: $$ P(T = t) = {\tau\choose t} s^{t}(1-s)^{(\tau - t)}$$
If we wanted to find the MLE for $\tau$ analytically, we'd have to take the derivative with respect to $\tau$, which would be a bloodbath, pretty much. Just think of all of those factorials. But remember that we were able to numerically estimate the value of $s$ earlier, simply by finding where the PMF is maximized. We can do that to estimate the value of $\tau$.
End of explanation
%%time
# Now let's code this up so we can simulate it for various values of s, tau_1, and tau_2
# fix t1 and t2--the values that we observed IRL
t1 = 100
t2 = 95
# we aren't going to evaluate this over every possible tau, because that would be wayyyyy to large
# let's estimate a window
stdv1 = ceil(sqrt(t1*(1-s)))
stdv2 = ceil(sqrt(t2*(1-s)))
num_stdv = 3
window_est_tau1 = [int(max(ceil(t1)-num_stdv*stdv1,0)*(1/s)), int((ceil(t1)+num_stdv*stdv1)*(1/s))]
window_est_tau2 = [int(max(ceil(t2)-num_stdv*stdv2,0)*(1/s)), int((ceil(t2)+num_stdv*stdv2)*(1/s))]
# calculate the probabilty that binomial distributions with tau1 = j and tau2 = i resulted in samples t1 and t2
# use mutiprocessing bc this is kinda slow
# function that we are evaluating
def evaluate_binomial(j_tau1_i_tau2):
    binomial_tau1 = binom(j_tau1_i_tau2[0],s)
    binomial_tau2 = binom(j_tau1_i_tau2[1],s)
    return (j_tau1_i_tau2, binomial_tau1.pmf(t1)*binomial_tau2.pmf(t2))
# create the processes
pool = Pool(processes=12)
tau2_and_tau1_list = []
probabilities_tau2_and_tau1_list = []
for j_tau1 in range(window_est_tau1[0],window_est_tau1[1]):
    for i_tau2 in range(window_est_tau2[0],window_est_tau2[1]):
        tau2_and_tau1_list.append((j_tau1,i_tau2))
# get the results
probabilities_tau2_and_tau1_list = pool.map(evaluate_binomial, tau2_and_tau1_list)
probabilities_tau2_and_tau1 = np.zeros((window_est_tau1[1],window_est_tau2[1]))
for probabilty_calcd in probabilities_tau2_and_tau1_list: 
     probabilities_tau2_and_tau1[probabilty_calcd[0]] = probabilty_calcd[1]
figure, ax = plt.subplots(1,1,figsize = (6,12))
plt.imshow(probabilities_tau2_and_tau1, cmap='Blues', interpolation='nearest', aspect = 'equal', origin = 'lower')
plt.plot([0,window_est_tau2[1]],[0,window_est_tau2[1]],'--',color = 'gray')
ax.fill_between([0,window_est_tau2[1]],0,window_est_tau1[0],alpha = 1,color = 'white')
ax.fill_between([0,window_est_tau2[0]],window_est_tau1[0],window_est_tau1[1],alpha = 1,color = 'white')
_ = plt.xlim([0,window_est_tau2[1]])
_ = plt.ylim([0,window_est_tau1[1]])
print("Estimate of the probabilty that tau_1 < tau_2 given that t_1 = {} and t_2 = {}: {}".format(
    t1,t2,np.triu(probabilities_tau2_and_tau1).sum()/probabilities_tau2_and_tau1.sum()))
Explanation: Finally: the question that we actually wanted to answer
Given $t_1$ and $t_2$ ($t_1 > t_2$), what is the probabilty that $\tau_1 > \tau_2$?
Now, the more samples we have (the more experiments we run) the better we can estimate $\tau$ . Unfortunately, in real life, using PowerTrack, we're not going to be able to take 100 different random samples of the total, since PowerTrack sampling is deterministic, but estimating the value of the parameter follows the same idea--find the value of the parameter that is most likely to give the observed result.
Now, here's where I'm venturing off the beaten path a little. I'm not that interested in making sure I have the best estimate of $\tau$, as long as I get the ordering ("is $\tau_1$ bigger than $\tau_2$?") correct.
I figure that the proportion of the likelihoods (likelihood that $\tau_1$ is greater than $\tau_2$, divided by the total likelihoods) is probably a good estimate of how likely it is that $\tau_1$ is actually greater than $\tau_2$.
Also note that I can't possibly test the likelihood of every possible value of $\tau$, since $\tau$ could theoretically go to infinity. Instead, I'm making an initial guess at $\tau$ (simply the sample size divided by the sample proportion), and only testing $\tau$ values within 3 standard deviations of that value. I'm not sure that's a particularly good way to do it, since the standard deviation is really about the deviation of the sample size, not the deviation of $\tau$. But it seems okay.
End of explanation
<END_TASK>