The Canonical Bayesian cancer question answered in PYMC

[zweli]

I came accross the canonical Bayesian cancer question and thought it would be fun to solve it using PYMC.

Here is the question:

And here is the solution in PYMC:

prevalence_cancer = 0.008
sensitivity = 0.98
specificity = 0.97

with pm.Model() as model:

    cancer = pm.Bernoulli("cancer", prevalence_cancer)
    
    test_positive = pm.Bernoulli("test_positive", 
                                 p=pm.math.switch(cancer, sensitivity, 1 - specificity), 
                                 observed=1
                                 ) 
    
    trace = pm.sample(draws=5000, tune=1000, chains=4, random_seed=1)

let’s look at the posterior:

or the summary:

This is what the solution would be if you did this by hand. I think this was a nice simple example for newbies to PYMC

[ricardoV94]

Yeah although for this I would just do grid analysis

[bob-carpenter]

That’s a nice simulation, but this one’s easy to do analytically—it’s why it’s the first example. Let z = 1 if the patient has cancer and let y = 1 if they get a positive test. Then

\begin{align} \Pr[z = 1 | y = 1] &= \dfrac{\Pr[z = 1, y = 1]}{\Pr[y = 1]} \\[8pt] &= \dfrac{\Pr[z=1, y=1]}{\Pr[z=1, y=1] + \Pr[z=0, y=1]} \\[8pt] &= \dfrac{\Pr[z=1] \cdot \Pr[y=1 | z=1]}{\Pr[z=1] \cdot \Pr[y=1 \mid z=1] + \Pr[z=0] \cdot \Pr[y=1 | z = 0]} \\[8pt] &= \frac{0.008 \cdot 0.98} {0.008 \cdot 0.98 + (1 - 0.008) \cdot (1 - 0.97)} \\[4pt] &= 0.21 \end{align}

First line is just definition of conditional probability, denominator in second from total law of probability, the third line is just the chain rule, and the rest is just plugging in numbers.