Bayesian modelling of two sequential tests

You need to parametrize those prior and conditional distributions (that is find a distribution and parameters that convey the prior and test information in a mathematical format). Assuming X is binary it might look something like:

with pm.Model() as model:
    X = pm.Beta("X", alpha=1, beta=10)  # 10% prior probability of X = 1
    testA = pm.Bernoulli(p=fA(X), observed=testA_result)
    testB = pm.Categorical(p=fB(X, testA), observed=testB_result))

You need to find out what the right prior information is for X, and what fA, and fB should be, those are the functions that give you the probability of the outcomes given the variable of interest X and the observed results testA (in the case of testB). That is something you have to figure out.

For example, perhaps the probability of testA being positive is well enough described by a logistic transformation of X, so:

def fA(x):
    return pm.math.sigmoid(x)

You also need to find how fB should convert A, and X into categorical probabilities.

Once you have that, you just have to call:

with model:
    posterior = pm.sample()

Which will give you the posterior probability of X, given your test results.

The hard work is in those fA and fB which have to reflect the information that the tests A and B give you in a mathematical format. That’s the modelling task, which can’t be automated. I suggest you try to start with just X and testA, and once you understand that you can add testB.

2 Likes