Binomial and negative binomial distributions

Today's post is prompted by this question from Reddit:

How do I calculate the distribution of the number of selections (with replacement) I need to make before obtaining k? For example, let's say I am picking marbles from a bag with replacement. There is a 10% chance of green and 90% of black. I want k=5 green marbles. What is the distribution number of times I need to take a marble before getting 5?

I believe this is a geometric distribution. I see how to calculate the cumulative probability given n picks, but I would like to generalize it so that for any value of k (number of marbles I want), I can tell you the mean, 10% and 90% probability for the number of times I need to pick from it.

Another way of saying this is, how many times do I need to pull on a slot machine before it pays out given that each pull is independent?

Note: I've changed the notation in the question to be consistent with convention.

In [1]:

from __future__ import print_function, division

import thinkplot
from thinkstats2 import Pmf, Cdf

from scipy import stats
from scipy import special

%matplotlib inline

Solution¶

There are two ways to solve this problem. One is to relate the desired distribution to the binomial distribution.

If the probability of success on every trial is p, the probability of getting the kth success on the nth trial is

PMF(n; k, p) = BinomialPMF(k-1; n-1, p) p

That is, the probability of getting k-1 successes in n-1 trials, times the probability of getting the kth success on the nth trial.

Here's a function that computes it:

In [2]:

def MakePmfUsingBinom(k, p, high=100):
    pmf = Pmf()
    for n in range(1, high):
        pmf[n] = stats.binom.pmf(k-1, n-1, p) * p
    return pmf

And here's an example using the parameters in the question.

In [3]:

pmf = MakePmfUsingBinom(5, 0.1, 200)
thinkplot.Pdf(pmf)

We can solve the same problem using the negative binomial distribution, but it requires some translation from the parameters of the problem to the conventional parameters of the binomial distribution.

The negative binomial PMF is the probability of getting r non-terminal events before the kth terminal event. (I am using "terminal event" instead of "success" and "non-terminal" event instead of "failure" because in the context of the negative binomial distribution, the use of "success" and "failure" is often reversed.)

If n is the total number of events, n = k + r, so

r = n - k

If the probability of a terminal event on every trial is p, the probability of getting the kth terminal event on the nth trial is

PMF(n; k, p) = NegativeBinomialPMF(n-k; k, p) p

That is, the probability of n-k non-terminal events on the way to getting the kth terminal event.

Here's a function that computes it:

In [4]:

def MakePmfUsingNbinom(k, p, high=100):
    pmf = Pmf()
    for n in range(1, high):
        r = n-k
        pmf[n] = stats.nbinom.pmf(r, k, p)
    return pmf

Here's the same example:

In [5]:

pmf2 = MakePmfUsingNbinom(5, 0.1, 200)
thinkplot.Pdf(pmf2)

And confirmation that the results are the same within floating point error.

In [6]:

diffs = [abs(pmf[n] - pmf2[n]) for n in pmf]
max(diffs)

Out[6]:

8.6736173798840355e-17

Using the PMF, we can compute the mean and standard deviation:

In [7]:

pmf.Mean(), pmf.Std()

Out[7]:

(49.998064403376738, 21.207570382894403)

To compute percentiles, we can convert to a CDF (which computes the cumulative sum of the PMF)

In [8]:

cdf = Cdf(pmf)
scale = thinkplot.Cdf(cdf)

And here are the 10th and 90th percentiles.

In [9]:

cdf.Percentile(10), cdf.Percentile(90)

Out[9]:

(26, 78)

Copyright 2016 Allen Downey

MIT License: http://opensource.org/licenses/MIT

Probability is hard: part four

This is the fourth part of a series of posts about conditional probability and Bayesian statistics.

• In the first article, I presented the Red Dice problem, which is a warm-up problem that might help us make sense of the other problems.
• In the second article, I presented the problem of interpreting medical tests when there is uncertainty about parameters like sensitivity and specificity.  I explained that there are (at least) four models of the scenario, and I hinted that some of them yield different answers.  I presented my solution for one of the scenarios, using the computational framework from Think Bayes.
• In the third article, I presented my solution for Scenario B.

Now, on to Scenarios C and D.  For the last time, here are the scenarios:

Scenario A: The patients are drawn at random from the relevant population, and the reason we are uncertain about t is that either (1) there are two versions of the test, with different false positive rates, and we don't know which test was used, or (2) there are two groups of people, the false positive rate is different for different groups, and we don't know which group the patient is in.

Scenario B: As in Scenario A, the patients are drawn at random from the relevant population, but the reason we are uncertain about t is that previous studies of the test have been contradictory. That is, there is only one version of the test, and we have reason to believe that t is the same for all groups, but we are not sure what the correct value of t is.

Scenario C: As in Scenario A, there are two versions of the test or two groups of people. But now the patients are being filtered so we only see the patients who tested positive and we don't know how many patients tested negative. For example, suppose you are a specialist and patients are only referred to you after they test positive.

Scenario D: As in Scenario B, we have reason to think that t is the same for all patients, and as in Scenario C, we only see patients who test positive and don't know how many tested negative.

And the questions we would like to answer are:

1. What is the probability that a patient who tests positive is actually sick?
2. Given two patients who have tested positive, what is the probability that both are sick?

I have posted a new Jupyter notebook with a solution for Scenario C and D.  It also includes the previous solutions for Scenarios A and B; if you have already seen the first two, you can skip to the new stuff.

You can read a static version of the notebook here.

OR

You can run the notebook on Binder.

If you click the Binder link, you should get a home page showing the notebooks in the repository for this blog.  Click on test_scenario_cd.ipynb to load the notebook for this article.

Comments

As I said in the first article, this problem turned out to be harder than I expected.  First, it took time, and a lot of untangling, to identify the different scenarios.  Once we did that, figuring out the answers was easier, but still not easy.

For me, writing a simulation for each scenario turned out to be very helpful.  It has the obvious benefit of providing an approximate answer, but it also forced me to specify the scenarios precisely, which is crucial for problems like these.

Now, having solved four versions of this problem, what guidance can we provide for interpreting medical tests in practice?  And how realistic are these scenarios, anyway?

Most medical tests are more sensitive than 90% and yield false positives less often than 20-40% of the time.  But many symptoms that indicate disease are neither sensitive nor specific.  So this analysis might be more applicable to interpreting symptoms, rather than test results.

I think all four scenarios are relevant to practice:

• For some signs of disease, we have estimates of true and false positive rates, but when there are multiple inconsistent estimates, we might be unsure which estimate is right.  For some other tests, we might have reason to believe that these parameters are different for different groups of patients; for example, many antibody tests are more likely to generate false positives in patients who have been exposed to other diseases.
• In some medical environments, we could reasonably treat patients as a random sample of the population, but more often patients have been filtered by a selection process that depends on their symptoms and test results.  We considered two extremes: in Scenarios A and B, we have a random sample, and in C and D we only see patients who test positive.  But in practice the selection process is probably somewhere in between.

If the selection process is unknown, and we don't know whether the parameters of the test are likely to be different for different groups, one option is to run the analysis for all four scenarios (or rather, for the three scenarios that are different), and use the range of the results to characterize uncertainty due to model selection.

UPDATES

[May 11, 2016] In this comment on Reddit, /u/ppcsf shows how to solve this problem using a probabilistic programming language based on C#.  The details are in this gist.

Probability is hard: part three

This is the third part of a series of posts about conditional probability and Bayesian statistics.

• In the first article, I presented the Red Dice problem, which is a warm-up problem that might help us make sense of the other problems.
  
• In the second article, I presented the problem of interpreting medical tests when there is uncertainty about parameters like sensitivity and specificity.  I explained that there are (at least) four models of the scenario, and I hinted that some of them yield different answers.  I presented my solution for one of the scenarios, using the computational framework from Think Bayes.

Now, on to Scenario B.  For those of you coming in late, here are the Scenarios again:

Scenario A: The patients are drawn at random from the relevant population, and the reason we are uncertain about t is that either (1) there are two versions of the test, with different false positive rates, and we don't know which test was used, or (2) there are two groups of people, the false positive rate is different for different groups, and we don't know which group the patient is in.

Scenario B: As in Scenario A, the patients are drawn at random from the relevant population, but the reason we are uncertain about t is that previous studies of the test have been contradictory. That is, there is only one version of the test, and we have reason to believe that t is the same for all groups, but we are not sure what the correct value of t is.

Scenario C: As in Scenario A, there are two versions of the test or two groups of people. But now the patients are being filtered so we only see the patients who tested positive and we don't know how many patients tested negative. For example, suppose you are a specialist and patients are only referred to you after they test positive.

Scenario D: As in Scenario B, we have reason to think that t is the same for all patients, and as in Scenario C, we only see patients who test positive and don't know how many tested negative.

And the questions we would like to answer are:

1. What is the probability that a patient who tests positive is actually sick?
2. Given two patients who have tested positive, what is the probability that both are sick?

I have posted a new Jupyter notebook with a Solution for Scenario B.  It also includes the previous solution for Scenario A; if you read it before, you can skip to the new stuff.

You can read a static version of the notebook here.

OR

You can run the notebook on Binder.

If you click the Binder link, you should get a home page showing the notebooks and Python modules in the repository for this blog.  Click on test_scenario_b.ipynb to load the notebook for this article.

I'll give you a chance to think about Scenarios C and D, and I'll post my solution in the next couple of days.

Enjoy!

Probability is hard, part two

If you read yesterday's post, you know that my colleague Sanjoy Mahajan and I have been working on a series of problems related to conditional probability and Bayesian statistics.  In the previous article, I presented the Red Dice problem, which is relatively simple.  I posted it here because it presents four different versions of the problem, which we'll need for the next step, and it demonstrates the computational tools I use in my solution.

Now I'll present the problem we are really working on: interpreting medical tests.  Suppose we test a patient to see if they have a disease, and the test comes back positive. What is the probability that the patient is actually sick (that is, has the disease)?

To answer this question, we need to know:

• The prevalence of the disease in the population the patient is from. Let's assume the patient is identified as a member of a population where the known prevalence is p.
• The sensitivity of the test, s, which is the probability of a positive test if the patient is sick.
• The false positive rate of the test, t, which is the probability of a positive test if the patient is not sick.

Given these parameters, it is straightforward to compute the probability that the patient is sick, given a positive test.  This problem is a staple in probability classes.

But here's the wrinkle.  What if you are uncertain about one of the parameters: p, s, or t?  For example, suppose you have reason to believe that t is either 0.2 or 0.4 with equal probability.  And suppose p=0.1 and s=0.9.

The questions we would like to answer are:

1. What is the probability that a patient who tests positive is actually sick?
2. Given two patients who have tested positive, what is the probability that both are sick?

As we did with the Red Dice problem, we'll consider four scenarios:

Scenario A: The patients are drawn at random from the relevant population, and the reason we are uncertain about t is that either (1) there are two versions of the test, with different false positive rates, and we don't know which test was used, or (2) there are two groups of people, the false positive rate is different for different groups, and we don't know which group the patient is in.

Scenario B: As in Scenario A, the patients are drawn at random from the relevant population, but the reason we are uncertain about t is that previous studies of the test have been contradictory. That is, there is only one version of the test, and we have reason to believe that t is the same for all groups, but we are not sure what the correct value of t is.

Scenario C: As in Scenario A, there are two versions of the test or two groups of people. But now the patients are being filtered so we only see the patients who tested positive and we don't know how many patients tested negative. For example, suppose you are a specialist and patients are only referred to you after they test positive.

Scenario D: As in Scenario B, we have reason to think that t is the same for all patients, and as in Scenario C, we only see patients who test positive and don't know how many tested negative.

I have posted a Jupyter notebook with solutions for the basic scenario (where we are certain about t) and for Scenario A.  You might want to write your own solution before you read mine.  And then you might want to work on Scenarios B, C, and D.

You can read a static version of the notebook here.

OR

You can run the notebook on Binder.

If you click the Binder link, you should get a home page showing the notebooks and Python modules in the repository for this blog.  Click on test_scenario_a.ipynb to load the notebook for this article.

I'll post my solution for Scenarios B, C, and D next week.

Probability is hard

For more than a month, my colleague Sanjoy Mahajan and I have been banging our heads on a series of problems related to conditional probability and Bayesian statistics.  We knew when we started that this material is tricky, as demonstrated by veridical paradoxes like the Monty Hall problem, the Girl Named Florida, and so on.  But even though we were prepared, we have been surprised, continually, by how long it is taking and how effectively we have confused ourselves and each other.

A few times we have hit a brick wall on a hard problem and made a strategic retreat by working on a simpler problem.  At this point, we have retreated all the way to what I'll call the Red Dice problem, which goes like this:

Suppose I have a six-sided die that is red on 2 sides and blue on 4 sides, and another die that's the other way around, red on 4 sides and blue on 2. 

I choose a die at random and roll it, and I tell you it came up red. What is the probability that I rolled the second die (red on 4 sides)?  And if I do it again, what's the probability that I get red again?

There are several variations on this problem, with answers that are subtly different.  I explain the variations, and my solution, in a Jupyter notebook:

You can read a static version of the notebook here.

OR

You can run the notebook on Binder.

If you click the Binder link, you should get a home page showing the notebooks and Python modules in the repository for this blog.  Click on red_dice.ipynb to load the notebook for this article.

Once we have settled the Red Dice problem, we will get back to the original problem, which relates to interpreting medical tests (a classic example of Bayesian inference that, again, turns out to be more complicated than we thought).


UPDATE:  After reading the notebook, some readers are annoyed with me because Scenarios C and D are not consistent with the way I posed the question.  I'm sorry if you feel tricked -- that was not the point!  To clarify, Scenarios A and B are legitimate interpretations of the question, as posed, which is deliberately ambiguous.  Scenarios C and D are exploring a different version because it will be useful when we get to the next problem.

Stats can't make modeling decisions

Here's a question that appeared recently on the Reddit statistics forum:

If effect sizes of coefficient are really small, can you interpret as no relationship?  Coefficients are very significant, which is expected with my large dataset. But coefficients are tiny (0.0000001). Can I conclude no relationship? Or must I say there is a relationship, but it's not practical?

I posted a response there, but since we get questions like this a lot, I will write a more detailed response here.

First, as several people mentioned on Reddit, you have to distinguish between a small coefficient and a small effect size.  The size of the coefficient depends on the units it is expressed in.  For example, in a previous article I wrote about the relationship between a baby's birth weight and its mother's age ("Are first babies more likely to be light?").  With weights in pounds and ages in years, the estimated coefficient is about 0.017 pounds per year.

At first glance, that looks like a small effect size.  But the average birth weight in the U.S. is about 7.3 pounds, and the range from the youngest to the oldest mother was more than 20 years.  So if we say the effect size is "about 3 ounces per decade", that would be easier to interpret.  Or it might be even better express the effect in terms of percentages; for example, "A 10 year increase in mother's age is associated with a 2.4% increase in birth weight."

So that's the first part of my answer:

Expressing effect size in practical terms makes it easier to evaluate its importance in practice.

The second part of my answer addresses the question, "Can I conclude no relationship?"  This is a question about modeling, not statistics, and

Statistical analysis can inform modeling choices, but it can't make decisions for you.

As a reminder, when you make a model of a real-world scenario, you have to decide what to include and what to leave out. If you include the most important things and leave out less important things, your model will be good enough for most purposes.

But in most scenarios, there is no single uniquely correct model. Rather, there are many possible models that might be good enough, or not, for various purposes.

Based on statistics alone, you can't say whether there is, or is not, a relationship between two variables.  But statistics can help you justify your decision to include a relationship in a model or ignore it.

The affirmative

If you want to argue that an effect SHOULD be included in a model, you can justify that decision (using classical statistics) in two steps: 

1) Estimate the effect size and use background knowledge to make an argument about why it matters in practical terms.  For example, a 3 ounce difference in birth weight might be associated with real differences in health outcomes, or not.

AND

2) Show that the p-value is small, which at least suggests that the observed effect is unlikely to be due to chance.  (Some people will object to this interpretation of p-values, but I explain why I think it is valid in "Hypothesis testing is only mostly useless").

In my study of birth weight, I argued that mother's age should be included in the model because the effect size was big enough to matter in the real world, and because the p-value was very small.

The negative

If you want to argue that it is ok to leave an effect out of a model, you can justify that decision in one of two ways:

1) If you apply a hypothesis test and get a small p-value, you probably can't dismiss it as random.  But if the estimated effect size is small, can use background information to make an argument about why it is negligible.

OR

2)  If you apply a hypothesis test and get a large p-value, that suggests that the effect you observed could be explained by chance.  But that doesn't mean the effect is necessarily negligible.  To make that argument, you need to consider the power of the test.  One way to do that is to find the smallest hypothetical effect size that would yield a high probability of a significant test.  Then you can say something like, "If the effect size were as big as X, this test would have a 90% of being statistically significant. The test was not statistically significant, so the effect size is likely to be less than X.  And in practical terms, X is negligible."

The Bayesian

So far I have been using the logic of classical statistics, which is problematic in many ways.

Alternatively, in a Bayesian framework, the result would be a posterior distribution on the effect size, which you could use to generate an ensemble of models with different effect sizes. To make predictions, you would generate predictive distributions that represent your uncertainty about the effect size. In that case there's no need to make binary decisions about whether there is, or is not, a relationship.

Or you could use Bayesian model comparison, but I think that a mostly misguided effort to shoehorn Bayesian methods into a classical framework.  But that's a topic for another time.

Q&A

Update May 4, 2016: I got a few questions about this article that I thought I should answer here.

Q: Doesn't N matter too? More likely to find significance in larger samples, making justification that much more important.

A: For the arguments I outlined, we don't need to know N directly.  You are right that if N is very large, an effect might be statistically significant even if it is very small.  But then you could apply Negative Argument #1.

If N is very small, an effect might not be statistically significant even if it is substantial.  In that case you wouldn't be able to make a strong argument either way.  The affirmative argument would fail because the apparent effect could plausibly be explained by chance.  The negative argument would fail because the test was underpowered (specifically, in Negative Argument #2, X would be big).

Q: For the posterior to be useful for decision making, you need to know that the model is causally correct as well, don't you?

A: Good question!  It depends on what kind of decision-making you are talking about.

For example, suppose you find that preschool education predicts future earnings.  The effect might be directly causal, or it might be that children who get preschool education have other advantages.

If the task is to predict future earnings, you would probably want to include preschool education in the model, and it would probably help, causal or not.

But suppose you are considering an intervention, like universal preschool education.  In that case, you definitely want to know whether the effect is causal.  If it isn't, the intervention might do little or no good.


Update May 5, 2016.  In response to another question, I wrote, "For any (non-trivial) real-world scenario, there is no one unique correct model; rather, there are many models that might be good enough (or not) for various purposes. Model choice can be informed by quantitative factors, but there might be several contradictory factors, as well as value judgments."

Someone asked me to give examples of contradictory criteria and value judgments.  Here's my reply:

Let's stick with predicting future earnings, and let's say there are about 10 predictive factors you are considering, like SAT scores, high school grades, parent's socioeconomic status, etc.

With just 10 factors, there are more than 1000 models to choose from. For each candidate model, you might consider these criteria:

1) How good the predictions are.  This one is obviously, but there are several ways to define it, depending on whether you want to minimize absolute error, relative error, mean squared error, or some other cost function.

2) How many factors are in the model. You might prefer a simpler model, but there are different ways you might define "simple".

3) How early different factors can be measured. If elementary school grades predict almost as well as high school grades, you might prefer elementary school grades because they are available earlier.

4) How easily different factors can be measured. You might prefer a model that runs on cheap data, even if its not quite as good as a more expensive model.

5) How interpretable the model is. If you are trying to explain something about the factors that contribute to earnings, you might prefer a model that makes sense to people (although you might make a bad choice if you let your preconceptions drive the bus).

6) How causal the model is. If you have background knowledge about which factors are more likely to be causal, you might want to focus on those factors, depending on the purpose of the model.

And I could go on. But I don't think I'm saying anything truly profound, just that there is no objective, uniquely correct way to navigate tradeoffs like this.

(If you like, this is basically what Kuhn said in "Objectivity, Value Judgment, and Theory Choice". He talked about theory choice rather than model choice, but I think that's the same thing.)

Bayes on Jupyter

I am working on an updated version of my workshop, Bayesian Statistics Made Simple, now using Jupyter notebooks (formerly known as IPython). It's still a work in progress, but you can see a draft of my slides here:

 

If you want to run the code, you can run the notebook in your browser by hitting this button



You should see a home page like this:

If you want to try the exercises, open workshop01.ipynb.  If you just want to see the answers, open workshop01_soln.ipynb.

Either way, you should be able to run the notebooks in your browser and try out the examples.

If you run into any problems, let me know.  Comments and suggestions are welcome.

Special thanks to the generous people who run Binder, which makes it easy to share and reproduce computation.  You can watch their video here: