A handy trick to count uncountable things — from Instagram clout, to penguin population and Man City fans in North London

We have all been there …

With hopeful eyes, we framed the Christmas spirit through our iPhone 15 lens, capturing the dinner table with an artist’s persistence and a poet’s soul, each photo a festive vignette brimming with friends, the shimmering glaze of the turkey, and the twinkle of ornaments. We meticulously pair every snapshot with an upbeat caption, carefully woven from the threads of our thoughts, hoping each post would jingle all the way through Instagram’s bustling holiday traffic.

And yet, for all the hard work to unveil our holiday spirit to the digital world, we receive a mere 15 likes on our Instagram posts (20 if you count the x-posts on Facebook)

Maybe it is unrealistic to expect a social media uproar if you have a mere 300 Instagram followers, but surely these heartfelt Christmas posts must compel more than 5% of our friends to react, right? The problem is, your total follower count is not the same as your “true audience size” — the number of connections that are active and see your posts. You may have many followers, but the majority may not be active, or they may have too dense a follow graph that your posts realistically never reach them.

Then the question is, how can you measure your true social media reach, your Internet clout?

The Lincoln-Peterson index — or how to count uncountable things

Let’s say you make an Instagram post on a private account of 300 followers and receive 40 likes. You know your audience is somewhere between 40 and 300. If you have supreme confidence in your ability to craft Instagram posts, you may suppose your reaction rate be 100%, your audience a mere 40 users and it was Instagram who squeeze your reach.

But maybe your photos don’t look too good, or maybe you are not as witty as you think. Maybe you did reach your full 300 followers but you can only get 5% of them to gingerly give you a virtual moral support. There’s no way to know with one post. But if you have two posts, you can get a good idea, even if you don’t know how compelling each individual posts are.

Let’s say your second post has 60 likes and there are 15 mutual likes (people who like both posts). We can naively calculate our true audience size using the good old Set Intersection theory, knowing that it is very conservative estimation, because we potentially neglect a population that doesn’t like either of the Instagram posts.

Fortunately, there is a neat method called the Lincoln Index that allows us to estimate the total audience. Let n1, n2 be the number of likes on the first and second post respectively, and m the number of people who like both:

Our estimated reach according to the Lincoln index is 40 x 60/15 = 160. Notice that it is a lot larger than the naive estimate of 85, which exonerates Instagram from squeezing our social clout.

How does it work mathematically? — or the part where you should feel free to skip

Let us re-frame the problem: we want to choose n2 people to like the second post out of N people, where exactly m out of n2 are among n1 people who like the first post. It should be fairly straightforward that the probability to arrive at such choice is:

We can use the maximum likelihood method to find an estimate N. For the mathematically inclined, here is an outline of the proof:

• The trick is instead of trying to find N that maximizes the likelihood function:

• we find the largest 𝑁 that satisfies the inequality:

• it should be obvious, after a bit of arithmetic juggling, that this is equivalent to:

• And so the estimate N_hat that maximizes the likelihood function has values:

A simulation

I know not everyone is a fan of mathematical proofs so I will try to convince you with a small simulation.

Let’s say we know our social media reach to be 100 and we have 2 Instagram posts with the probability of like p1 and p2, respectively. We can simulate different scenarios of p1 and p2, observe the number of likes on each posts, then reapply the naive Method and the Lincoln methods to see how far each of them are from arriving at the value of 100 (ground truth).

We can see that at different values of p1, p2, the Lincoln estimates remain unbiased (the mean of the 10,000 simulations approximate the true population of 100), while the 95% confidence interval reduces as p1, p2 approximate 1. This makes sense i.e. if we are confident that the like rate of our posts is close to 100%, we should be fairly confident about our audience size estimate. The naive estimate, unfortunately, underestimates our social reach by a lot.

As usual, all good statistics come with some caveats:

• (1) We assume that between each post, the true audience size does not change. This may be true for a casual Instagram user, but probably unrealistic if you are an influencer with an ever increasing fan base.
• (2) We also assume there are overlapping mutual likes between the two posts. Some of our readers may have found a small cheat in my code:

While the mathematics remain valid, the estimate becomes fairly rough and biased at smaller sample size. Luckily, we have a neat modification, called the Chapman estimate, to address this issue:

So what is this counting method useful for? This statistical technique— estimating any population size where is impractical to count every individual is very popular in ecology, under a commonly known name mark-recapture method. It is also common to see this method in action in Tech companies, where we want to estimate the amount of bugs in a video game before release, or the amount of policy-violated contents in a social media platform.

I was once challenged in an interview to estimate the number of Man City fans in London. It would have been a fun weekend activity to do some sampling in Trafalgar Square, interviewing people and crunching numbers according to the Lincoln estimate. Although I suspect that the City fans are so few you may actually be able to just count them all!

If you enjoy this article, you may also enjoy my other article about interesting statistical facts and rules of thumbs

• Disney Movies were right — We are all special, and Statistically so
• A statistical rule to optimize your life: the Lindy’s Effect
• Rules of Three: Calculating the probability of events that have not yet occurred

For other deep dive analyses:

• How Bayesian Statistics convince me to hit the gym?
• Making big bucks with a data-driven sport betting strategy
• Love in the time of Pandemic: A Probabilistic approach to Dating

The entire code for this project (simulation and graph) can be found in my Github.

The statistical theory behind why your Instagram posts have so few likes was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.

Dating is hard — but maybe Probability Theory can make it slightly easier this cuffing season

Back in 2019, I wrote an article about Optimal Stopping theory and its application to dating.

Finding the optimal dating strategy for 2019 with probability theory

The lesson was simple for the hopeful romantic: to maximize your chance of finding the best match, you should get to know but not commit yourself to any of the first 37% of the people (you think) you will meet, and then immediately settle for the first person who’s better than the ones you saw before. Doing this will guarantee a 37% chance for you to find the Mathematically Destined One.

It has been 3 years since then and it was enough time for me to learn that reality is nothing like these mathematical formulae.

However, reality is harsh.

I graduated from college and moved to a new city for work. And let me tell you, even at the age when everything seems complain-worthy, from the long waits at Blue Bottle to the non-existence of public transportation in the SF Bay Area and the exorbitant price I pay for my tiny 1-bedroom, the lack of time and opportunity to date is still at the top of my gripes.

And then the pandemic hit and abruptly changed the definition of romance, from intimate experiences with dinner, a few drinks in a crowded bar or a goodnight kisses to awkward Zoom call, mask-wearing and the threat of catching a deadly disease. It’s virtually impossible to date a person, let alone 10 or 37.

And don’t even get me started on how hyper-competitive the dating market is in San Francisco or New York City. Even if (you think) you have found that significant other, chances are that the person is not available or does not reciprocate your feelings.

So the Math seems way out of touch with reality, what now?

We follow what good scientists would do: we acknowledge the shortcomings of our model, relax a few assumptions and come up with the new algorithm to better reflect this reality.

1. Finding the Many Ones — a top k approach

One way to make dating a tiny little bit easier for us is to ask ourselves: is the best match really that different from the “second best” or “third best”? After all, we are probably so infatuated by love that committing to a serious relationship with any of the top 3 or top 5 candidates will likely make very little difference to our lifetime utility.

So, let’s say, instead of finding the best match among 100 people, we set out to find the top k love candidates. We use the same approach: at first, we explore our compatibility with the first r candidates and make a mental note on the maximum score observed among them. Let’s call this threshold t. We then select each candidate from the interval {r + 1, …, n} with score higher than t until we select k persons or we run out of candidates.

A simulation:

Let us also start with running a simulation to build up some intuition around the optimal value of r. The setup is fairly simple with the following Python code:

Let us visualize the result when n = 100:

We can see that k = 1 turns this into the classic problem where we only pick the best person. The ideal exploration phase will thus last for around 36–37 people. Interestingly, we have to conclude our exploration phase a lot sooner to only around 23 or 24 people when we want to choose the top 3 candidates, and to 9 people when we want to cast our net wide for the top 10 candidates.

You may ask, the probability of success actually decreases as we increase the number of top partners we want to pick (around 15% when k= 3 and 4% when k = 10), so is it actually worse than our original strategy, where we are guaranteed a success rate of 37%?

• That is true, but the catch is, the original magic number 37 only guarantees that you will find the best one, but it says nothing about their availability. It’s true that as k increases, the probability of successfully choosing top k partners decreases, but the chance that at least one of them accept is significantly higher. After all, in love, it is only one that counts.

How does the Math actually work?

We can also run some math to find r.

Let S(n, k, r) be the event that we succeed, that is, we find the top k candidates with n by using our algorithm on the first r observations. We want to maximize the probability of success i.e. P(S(n, k, r)).

To do this, we have to use a small trick: we break down P(S(n, k, r)) into the sum of small components P(S(n, k, r, i)), where S(n, k, r, i) is the event that we succeed with i being the last selected candidate.

To break P(S(n, k, r, i)) even further, notice that for i to be the last selected candidates:

• It has to be one of the top k candidates, with probability k/n
• We were able to pick k-1 candidates among i-r-1 candidates:

• One of the first r people is the best among the first i-1 people, with probability r/(i-1)

Our equation becomes:

Analytically finding r to maximize this equation is extremely complicated and beyond the scope of this post. But, if you trust me, then as n → ∞, the optimal r approaches:

And if you need some convincing, this is what the r_best looks like between our simulated and theoretical results:

2. Handling rejection — when you are not The One for the other

Let’s say that we are still committed to finding the best person, and that person only, but we realize that there is a (pretty high) chance that the person is not available, e.g., they may not be attracted to you, they are already in a relationship, or … you happen to be in their exploration phase.

So let us tweak the problem a little bit and allows the candidate the right to refuse our offer of love. We assume that if given our proposal, any candidate will reciprocate our feelings with the same probability p, independent of their rank and the disposition of other people we interact with. This is obviously naive, but the assumption will do for now.

The math is extremely messy so I will save you from it. Basically, in this case, the strategy is also to explore r candidates first, where r = p^(1/(1-p)). For example, if we think the acceptance rate to be 0.5, then we should explore r = 0.5^(1/0.5) = 25% of our candidates and stop with the first available applicant who is better than those preceding them.

You can see that in the simulation below:

Depending on your success projection, you may be able to explore more or fewer candidates. At 25% successful rate, you can only explore around 15 candidates. If you live in the Bay Area where I think (I have no data to back my claim) the availability rate is much less than 10%, then you can only afford exploring 7 candidates. That is, quite a grim outlook.

Beyond Love and Dating:

At this point, you probably can tell that even with these revised assumptions, it’s impossible for us to capture the complexity of human interaction that goes into love and dating.

However, this does not mean that what everything you have read so far is just fun and games. What you have been introduced so far is the very backbone of Optimal Stopping Theory, where we need to choose a time to take a given action based on sequentially observed random variables in order to maximize an expected payoff. In fact, you can basically swap “finding a partner” with “hiring a data scientist,” “finding a person to sell my house to” or “option trading” and the math would be the same. In particular, the top k problem you have seen is a well-known ML-system problem called Online Data sampling, where you need to perform online data selection and sampling from a stream of data without revision.

But then again, none of them is as fun as dating, isn’t it?

If you enjoy this article, you may also enjoy my other article in a series about interesting statistical facts and rules of thumbs

• Disney Movies Were Right — We Are All Special, and Statistically So
• A statistical rule to optimize your life: the Lindy’s Effect

For other deep dive analyses:

• How Bayesian Statistics convince me to hit the gym?
• Making big bucks with a data-driven sport betting strategy

You can find all the codes to my article at my Github link.

Love in the time of Pandemic: A Probabilistic approach to dating was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.

A little bit of Statistical Theory will convince you that the idea of “the average Joe” is really just a big myth

I love Disney movies !!!

Besides all magical kingdoms and thrilling adventures, mystery and romance, the one thing Disney Movies are so good at is convincing me that, after all, I am special.

Nowadays, we live by comparisons to others. We want to live an exciting life, to be more than just an “average Joe” or “Jane,” and if you work in tech, to be that 10x engineer friend who make millions from a start-up exit and afford a vacation home in Tahoe. So among all the movies constantly reminding us of how insignificant we are in the monotonous capitalistic rat race, there is a deep sense of comfort on a Sunday evening when we rewatch movies like Toy Stories or Camp Rock. Isn’t it great to see Woody or Mitchell (Demi Lovato), through a series of uplifting songs, recognizing that they are special who they are, and they deserve as much love and attention just for being themselves? Then Monday hits and you tell yourself as you are brought back to the reality of the rat race: “That was fun, but it was movie after all.”

But, what if I tell you that, statistically, these movies are right, that it’s almost impossible for us to be “normal,” and “average”, and that each of us is really special who we are?

The misleading nature of the Normal Distribution

Most of us are familiar with the Normal distribution and the famous 68–95–99.7 rule, that is, 68%, 95%, and 99.7% of the population will lie within 1, 2 and 3 standard deviation from the mean. In other words, almost all observed data is extremely close to the mean value.

A popular example is height. The heights of adult men in America (>20 years old) are empirically normally distributed, with a mean of 5'7’” and a standard deviation of 3 inches. So it shouldn’t surprised you to see most people around you is between 5'2" and 6'2'’(2 std from the mean).

This would make you think: by the same logic, it should also applies to other attributes like weight, talents, education level, etc. right? It seems quite likely for any of us, including myself, to be average in everything, with really nothing special to show.

Well, there’s one thing you forgot to factor in: dimensionality.

The “blessing” of dimensionality

Let’s revisit the problem but in mathematical terms: suppose each person’s has n factor, each of which is represented by an independent, standard normal random variable. That is, each person’ attribute vector has the form (x1, x2, x2, …, xn) where each xi i.i.d Normal with mean 0 and variance 1. The expected value of each factor is 0, so you might expect that most points are piled up near the origin (0, 0, 0, …, 0), right?

One way to see how far we are from the average (0, 0, 0, …, 0) is to look at the Euclidean norm. I used the squared Euclidean norm for computational simplicity but the idea is the same. Particularly, we want to find the distribution of this distance, P(d²) and see where the distribution mass is maximized or where most of the points are located.

Since d² comes from a sum of square of n-i.i.d standard normal distribution, d² follows a Chi-squared distribution with n degrees of freedom. To study the behavior of the n-degree of freedom Chi-squared distribution, we can approximate it by a Normal distribution with mean n and variance 2n.

Since it is hard to think about how d² behaves especially when n increases, we can also change our variable and look at d²/n. Intuitively, we can think of multiple spheres with different radius value surrounding the origin. We particularly care about the spheres with radius n, and we want to see how the points are distributed, relative to this n-radius fictional spheres. It is easy to see find the distribution for d²/n as:

Since d²/n follows a Normal distribution, we can also apply the 68–95–99.7 rule:

This means that 95% of the points lies just between [1–2 * σ, 1+2*σ] from the imaginary sphere of radius n we mentioned earlier. As n increases, this spaces gets smaller and smaller i.e. the points lie much closer to the sphere. When n → ∞, all the points lie on the imaginary sphere.

This is surprising, because initially, we thought most of the points will be close to the origin. In fact, most of them are n-distanced away from the origin, in a n-dimensional soap bubble with radius n.

Initially, we thought most of the points will be close to the origin. In fact, most of them are n-distanced away from the origin, in a n-dimensional soap bubble with radius n.

A simulation

I know not everyone is a fan of mathematical proofs so I will try to convince you with a small simulation.

For different values of n, we sample a large amount of “people” — where each person comes from a n-dimensional Normal distribution. We then calculate the distance (squared Euclidian norm) from each person to the “average Joe” or the origin. The graph here shows the distribution of distance when n = 100.

The two red lines shows 2 standard deviation away from the mean and where where 95% of the points exist. In the graph below, we see that as n → ∞, the ratio between the boundaries and n approaches 1, suggesting most points lie in a very small area n-distance away from the origin.

Another way to look at this is shown on the right graph: by the time we reach 2000 dimensions, nearly all of the vectors sampled from the distribution have a magnitude near 1, and virtually no points are close to the origin.

Some caveats:

So I hope I have convinced you that: somebody exists who is more or less the average height, or weight. But if there are 100 or more roughly independently varying criteria to measure people by, then probably nobody exists who is average in every way; there are too many ways to be different.

However, this analysis relies on the assumption of independence. We need to come up with a lot of independent factors to convince ourselves that we are special. That is not necessarily the case for most of the common attributes, e.g. height and weight can be correlated through genetics, your ability to do various sports correlate through your physique, etc.

But, at least for now, whenever you saw another cheesy Disney movie that tell you to believe in special self, at least you know there is some statistical backing to that claim.

If you enjoy this article, you may also enjoy my other article in a series about interesting statistical facts and rules of thumbs

• A statistical rule to optimize your life: the Lindy’s Effect
• Rules of Three: Calculating the probability of events that have not yet occurred

For other deep dive analyses:

• How Bayesian Statistics convince me to hit the gym?
• Making big bucks with a data-driven sport betting strategy
• Love in the time of Pandemic: A Probabilistic approach to Dating

The entire code for this project (simulation and graph) can be found in my Github.

Disney Movies were right — we are all special, and statistically so was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.

Sometimes, I wake up and feel anxious to open my phone. Scrolling through all the social media news feed, glancing at the New York Times app, hoping around some Youtube new uploads, takes a lot of time.

The flood of information that swamps me daily produces more pain than gain. Every iota of information has been so well-optimized at ‘looking urgent/important’ that we have to constantly making choices in terms of what to process, knowing that little (if none) of them will remain relevant in the next 24 hours?

What if there existed a rule to identify what type of news or information that will still be relevant for a long time to come?

What if there was a guide you could use to filter the knowledge worth knowing from the rest of the noise?

There is — it’s called the Lindy Effect.

Lindy’s effect — or the value of Time as a judge of quality

The Lindy’s effect is an observation that:

The longer an artifact has been around, the longer it’s likely to stay around.

Nassim Nicholas Taleb, in his book, Antifragile, explains the concept quite cleanly: “If a book has been in print for forty years, I can expect it to be in print for another forty years. But,[…], if it survives another decade, then it will be expected to be in print another fifty years.”

The idea is, the future life expectancy of some non-perishable things like a book/technology/intellectual production is proportional to their current age; every additional period of survival implies a longer remaining life expectancy.

A mathematical motivation:

The Lindy’s effect starts with an idea that the lifetimes of some intellectual artifacts follow a Pareto distribution: a skewed distribution with heavy, “slowly decaying” tails (i.e. much of the data is in the tails).

You may have heard of the Pareto distribution as the famous 80/20 rule: there are a few items that are very common and contribute a lot to the “head” of the distribution e.g. 20% of the top-earners account for 80% of the national wealth; the top football teams represent a small fraction in the English leagues, yet they account for the bulk of TV royalties.

Modeling the lifetime of artifacts by the Pareto distribution means most items will last only for a short period of time, but there is a long tail of items that could last substantially longer. Lots of economic research has verified this property among intellectual artifacts.

We can write the probability distribution function of lifetime — PDF(t) of such an artifact under the Pareto form:

where t_m (the scale parameter) is the minimum possible value of t, and α is the shape parameter, also known as the tail index. Notice that t_m is necessarily positive i.e. the lifetime of an artifact has to be > 0.

Now, we want to find the expected lifetime E(t), conditional on it having lasted a period of time t_0. In other words, we want to find E(t|t>t_0).

and so the expected additional life after t_0, or E (t-t_0|t>t_0) is:

While the math can look intimidating, the result is surprisingly clean. The conditional expected additional lifetime is simply the time that has passed (t_0) multiplied by a proportionality constant 1/(α-1).

• The closer α is to 1, the larger the proportionality constant, and the larger expected additional lifetime we can hope for our artifact. If α = 2, the proportionality constant is 1. That is, the expected additional life equals the life seen so far.
• For example, when α = 2, once t_0 is reached, we expected the future lifetime is also t_0; if that point is reached making the total lifetime so far 2t_0, then the new expected future lifetime is 2t_0; and so on.

A simulation

I know not everyone is a fan of mathematical proofs so I will try to convince you with a small simulation.

Let’s assume that 2021 is a 20th anniversary of your favorite book and you want see how many additional years the book have before completely falling into oblivion.

What we will do is to simulate many, many scenario of lifetime from P(t-t_0|t>20, α) for different values of α, and see how the distributions behave.

Now, actually performing this simulation requires a technique called Inverse Random Sampling. Its details are beyond the scope of this post, but if you are interested, I wrote a small note explaining the technique .

You can see that as α increases in value, the conditional distribution becomes “heavier” around the head, thus resulting in a smaller conditional expectation of additional years remaining. A book that has been in print for 20 years is expected to last another 40 years if α=1.5, 20 if if α=2, and 10 if α=3.

One interesting result is to compare how the conditional expectation of remaining years differ between a Lindy’s good and a human life.

With some quick research, I found that US males lifetime distribution follows a Weibull distribution with a mean expectancy of 70. Given this information, we can then do some simulation to see what the conditional probability of remaining lifetime after one has celebrated their 20th birthday:

What we learn is technologies have different survival patterns than people.

• Given the first 20 years, the expected remaining year for a Lindy’s good is 20, while we do expect an average US male to live for another 50 years (which adds up to the expectancy of 70).
• However, the older a person is, the fewer expected years he has left. That is because human lifetimes follow thin-tailed distributions — much thinner than the Pareto distribution’s tail. The Pareto distribution’s heavy tail is the reason why the longer a technology has survived, the longer we can expect it continue surviving. In Taleb’s words, technology “ages in reverse.”

Applying Lindy’s effect to optimize your life:

In essence, Lindy’s effect teaches us a simple principle: if you cannot decide, then let time be the natural filter for contents:

• I have uninstalled all the news app from my phone, while subscribing to weekly recap newsletters. If an event is important and newsworthy, then it should remain so beyond the first 24 hours.
• Books: not knowing which books to choose between all presented in Bill Gates’ reading notes or NPR’s Book Concierge? Read the Classics. The Bible, Machiavelli’s The Prince, Marcus Aurelius’s Meditations, the insights of these books have tested against changing time and they most likely outlive most of this year’s NYT’s bestsellers.
• ML algorithms: Linear regression via least squared estimation, invented by Gauss in the 1800s, is still by far the most popular Machine Learning techniques among data scientists, according to a 2016 Kdnuggets survey. If I have to place a bet between least squared linear regression and deep neural nets and choose which algorithm to still be used in the year 2200+, my bet is definitely with the former.

Some caveats:

As we have seen, Lindy’s effect is very neat probability statement. But they are exactly that, probability statement, and you shouldn’t take it too seriously:

• We don’t know for sure which item exhibits Lindy’s effect. The condition for it is rather strict. Not everything follows Pareto distribution: human lifetime certainly does not, neither does atomic decay. On top of that, even if an artifact’s lifetime is distributed in a Pareto manner, we have seen that there’s a big difference in additional lifetime expectation if we just move the shape parameter a little bit.
• Note that I computed E(t|t>t_0), i.e. it only conditions on knowing that t>t_0. If you have additional information, such as evidence that a technology is in decline, then you need to condition on that information.

The Lindy Effect is thus a trivial heuristics, the usefulness of which decreases as you acquire additional relevant information.

But you know, to filter books and news to read, maybe trivial heuristics is good enough.

If you enjoy this article, you may also enjoy my other article in a new series about interesting statistical facts and rules of thumbs

• Disney Movies Were Right — We Are All Special, and Statistically So
• Rules of Three: Calculating the probability of events that have not yet occurred

For other deep dive analyses:

• How Bayesian Statistics convince me to hit the gym?
• Making big bucks with a data-driven sport betting strategy
• A Bayesian quest to find God

The entire code for this project can be found in my Github.

A statistical rule to optimize your life: the Lindy’s Effect was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.

How confident are we that blue states California, New York, etc. do not turn red in the next election?

Author’s note: I’m experimenting with a shorter form of articles, which aims to introduce interesting statistical facts and rules of thumbs. Let me know how you like them. For more in-depth analysis:

• How Bayesian Statistics convince me to hit the gym?
• Making big bucks with a data-driven sport betting strategy
• A Bayesian quest to find God

Another election cycle, another time we are bewildered by the overall results. There was no blowout for any side. The election was a nail biter. Both Democrats and Republicans made unexpected inroads in several races.

Facing the rising uncertainty each election cycles, a friend of mine acclaimed:

The one thing we know for sure is that California has never, and will never turn Red.

Well, for one, it is historically incorrect — California had been consistently voting Red before the ’92 election¹. My friend would have to qualify his observation to “elections that he observed [since 1992]”, or if he changed the place to District of Columbia. However, he will still be less confident about his projection if he knew that even in this election, a third of the 10 million votes in California was casted for President Trump². Just because we [seemingly] observe a place to be a Democratic stronghold does not mean that it cannot turn Red in the future and vice versa.

The type of projections my friend made is a very common mental induction that we all exercise on a daily basis. Just as airplanes should not crash, wisdom teeth removal should not cause nerve damage and a reigning Premier League champion should not be upset by a team wallowing in the fifth tier of English football. In reality, however, these things do happen. Just because I have not ran into any accidents on a motorbike for the last 10 years does not mean the risk of accident is 0.

The question is: can we quantify the probability of something that has yet to occur? Ideally, we would want to assign a probability to the rare event, but the numerators for these events are 0, and our hands are tied. A convenient way to quantify the amount of risk in these cases is called the Rule of Three:

If a certain event did not occur in a sample with n subjects, the interval from 0 to 3/n is a 95% confidence interval for the rate of occurrences in the population.

One thing to note here is that the smaller the sample observed, the wider the confidence interval, meaning the upper limit of a confidence interval from a small sample is greater than that from a large sample. But it would be wrong to say that the true probability of an adverse event occurring is larger in a small series.

How does it work?

The proof of this rule is quite simple, denoting the number of events by X and the probability that we observe an adverse event by p (p is close to 0), we want to find the values of the parameter p of a binomial distribution of n observation that give Pr(X = 0) ≤ 0.05.

And so we need to solve for p such that:

Then we can approximate the entity ln(1-p) by using the Maclaurin series, which is the evaluation of the Taylor expansion of function f(x) at a = 0

Applying this to ln(1-p), we get:

Thus, we have -np ≤ -2.996 or p ≤ 3/n.

A quick simulation:

While the math is clear and simple, to demonstrate that it actually works, we can try to run a simple simulation:

• Imagine us testing for defects in a production line. Let’s say we know the probability of defect is p = 0.001
• We test every products that come off the line, one by one, until we come across a defect (test number n+1).
• We then make a confidence interval bounded by 0 and 3/n and make note of it. We make check whether the confidence interval contains the true probability of defect p = 0.001.
• We repeat this process for a lot of time e.g 100,000 iterations

This is a snippet of the R code:

As you can see, in the long run, about 95% of the interval [0, 3/n] such intervals should contain the true underlying probability of defect.

Back to the original question:

So let us come back to the original questions: let’s forget all the elections before 1992. Since we have not observed California voting Red for the last 8 elections, using the rule of three, we can deduce that the 95% confidence interval of California going red is between 0 and 3/8 = 37.5%. For District of Columbia, it would be [0, 3/15= 20%].

To be clear, this is a fun example, and to make it work, I have to make some sweeping assumptions here, none of which is probably true.

• (1) We assume that each election result is independent of one another. This is largely untrue, as we know an election result can at least influences district mapping or motivates the turnout of the next one, etc.
• (2) There is an underlying, fixed probability of a state voting Red, this is probably untrue, as we know that the state demography, economic outlook and people’s sentiments, which largely influence party affiliation and how the result come about, change all the time.
• (3) For California, we ignore all elections before 1992. If we included data before, where the rare event did occur [but it was not that rare since California was considered a Red state before 1990s], the rule of three should not have been used [D.C does not have this problem however]

An actual use case:

So what is the Rule of Three actually used for? The very popular use case is in pre-licensure clinical trials where the adverse event rate is very rare e.g. vaccine development program. Usually, a sample size in the range of 10,000 subjects can be considered adequate for establishing the protective efficacy of a new vaccine, but not sufficient for detecting the rarer events or complications.

If no event of any particular sort has been observed during a pre-licensure clinical program involving 10,000 individuals, it can be estimated that this event rate has an upper limit of 3 per 10,000. The actual complication rates, those occurring at a lower frequency than the vaccine-targeted disease, can only be detected by some reliable post-marketing surveillance systems.

So the next time your dentist says he never fails a procedure, you know what to do: ask him how many patients he has served, used the rule of three to find the confidence interval for his failure rate, and see how comfortable you are with the upper bound.

You can find the code for the article at my Github page.

References:

[1] Why Is California a Blue State?

[2] How this year’s ballot measures show California’s growing political divide

The Rule of Three: calculating the probability of events that have not yet occurred was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.

Arguably the most nervy moment of all during a football match, the penalty shot is psychological warfare: Two men, one ball and a clear chance to score a goal. A penalty is supposed to be finished at the professional level: period. Anything else is a failure.

Such a task should only be required of players made with the iciest of veins.

There are great penalty-takers, those we can always count on to deliver the ball to the back of the net with a staggering metronomic efficiency. And there are those whose dismal effort, trademarked tip-toe run up or near-orbital strikes are only source for our deep regrets or relentless ridicule.

You would think taking penalties is easy. Well, maybe spend 5 minutes on this video.

So who is the true Man of Steel?

You would imagine that this would be a simple secondary Math problems. We collect the records of all spot kicks, calculate the corresponding ratios of conversion for all these players, and TADA!

But let us look at real examples and see how this simple logic falls through !!!

On one hand, we have Cristiano Ronaldo, five Ballons d’Or under his belt, whose nickname ‘Penaldo’ stands for his unwavering self-composure in front of the goal with an impressive records of 85 goals out of 98 penalty kicks (86.7%). On the other hand, we have a young talent, Nicolas Pepe at Arsenal, whose recent goal against Aston Villa made his tally totaling 10 successes out of 11 (90.9%). Is Pepe better/more reliable at penalties than Ronaldo?

Of course, saying anything like that would be ridiculous !!!

While the young player has slightly higher proportions of scoring, it’s not a lot of evidence: a typical player would score almost 80% of the time, and Pepe’s 90.9% success rate could be due to luck. Ronaldo, on the other hand, has a lot of evidence to show that he is above an average player.

How do we systematically compare and make meaningful judgement about ratios like 85/98 vs 10/11?

Bayesian Statistics to the rescue:

And it’s not in just football that we find this dilemma. We are bound to get tripped up working with successes/totals ratios due to the uncertainty in low counts. The simple solution is to set some thresholds and filter out those players who do not meet the minimum requirement. But it is far from ideal: we are throwing away what could be useful information.

Let me take a step back and re-frame the questions: we want to somewhat gauge the innate penalty-taking ability of a player. His historical record is a mere estimate of this ability. The task is to more accurately estimate his ability given his historical record.

I hope you can see how this naturally leads us to a classic Bayesian setting: we can start off with some prior expectation of the player’s ability, and by updating it with new evidence — that is, with his actual scoring record, we get arrive at a more practical estimation of his true ability.

For a good review of the Bayesian framework, you can have a look at my other articles:

• How Bayesian statistics convinced me to hit the gym
• A Bayesian quest to find God

Prior expectation and Empirical Bayes:

I’m stepping on a landmine right here. Although the Bayesian setting makes intuitive sense, it is almost impossible to get people agree on a reasonable prior to use for each setting.

How might one go about specifying their prior belief? Especially if you don’t have an informed guess? One way might be to assign a “flat prior” — a prior where you specify that every value of your parameter is equally likely. Another way is to use Empirical Bayes, where your data is used to estimate the prior. Here’s the kicker: instead of having a fixed notion of the prior ahead of time, I will save myself from all the trouble by taking a look at the data after collection and estimate the prior distributions accordingly.

There are advantages from this approach. It is simple, and it takes away (well, sort of) the subjectivity that can creeps in should I have to choose my own priors. Personally, I found that this approach helps me achieve a degree of stabilization of results. A flat prior or something out of line with the data could lead to strange results, e.g negative posterior probabilities. And believe me, this is not uncommon. I’m more likely to get sane answers with empirical Bayes.

Listen to what famous statistician Bradley Efron has to say about this method:

Ever since my first exposure in the 1970’s, via Robbins and Stein, empirical Bayes has been the most riveting topic for me. It still seems like magic sometimes. [1]

David Robinson — chief Data Scientist at DataCamp, also offers some good insights here:

It basically comes down to how many observations we have: if we have a lot, we can get a good estimate that doesn’t depend much on any one individual. Empirical Bayes is an approximation to more exact Bayesian methods- and with the amount of data we have, it’s a very good approximation.[2]

I know, it feels like cheating, and rightly so. Using actual data to estimate the prior on top of subsequent use for the inference is criticized by many who consider the prior information exogenous to observations, as it can underestimate modeling errors. Any data is going to be complacent with a model which used the same data to specify some of its features.

Finding the priors — The Beta Distribution:

Let us take a look at the data.

I have scraped 13 year-worth of penalty kicks data from 6 different leagues across Europe. This is how it looks:

As explained earlier, to come up with a good approximation for the priors for the success ratio, we will take a look at its actual distributions.

Here, our task is to find a distribution that characterizes our empirical penalties conversion rate: that is, (1) it takes in a domain from 0 to 1 (no negative probability), (2) it needs to be flexible, in the sense that it could represent what we see in our histogram, and (3) if it makes our Bayesian calculation easier when we combine with the likelihood from the binomial distribution — which describes the number of successes within n trials (i.e: number of goals out of many penalties kicks), all the better.

Lo and behold, the Beta distribution.

The Beta distribution takes on the domain [0,1], so we know it cannot do worse than the flat prior. It is highly flexible within the domain, being characterized by two positive shape parameters, denoted by α and β (see figure above).

The best part about the Beta distribution is that it is a conjugate prior to the binomial distribution, and this works wonder to our calculation. We can start off with Beta distribution B(α , β) as our prior, updates with information from the Binomial distribution with s successes (scores) and f failures (misses), and arrive at the posterior distribution which is also a Beta distribution of a precise form B(α + s, β + f). The mean of such distribution is (α + s)/(α + β +s+ f)

Simply Magical !!!

I include some math at the end of this post just in case you are not convinced !

Let’s go on and find our best penalty taker:

Estimate the priors for our actual data:

Now that we have a clear framework for our analysis, it’s time to apply it to our real dataset.

First, let’s fit a Beta distribution to our conversion rate dataset. This could be done quite simply with R.

Using this method, I get an estimation of α = 9.29 and β = 2.93. Let see how well this estimation fits our data.

Obviously, it is not perfect, but something we can work with for now.

Updating each player statistics using our empirical priors:

We can now go on and find the posterior distributions for each player’s conversion rate using the formula I mentioned. That is, a Beta distribution of the form B(α + s, β + f), where s and f are the numbers of successes and failures, respectively. The posterior mean will be (α + s)/ (α + β +s+ f).

For example, the posterior distribution mean for Ronaldo and Pepe’s penalty conversion rate will be:

We can plot their posteriors distributions (with some other players’ ) as follows:

We can see that Ronaldo’s strong records results in very high posterior mean and a narrow credible interval, while Pepe’s posterior mean conversion has a much larger credible intervals, as he has less data to stand for. Meanwhile, Griezmann and Sanchez are just not that good when it comes to penalty.

The top 10 penalty takers:

After applying the empirical Bayes estimate and beta-binomial posterior calculation, let me introduce the top 10 penalty takers over the last 13 years.

So it turns out that the best penalty takers, at least by this method, are Salihovic, and G.Rossi, followed very closely by Boudebouz, Mark Noble and Lewandowski. The famous penalty takers: Lampard, Gerrard, Harry Kane are not as good as they are claimed to be. Interestingly, we can see that for these players, their conversion rates shrink after we apply this procedure. We can plot the original and the adjusted ratios on a plot to assess the adjustment. This section is inspired by an amazing great article by David Robinson [2].

The horizontal dashed line shows the mean of the prior distribution — what we should guess a player’s success rate is if we don’t have any information about the player. Then, we have the diagonal line y = x, showing the data points in which the actual goal scoring average equals the posterior adjustment.

There are two interesting observations. First, we can see that the points above the prior mean tend to shrink down, while the points below get inflated up towards the mean. Second, the points that was least affected by the shrinkage are points associated with a lot of observations. And when you think about, it makes perfect sense: if we have a lot of evidence through observation, we can rely on the actual estimation; otherwise, there is a lot of value in believing that our data tends towards the population mean.

Conclusion:

So there you go, I have just introduced the concept of empirical Bayesian estimation: we can leverage the overall distribution of our data as a prior to arrive at a more realistic estimate of each average. Then, we use the powerful concept of conjugate priors to bypass the complications of calculating the posterior distributions in our procedure. The beta-binomial combination are so popular that you can apply it to a host of scenarios: basketball, A/B testing with clickthrough, school canteen rating, etc.

For the purpose of this article, I have conveniently ignored many complications. For one, we can tell that a beta distribution for our prior is barely a great fit. In fact, we can fit a bimodal distribution with 2 Beta distributions, assigning every player in two different groups according to the likelihood that they belong to each distribution.

You can checkout my code that address this issue.

Another complication is that the penalty conversion rate changes over time and differ across different leagues. As far as we consider all these different factors, we take many steps away from our simple Beta-Binomial distribution and enter the realm of Bayesian hierarchical model. And that’s where the fun begins !!!

Appendix:

The code for this article and more can be found here.

This is a short proof for the Beta-Binomial conjugate pairing.

The setting is simple: we observe a number of success counts s that we believe follows a Binomial distribution with sample size n and success rate r. Our job is to find the posterior distribution of r given s.

Time for us to call upon the famous Bayesian formula — only for continuous variable:

The likelihood function p(s|r) and the prior p(r) follow well-defined Binomial and Beta distributions.

Then, our job is just plugging in and finding the posterior distribution.

The result is a Beta distribution of the form B(α + s, β + f)

[1] Statistic View. (2015, June 1). “Empirical Bayes has been the most riveting topic for me. It still seems like magic sometimes: An interview with Bradley Efron.”

[2] Robinson, D. (2015, October 1). “Understanding empirical Bayes estimation (using baseball statistics).”

Men of Steel: Finding the best penalty-takers with empirical Bayes estimation was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.