covidsim.py

import numpy.random as npr

# Monte Carlo simulation code to estimate the probability that I had COVID
# back in February 2020. My flatmate had COVID almost certainly (he still doesn't
# have his sense of smell back and also got a positive antibody test about 
# 6 months later). I had more ambiguous symptoms but was ill with something
# COVID like at the same time. How likely was it that I had COVID vs something
# unrelated?

# Please note that I am not a specialist in disease transmission and it is
# entirely plausible that I have made serious errors in the course of writing
# this. ALso please note that even if I/you had COVID one should not assume
# that results in immunity, especially a year on.

# Simulation paramters.

# SIZE is number of samples to take. Quite high due to considering relatively rare events.
SIZE = 10 ** 8

# Random seed to use. Doesn't matter very much just fixing for consistency.
SEED = 0

# Henceforth some slightly made up probabilities. They're probably wrong but
# I've tried to pick numbers that are at least justifiable and where uncertain
# have tried to err on the side that makes it less likely that I had COVID.

# These numbers might be totally off it's hard to say but the relative chances
# seems to be what's important and I think it's fair to say that M was an order
# of magnitude more likely than me to catch it (given that he's a doctor and
# was also seeing a lot more people)
M_CAUGHT_OUTSIDE = 0.001
I_CAUGHT_OUTSIDE = 0.0001

# This is based on https://thetab.com/uk/2020/10/07/this-is-how-likely-you-are-to-catch-coronavirus-at-uni-if-your-flat-mate-tests-positive-177632
# which is probably wrong but if anything based on the rate at which M and I
# gave each other whatever we happened to catch it's probably too low if anything.
TRANSMISSION = 0.38

# IDK people's guesses for this seem all over the place but this doesn't seem
# unreasonable.
SYMPTOMATIC_GIVEN_COVID = 0.8

# I got ill a lot in that flat, so getting ill in roughly that time period is
# pretty plausible.
I_CAUGHT_SOMETHING_ELSE = 0.3

# This is based on severity of symptoms I had: i am ill that badly about once a
# year at most - My illness patterns tend to be little and often - so it's at
# least surprising to be that ill from something that wasn't COVID.
SYMPTOMATIC_GIVEN_OTHER = 0.1

def bernoulli(p):
    """Generate samples which are each true with probability ``p``"""
    return npr.random(SIZE) <= p


if __name__ == '__main__':
    npr.seed(SEED)

    m_caught_covid_outside = bernoulli(M_CAUGHT_OUTSIDE)
    i_caught_covid_outside = bernoulli(I_CAUGHT_OUTSIDE)

    # The event that whichever of us caught COVID passed it to the other.    
    transmit = bernoulli(TRANSMISSION)

    # I had covid if I got it outside or if M had it and passed it on to me
    i_had_covid = i_caught_covid_outside | (m_caught_covid_outside & transmit)

    # Same as above but the other way around.
    m_had_covid = m_caught_covid_outside | (i_caught_covid_outside & transmit)

    # Maybe I had the flu or something? Who knows?
    i_had_something_else = bernoulli(I_CAUGHT_SOMETHING_ELSE)

    # Either I had COVID and got COVID symptoms from it, or I got something
    # else and had COVID symptoms from it. Possibly both, which would be bad
    # luck!
    i_had_covid_symptoms = (i_had_covid & bernoulli(SYMPTOMATIC_GIVEN_COVID)) | (i_had_something_else & bernoulli(SYMPTOMATIC_GIVEN_OTHER))

    # We know within rounding error of certainty that M had COVID (positive
    # antibody result and persistent anosmia). Also I had the observed COVID
    # symptoms. So this is the event that we are conditioning on.
    real_life = i_had_covid_symptoms & m_had_covid

    # In the observed scenario, how often did I have COVID?
    results = i_had_covid[real_life]

    # Weight the estimate towards zero (this is an improper prior but honestly at these sample sizes it doesn't really matter what prior I pick)
    print(f"In {len(results)} times I had COVID {results.sum()} times, for an estimated chance of {(results.sum()) / (len(results) + 1)}")