It’s been a while since anyone said “killer app” without irony, so let me remind you that a killer app is software “so necessary or desirable that it proves the core value of some larger technology,” quoth Wikipedia. For example, most people didn’t have much use for the internet until the world wide web was populated with useful content and the first generation of browsers made it easy to access.
So what is the Bayesian killer app? That is, for people who don’t know much about Bayesian methods, what’s the application that demonstrates their core value? I have a nomination: Thompson sampling, also known as the Bayesian bandit strategy, which is the foundation of Bayesian A/B testing.
I’ve been writing and teaching about Bayesian methods for a while, and Thompson sampling is the destination that provides the shortest path from Bayes’s Theorem to a practical, useful method that is meaningfully better than the more familiar alternative, hypothesis testing in general and Student’s t test in particular.
So what does that path look like? Well, funny you should ask, because I presented my answer last November as a tutorial at PyData Global 2022, and the video has just been posted:
This tutorial is a hands-on introduction to Bayesian Decision Analysis (BDA), which is a framework for using probability to guide decision-making under uncertainty. I start with Bayes’s Theorem, which is the foundation of Bayesian statistics, and work toward the Bayesian bandit strategy, which is used for A/B testing, medical tests, and related applications. For each step, I provide a Jupyter notebook where you can run Python code and work on exercises. In addition to the bandit strategy, I summarize two other applications of BDA, optimal bidding and deriving a decision rule. Finally, I suggest resources you can use to learn more.
Outline
Problem statement: A/B testing, medical tests, and the Bayesian bandit problem
Prerequisites and goals
Bayes’s theorem and the five urn problem
Using Pandas to represent a PMF
Estimating proportions
From belief to strategy
Implementing and testing Thompson sampling
More generally: two other examples of BDA
Resources and next steps
Prerequisites
For this tutorial, you should be familiar with Python at an intermediate level. We’ll use NumPy, SciPy, and Pandas, but I’ll explain what you need to know as we go. You should be familiar with basic probability, but you don’t need to know anything about Bayesian statistics. I provide Jupyter notebooks that run on Colab, so you don’t have to install anything or prepare ahead of time. But you should be familiar with Jupyter notebooks.
My audience skews left; that is, the people who read my blog are more liberal, on average, than the general population. For example, if I surveyed my readers and asked where they place themselves on a scale from liberal to conservative, the results might look like this:
To be clear, I have not done a survey and this is fake data, but if it were real, we would conclude that my audience is more liberal, on average, than the general population. So in the normal use of the word skew, we might say that this distribution “skews to the left”.
But according to statisticians, that would be wrong, because within the field of statistics, skew has been given a technical meaning that is contrary to its normal use. Here’s how Wikipedia explains the technical definition:
positive skew: The right tail is longer; the mass of the distribution is concentrated on the left of the figure. The distribution is said to be right-skewed, right-tailed, or skewed to the right, despite the fact that the curve itself appears to be skewed or leaning to the left; right instead refers to the right tail being drawn out and, often, the mean being skewed to the right of a typical center of the data. A right-skewed distribution usually appears as a left-leaning curve.
https://en.wikipedia.org/wiki/Skewness
By this definition, we would say that the distribution of political alignment in my audience is “skewed to the right”. It is regrettable that the term was defined this way, because it’s very confusing.
Recently I ran a Twitter poll to see what people think skew means. Here are the results:
Interpreting these results is almost paradoxical: the first two responses are almost equally common, which proves that the third response is correct. If the statistically-literate people who follow me on Twitter don’t agree about what skew means, we have to treat it as ambiguous unless specified.
The comments suggest I’m not the only one who thinks the technical definition is contrary to intuition.
This has always been confusing for me, since the shape of a right-skewed distribution looks like it’s “leaning” to the left…
I learnt it as B, but there’s always this moment when I consciously have to avoid thinking it’s A.
This is one of those things where once I learned B was right, I hated it so much that I never forgot it.
It gets worse
If you think the definition of skew is bad, let’s talk about bias. In the context of statistics, bias is “a systematic tendency which causes differences between results and fact”. In particular, sampling bias is bias caused by a non-representative sampling process.
In my imaginary survey, the mean of the sample is less than the actual mean in the population, so we could say that my sample is biased to the left. Which means that the distribution is technically biased to the left and skewed to the right. Which is particularly confusing because in natural use, bias and skew mean the same thing.
So 20th century statisticians took two English words that are (nearly) synonyms, and gave them technical definitions that can be polar opposites. The result is 100 years of confusion.
For early statisticians, it seems like creating confusing vocabulary was a hobby. In addition to bias and skew, here’s a partial list of English words that are either synonyms or closely related, which have been given technical meanings that are opposites or subtly different.
accuracy and precision
probability and likelihood
efficacy and effectiveness
sensitivity and specificity
confidence and credibility
And don’t get me started on “significance”.
If you got this far, it seems like you are part of my audience, so if you want to answer a one-question survey about your political alignment, follow this link. Thank you!
My poll and this article were prompted by this excellent video about the Central Limit Theorem:
Around the 7:52 mark, a distribution that leans left is described as “skewed towards the left”. In statistics jargon, that’s technically incorrect, but in this context I think it’s is likely to be understood as intended.
This recent article in the Washington Post reports that “a police department in Maryland is training officers to spot the signs of driving high by watching people toke up in a tent”. The story features Lt. John O’Brien, who is described as a “trained drug recognition expert”. It also quotes a defense attorney who says, “There are real questions about the scientific validity of what they’re doing.”
As it happens, the scientific validity of Drug Recognition Experts is one of the examples in my forthcoming book, Probably Overthinking It. The following is an excerpt from Chapter 9: “Fairness and Fallacy”.
In September 2017 the American Civil Liberties Union (ACLU) filed suit against Cobb County, Georgia on behalf of four drivers who were arrested for driving under the influence of cannabis. All four were evaluated by Officer Tracy Carroll, who had been trained as a “Drug Recognition Expert” (DRE) as part of a program developed by the Los Angeles Police Department in the 1970s.
At the time of their arrest, all four insisted that they had not smoked or ingested any cannabis products, and when their blood was tested, all four results were negative; that is, the blood tests found no evidence of recent cannabis use.
In each case, prosecutors dismissed the charges related to impaired driving. Nevertheless, the arrests were disruptive and costly, and the plaintiffs were left with a permanent and public arrest record.
At issue in the case is the assertion by the ACLU that, “Much of the DRE protocol has never been rigorously and independently validated.”
So I investigated that claim. What I found was a collection of studies that are, across the board, deeply flawed. Every one of them features at least one methodological error so blatant it would be embarrassing at a middle school science fair.
As an example, the lab study most often cited to show that the DRE protocol is valid was conducted at Johns Hopkins University School of Medicine in 1985. It concludes, “Overall, in 98.7% of instances of judged intoxication the subject had received some active drug”. In other words, in the cases where one of the Drug Recognition Experts believed that a subject was under the influence, they were right 98.7% of the time.
That sounds impressive, but there are several problems with this study. The biggest is that the subjects were all “normal, healthy” male volunteers between 18 and 35 years old, who were screened and “trained on the psychomotor tasks and subjective effect questionnaires used in the study”.
By design, the study excluded women, anyone older than 35, and anyone in poor health. Then the screening excluded anyone who had any difficulty passing a sobriety test while they were sober — for example, anyone with shaky hands, poor coordination, or poor balance.
But those are exactly the people most likely to be falsely accused. How can you estimate the number of false positives if you exclude from the study everyone likely to yield a false positive? You can’t.
Another frequently-cited study reports that “When DREs claimed drugs other than alcohol were present, they [the drugs] were almost always detected in the blood (94% of the time)”. Again, that sounds impressive until you look at the methodology.
Subjects in this study had already been arrested because they were suspected of driving while impaired, most often because they had failed a field sobriety test.
Then, while they were in custody, they were evaluated by a DRE, that is, a different officer trained in the drug evaluation procedure. If the DRE thought that the suspect was under the influence of a drug, the suspect was asked to consent to a blood test; otherwise they were released.
Of 219 suspects, 18 were released after a DRE performed a “cursory examination” and concluded that there was no evidence of drug impairment.
The remaining 201 suspects were asked for a blood sample. Of those, 22 refused and 6 provided a urine sample only.
Of the 173 blood samples, 162 were found to contain a drug other than alcohol. That’s about 94%, which is the statistic they reported.
But the base rate in this study is extraordinarily high, because it includes only cases that were suspected by the arresting officer and then confirmed by the DRE. With a few generous assumptions, I estimate that the base rate is 86%; in reality, it was probably higher.
To estimate the base rate, let’s assume:
All 18 of the suspects who were released were, in fact, not under the influence of a drug, and
The 28 suspects who refused a blood test were impaired at the same rate as the 173 who agreed, 94%.
Both of these assumptions are generous; that is, they probably overestimate the accuracy of the DREs. Even so, they imply that 188 out of 219 blood tests would have been positive, if they had been tested. That’s a base rate of 86%.
Because the suspects who were released were not tested, there is no way to estimate the sensitivity of the test, but let’s assume it’s 99%, so if a suspect is under the influence of a drug, there is a 99% chance a DRE would detect it. In reality, it is probably lower.
With these generous assumptions, we can use the following table to estimate the sensitivity of the DRE protocol.
Suspects
Prob Positive
Cases
Percent
Impaired
86
0.99
85.14
93.8
Not impaired
14
0.40
5.60
6.2
With 86% base rate, we expect 86 impaired suspects out of 100, and 14 unimpaired. With 99% sensitivity, we expect the DRE to detect about 85 true positives. And with 60% specificity, we expect the DRE to wrongly accuse 5.6 suspects. Out of 91 positive tests, 85 would be correct; that’s about 94%, as reported in the study.
But this accuracy is only possible because the base rate in the study is so high. Remember that most of the subjects had been arrested because they had failed a field sobriety test. Then they were tested by a DRE, who was effectively offering a second opinion.
But that’s not what happened when Officer Tracy Carroll arrested Katelyn Ebner, Princess Mbamara, Ayokunle Oriyomi, and Brittany Penwell. In each of those cases, the driver was stopped for driving erratically, which is evidence of possible impairment. But when Officer Carroll began his evaluation, that was the only evidence of impairment.
So the relevant base rate is not 86%, as in the study; it is the fraction of erratic drivers who are under the influence of drugs. And there are many other reasons for erratic driving, including distraction, sleepiness, and the influence of alcohol. It’s hard to say which explanation is most common. I’m sure it depends on time and location. But as an example, let’s suppose it is 50%; the following table shows the results with this base rate.
Suspects
Prob Positive
Cases
Percent
Impaired
50
0.99
49.5
71.2
Not impaired
50
0.40
20.0
28.8
With 50% base rate, 99% sensitivity, and 60% specificity, the predictive value of the test is only 71%; under these assumptions, almost 30% of the accused would be innocent. In fact, the base rate, sensitivity, and specificity are probably lower, which means that the value of the test is even worse.
The suit filed by the ACLU was not successful. The court decided that the arrests were valid because the results of the field sobriety tests constituted “probable cause” for an arrest. As a result, the court did not consider the evidence for, or against, the validity of the DRE protocol. The ACLU has appealed the decision.
I’m happy to report that copyediting of Modeling and Simulation in Python is done, and the book is off to the printer! Electronic versions are available now from No Starch Press; print copies will be available in May, but you can pre-order now from No Starch Press, Amazon, and Barnes and Noble.
To celebrate, I just published one of the case studies from the end of Part I, which is about simulating discrete systems. The case study explores a classic question from queueing theory:
Suppose you are designing the checkout area for a new store. There is room for two checkout counters and a waiting area for customers. You can make two lines, one for each counter, or one line that serves both counters.
In theory, you might expect a single line to be better, but it has some practical drawbacks: in order to maintain a single line, you would have to install rope barriers, and customers might be put off by what seems to be a longer line, even if it moves faster.
So you’d like to check whether the single line is really better and by how much.
Simulation can help answer this question. The following figure shows the three scenarios I simulated:
The leftmost diagram shows a single queue (with customers arriving at rate 𝜆) and a single server (with customers completing service at rate 𝜇).
The center diagram shows a single queue with two servers, and the rightmost diagram shows two queue with two servers.
This figure shows the time customers are in the system, including wait time and service time, as a function of the arrival rate. The orange line shows the average we expect based on analysis; the blue dots show the result of simulations.
This comparison shows that the simulation and analysis are consistent. It also demonstrates one of the features of simulation: it is easy to quantify not just the average we expect but also the variation around the average.
That capability turns out to be useful for this problem because, as it turns out, the difference between the one-queue and two-queue scenarios is small compared to the variation, which suggests the advantage would be unnoticed in practice.
I conclude:
The two configurations are equally good as long as both servers are busy; the only time two lines is worse is if one queue is empty and the other contains more than one customer. In real life, if we allow customers to change lanes, that disadvantage can be eliminated.
From a theoretical point of view, one line is better. From a practical point of view, the difference is small and can be mitigated. So the best choice depends on practical considerations.
On the other hand, you can do substantially better with an express line for customers with short service times. But that’s a topic for another case study.
Sadly, today is my last day at DrivenData, so it’s a good time to review one of the projects I’ve been working on, using probabilistic predictions from Zamba to find animals in camera trap videos, like this:
Most recently, I’ve been working on calibrating the predictions from convolutional neural networks (CNNs). I haven’t written about it, but at ODSC East 2023, I’m giving a talk about it:
Don’t Discard That Pangolin, Calibrate Your Deep Learning Classifier
Suppose you are an ecologist studying a rare species like a pangolin. You can use motion-triggered camera traps to collect data about the presence and abundance of species in the wild, but for every video showing a pangolin, you have 100 that show other species, and 100 more that are blank. You might have to watch hours of video to find one pangolin.
Deep learning can help. Project Zamba provides models that classify camera trap videos and identify the species that appear in them. Of course, the results are not perfect, but we can often remove 80% of the videos we don’t want while losing only 10-20% of the videos we want.
But there’s a problem. The output from deep learning classifiers is generally a “confidence score”, not a probability. If a classifier assigns a label with 80% confidence, that doesn’t mean there is an 80% chance it is correct. However, with a modest number of human-generated labels, we can often calibrate the output to produce more accurate probabilities, and make better predictions.
In this talk, I’ll present use cases based on data from Africa, Guam, and New Zealand, and show how we can use deep learning and calibration to save the pangolin… or at least the pangolin videos. This real-world problem shows how users of ML models can tune the results to improve performance on their applications.
The ODSC schedule isn’t posted yet, but I’ll fill in the details later.
Way back in 2013, I wrote this blog post explaining why you should never use a statistical test to check whether a sample came from a Gaussian distribution. I argued that data from the real world never come from a Gaussian distribution, or any other simple mathematical model, so the answer to the question is always no. And there are only two possible outcomes from the test:
If you have enough data, the test will reject the hypothesis that the data came from a Gaussian distribution, or
If you don’t have enough data, the test will fail to reject the hypothesis.
Either way, the result doesn’t tell you anything useful.
In this article, I will explore a particular example and demonstrate this relationship between the sample size and the outcome of the test. And I will conclude, again, that
Choosing a distribution is not a statistical question; it is a modeling decision. No statistical test can tell you whether a particular distribution is a good model for your data.
I’ll start by generating a sample that is actually from a lognormal distribution, then use the sample mean and standard deviation to make a Gaussian model. Here’s what the empirical distribution of the sample looks like compared to the CDF of the Gaussian distribution.
It looks like the Gaussian distribution is a pretty good model for the data, and probably good enough for most purposes.
According to the Anderson-Darling test, the test statistic is 1.7, which exceeds the critical value, 0.77, so at the 5% significance level, we can reject the hypothesis that this sample came from a Gaussian distribution. That’s the right answer, so it might seem like we’ve done something useful. But we haven’t.
Sample size
The result from the A-D test depends on the sample size. The following figure shows the probability of rejecting the null hypothesis as a function of sample size, using the lognormal distribution from the previous section.
When the sample size is more than 200, the probability of rejection is high. When the sample size is less than 100, the probability of rejection is low. But notice that it doesn’t go all the way to zero, because there is always a 5% chance of a false positive.
The critical value is about 120; at that sample size, the probability of rejecting the null is close to 50%.
So, again, if you have enough data, you’ll reject the null; otherwise you probably won’t. Either way, you learn nothing about the question you really care about, which is whether the Gaussian model is a good enough model of the data for your purposes.
That’s a modeling decision, and no statistical test can help. In the original article, I suggested some methods that might.
A recent question on Reddit asked about using resampling with logistic regression. The responses suggest two ways to do it, one parametric and one non-parametric. I implemented both of them and then invented a third, which is hybrid of the two.
Different ways of computing sampling distributions – and the statistics derived from them, like standard errors and confidence intervals – yield different results. None of them are right or wrong; rather, they are based on different modeling assumptions.
In this example, it is easy to implement multiple models and compare the results. If they were substantially different, we would need to think more carefully about the modeling assumptions they are based on and choose the one we think is the best description of the data-generating process.
But in this example, the differences are small enough that they probably don’t matter in practice. So we are free to choose whichever is easiest to implement, or fastest to compute, or convenient in some other way.
It is a common error to presume that the result of an analytic method is uniquely correct, and that results from computational methods like resampling are approximations to it. Analytic methods are often fast to compute, but they are always based on modeling assumptions and often based on approximations, so they are no more correct than computational methods.
The Centers for Disease Control and Prevention states on its website that “in the United States, cigarette smoking causes about 90% of lung cancers.” If S is the event “smokes cigarettes” and L is the event “has lung cancer,” then the probability 0.90 is expressed in probability notation as
P(S and L).
P(S | L).
P(L | S).
Let’s consider a scenario that’s not exactly what the question asks about, but will help us understand the relationships among these quantities:
Suppose 20% of people smoke, so out of 1000 people, we have 200 smokers and 800 nonsmokers.
Suppose 1% of nonsmokers get lung cancer, so out of 800 nonsmokers, there would be 8 cases of lung cancer, 0 caused by smoking.
And suppose 20% of smokers get lung cancer, 19% caused by smoking and 1% caused by something else (same as the nonsmokers). Out of 200 smokers, there would be 40 cases of lung cancer, 38 caused by smoking.
In this scenario, there are a total of 48 cases of lung cancer, 38 caused by smoking. So smoking caused 38/48 cancers, which is 79%.
P(S and L) is 40 / 1000, which is 4%.
P(S | L) = 40 / 48, which is 83%.
P(L | S) = 40 / 200, which is 20%.
From this scenario we can conclude:
The percentage of cases caused by smoking does not correspond to any of the listed probabilities, so the answer to the question is “None of the above”.
In order to compute these quantities, we need to know the percentage of smokers and the risk ratio for smokers vs nonsmokers.
In reality, the relationships among these quantities are complicated by time: the percentage of smokers changes over time, and there is a long lag between smoking and cancer diagnosis.
How can I find the least frequent value (antimode) between 2 modes in a bimodal distribution?
I’m only mildly self taught in anything in statistics so please be patient with my ignorance. I’ve found so little info on a Google search for “antimode” that I thought it was a word made up by the author of a random article.
Then I found one tiny mention on the Wikipedia page for “Multimodal distribution” but no citation or details beyond that it’s the least frequent value between modes.
What data do I need in order to find this number and what is the formula?
This site had a short mention of it and in their example listed:
Mode A: 33.25836
Mode B: 71.55446
Antimode: 55.06092
But I can’t seem to reverse engineer it with just this data.
Here’s the reply I wrote:
With continuous data, there is no off-the-shelf formula to compute modes or antimodes. You have to make some modeling decisions.
One option is to use kernel density estimation KDE. Adjust the parameters until, in your judgment, the result is a good representation of the distribution, and then you can read off maxima and minima.
“Tell me if you agree or disagree with this statement: Most men are better suited emotionally for politics than are most women.”
That’s one of the questions on the General Social Survey. In 1974, when it was first asked, 47% of respondents said they agreed. In 2018, the most recent survey where this question was asked, that percentage was down to 14%.
Last week I posed the same question to an audience at PyData NYC. Of 56 people who responded, none of them agreed with the statement (0%).
However, if we take the question literally, the correct answer is “I don’t know”. In fact, it is almost certain that either
Most men are better suited emotionally for politics than are most women, or
Most women are better suited emotionally for politics than are most men.
The reason is that saying that “most A are better than most B” is equivalent to saying that the median in group A is higher than the median in group B.
Here’s the proof: If one median is higher than the other, we can find a value between them where most members of group A are higher (by the definition of median) and most members of group B are lower. So unless the medians are identical, either “most A are better than most B” or the other way around.
And what does it mean to be “suited emotionally for politics”? I don’t know, but if we had to quantify it, it would be some combination of emotional factors. Men and women differ in many ways, on average; the difference is small for some factors and bigger for others. So it is likely that a combination of factors also differs between the groups. And even if the difference is very small, the medians are unlikely to be identical.
So, if we take the question literally, the best answer is “I don’t know”. But it’s clear that people don’t take the question literally. Instead, they are probably answering something like “Do you think men are so much better suited emotionally for politics that they should hold a large majority of positions of power?”
