How do you work out Conditional Probabilities in Mathematica. Is it possible?
Asked Answered
U

2

13

Can Mathematica do Bayes Rule conditional probability calculations, without doing the calculation manually? If so how?

I have been searching both the Mathemtaica doco and the web for a hint but cannot find anything. I am not after how to do Bayes Rule manually via Mathematica, I want to know if there is a way to define the conditional probabilities and calculate other ones automagically.

So to use the toy example assuming Bernoulli distributions

P(Cancer+) = 0.01
P(Cancer-) = 0.99

P(Test+|Cancer+) = 0.9
P(Test-|Cancer+) = 0.1
P(Test+|Cancer-) = 0.2
P(Test-|Cancer-) = 0.8

Is it possible to work out

P(Cancer+|Test+) = 0.0434

So using the below.

Print["P(C+) = ", PCancerT=BernoulliDistribution[0.01]];
Print["P(C-) = ", PCancerF=BernoulliDistribution[0.99]];
Print[]
Print["P(T+|C+) = ", PTestTGivenCancerT=BernoulliDistribution[0.9]];
Print["P(T-|C+) = ", PTestFGivenCancerT=BernoulliDistribution[0.1]];
Print["P(T+|C-) = ", PTestTGivenCancerF=BernoulliDistribution[0.2]];
Print["P(T-|C-) = ", PTestFGivenCancerF=BernoulliDistribution[0.8]];
Print[]
Print["P(T+,C+) = ", PTestTAndCancerT = Probability[vCT&&vTTCT,{vCT\[Distributed]PCancerT,vTTCT\[Distributed]PTestTGivenCancerT}]];
Print["P(T-,C+) = ", PTestFAndCancerT = Probability[vCT&&vTFCF,{vCT\[Distributed]PCancerT,vTFCF\[Distributed]PTestFGivenCancerT}]];
Print["P(T+,C-) = ", PTestTAndCancerF = Probability[vCF&&vTTCF,{vCF\[Distributed]PCancerF,vTTCF\[Distributed]PTestTGivenCancerF}]];
Print["P(T-,C-) = ", PTestFAndCancerF = Probability[vCF&&vTTCF,{vCF\[Distributed]PCancerF,vTTCF\[Distributed]PTestFGivenCancerF}]];
Print[]
Print["P(C+|T+) = ?"];
Print["P(C+|T-) = ?"];
Print["P(C-|T+) = ?"];
Print["P(C-|T-) = ?"];

I can work out the joint probabilities by defining all the probability tables manually, but is there a way to get Mathematica to do the heavy lifting? Is there a way to define and calculate these kind of conditional probabilities?

Many thanks for any assistance, even it its “You can’t... stop trying” :)

PS : was this an attempt at doing something along these lines? Symbolic Conditional Expectation in Mathematica

Uphroe answered 25/10, 2011 at 9:54 Comment(0)
O
13

Actually... I worked this out symbolically in the past, and it covers a lot of simple (unchained) probabilities. I guess it wouldn't be that hard to add chaining(see below). You're welcome to reply with augmentation. The symbolic approach is far more flexible than working with Bernoulli distributions and creating a proc for Bayes theorem and thinking about the right way to apply it every time.

NOTE: The functions are not bound, like in the post above ((0 < pC < 1) && (0 < pTC < 1) && (0 < pTNC < 1)) because sometimes you want "unweighted" results, which produce numbers outside of 0-1 range, then you can bring back into the range by dividing by some normalizing probability or product of probabilities. If you do want to add bounds for error checking, do this:
P[A_ /;0<=A<=1] := some_function_of_A;

use Esc+cond+Esc to enter \\[Conditioned] symbol in Mathematica.

Remove[P];
Unprotect@Intersection;
Intersection[A_Symbol, B_Symbol] := {A, B}
Intersection[A_Not, B_Symbol] := {A, B}
Intersection[A_Symbol, B_Not] := {A, B}
P[Int_List/; Length@Int == 2] := P[Int[[2]] \[Conditioned] Int[[1]]] P[Int[[1]]]
   (*//  P(B) given knowledge of P(A)  //*)
P[B_, A_] := If[NumericQ@B, B, 
                P[B \[Conditioned] A] P[A] + P[B \[Conditioned] Not@A] P[Not@A]]
P[Not@B_, A_: 1] := If[NumericQ@A, 1 - P[B], 1 - P[B, A]]
P[A_ \[Conditioned] B_] := P[A \[Intersection] B]/P[B, A]
P[Not@A_ \[Conditioned] B_] := 1 - P[A \[Conditioned] B];

You then use it as such:

P[Cancer]=0.01;

Don't need "not cancer" since P[!Cancer] yields 0.99 (Esc+not+Esc types a very pretty logical not symbol, but Not[A], !A or \[Not]A work just fine too)

P[Test \[Conditioned] Cancer] = 0.9
P[Test \[Conditioned] ! Cancer] = 0.2

again: P[!Test \\[Conditioned] Cancer] will be 1-P[Test \\[Conditioned] Cancer] by definition, unless you override it.

Now let's query this model:

P[Test, Cancer]
P[!Test, Cancer]

returns

0.207
0.793

and

P[Cancer \[Conditioned] Test]
P[!Cancer \[Conditioned] Test]
P[Cancer \[Conditioned] !Test]
P[!Cancer \[Conditioned] !Test]

returns

0.0434783
0.956522
0.00126103
0.998739

I guess it would be a nice idea to define P(B|A1,A2,A3,...,An), anyone up for coding the chain rule using NestList or something like it? I didn't need it for my project, but it wouldn't be that difficult to add, should someone need it.

Octangle answered 4/12, 2011 at 20:23 Comment(9)
Thanks. I had another account here before, but can't remember what it was anymore. Staring to build reputation from ground up. Thanks for your encouragement @Mr.WizardOctangle
My apologies if I should have remembered your name. If you find the old account, you should be able to get an account merge by flagging one of your posts for diamond-moderator attention and requesting it.Toothy
It was Gr3gK1 stackoverflow.com/users/959803/gr3gk1 See if you can do something, I'd really appreciate it!Octangle
Gregory, you can do it yourself I believe, unless new-user privileges have changed. You should see a link in gray just below your post above that reads "flag" -- click it, and choose "it needs ♦ moderator attention" then "other." In the box, request an account merge and include "users/959803/gr3gk1" as above. If you have trouble I will do this for you, but it makes more sense coming from you.Toothy
Requested! Thanks for your help! You're a real asset to this site!Octangle
@GregoryKlopper You might want to take a stab at this question #8282725Henricks
Wow! Nice. Sorry for the delay.Uphroe
This is awesome. I'm surprised however that Mathematica doesn't have something like this built in?Myrtie
Pretty old thread but I'm at this point too. Have used your solution @GregoryKlopper and it works particularly well for basic conditions like P(A|B). I'm completely new to Mathematica and would like to know how to extend it to P(A|B,C,...) and larger CPTs (trying to use it with a given Bayes Network). Any idea?Fibroid
K
5

I wouldn't complicate the issue with Print statements and BernoulliDistributions. You know the probabilities, so the simplest thing to do is to calculate them directly, but perhaps using vectors to get P(B), and using the fact that pr(cancer) = 1-pr(not cancer) and so on.

Bayes' Theorem states that P(A|B)=(P(A ⋂ B))/(P(B))

The intersection is calculated as the conditional probability (test given cancer) times the probability of cancer.

So something like the following should work:

conditionalProb[pC_, pTC_, pTNC_] /; 
 (0 < pC < 1) && (0 < pTC < 1) && (0 < pTNC < 1) :=
 (pTC * pC)/({pTC, pTNC}.{pC, 1 - pC})

conditionalProb[0.01, 0.9, 0.2]

0.0434783

And yes, the Probability functionality in version 8 does allow you to calculate conditional probabilities "automagically", but for a problem like this with Bernoulli-distributed events, it's overkill.

Kurr answered 25/10, 2011 at 11:26 Comment(4)
Thankyou, the Print was just an attempt hopefully to make clearer when executed in mma to get across what I was to trying to achieve, maybe it backfired. As as the answer I know I can easily work out this particular answer with a simple Bayws rule function. I guess what I was asking does mma have the facility to do more complex inference without having to explicitly define the calculation. So for example being able to work out P(Cancer-|Test+, Test-, Test-). Where the probability of Cancer+/- changes slightly with each consecutive test result.Uphroe
Also possibly using different distributions.Uphroe
Yes I am aware about the probability function but have not been able to work out how to get it do the above calculation. I may be missing something very fundamental.Uphroe
@Bart, the Probability function should do what you want for other distributions, but the more direct approach in my function is clearer for the case you have of two events where the values are in essence Boolean. Changing the probability of cancer works nicely with my function: just change the first parameter, e.g. conditionalProb[0.03,0.9,0.2] .Kurr

© 2022 - 2024 — McMap. All rights reserved.