Calculating probabilities symbolically

Question

Is there a way to solve for statistical quantities analytically/symbolically in Mathematica?

example 1:
Lets say that I want to do a calculation that requires Bayes theorem.

I know p(a), p(b) and p(a | b). I want to calculate p(b | a).

Now, for this simple example, I know the solution, so I could write the relationship myself, but what about more complicated examples? Surely Mathematica is able to solve probability problems like it can solve algebraic ones? I want to stay away from simulation. I want Mathematica to do the five pages of math that I normally do on paper. Can Mathematica do this for me?

I'd like to type something like

Solve[p(a) == a1 && p(b) == a2 && p(a|b) == a3, p(b|a)]

and get

p(b|a) -> (a3*a2)/a1

example 2:
I want to use the law of total probability:

I know p(a|b), p(a|!b), and p(b) and I want to calculate p(a).

I'd like to type something like

Solve[p(a|b) == a1 && p(a|!b) == a2 && p(b) == a3, p(a)]

and get

p(a) -> a1*a3 + a2*(1-a3)

Answer 1

Here is code for a function bSolve which is to be used essentially like Solve. Given some equations as its first argument and some variables (or in fact more complicated expressions) as its second argument, it finds the values of the variables, eliminating all occurrences of p[...]. The probability of a is denoted by p[a] and the conditional probability of a given b is denoted p[a, b]. It wouldn't be hard to have support for \[Conditioned] if you want.

This code works by writing all probabilities in terms of the 2^n joint probabilities p[a&&b&&...], p[!a&&b&&...], p[a&&!b&&...] etc. (where n is the number of different events). Then we write all relevant equations, eliminate the p[...] and are left with the values of g[...] requested by the user.

p::usage = 
  "In bSolve, p[a] denotes the probability of a, and p[a, b] denotes \
the conditional probability of a given b." ;
bSolve::usage = 
  "bSolve[eqs, vars] tries to solve a system of equations eqs on \
probabilities p[a] or conditional probabilities p[a, b] \
(conventionally denoted p(a|b)) for the variables vars.";
Begin["bSolve`"];
Unprotect[bSolve, p];
Clear[bSolve, p];
bSolve[eqs : Except[_List], goal_] := bSolve[{eqs}, goal];
bSolve[eqs_, goal : Except[_List]] := bSolve[eqs, {goal}];
bSolve[eqs_List, goal_List] := (
   goalvars = g /@ Range[Length[goal]];
   alleqs = Join[eqs, Thread[goalvars == goal]];
   argsofp = Cases[alleqs, p[a___] :> a, Infinity];
   vars = 
    List @@ Union@(d @@ argsofp //.
          d[a___, (Not | And | Or | Xor)[b___], c___] :> d[a, b, c]);
   qtab = Table[q[i - 1], {i, 2^Length[vars]}];
   ptab = 
    p /@ Flatten@Outer[And, Sequence @@ Table[{a, Not[a]}, {a, vars}]];
   sols = 
    Solve[Eliminate[
      Join[{Total[qtab] == 1}, 
       alleqs /. {p[x_, y_] :> p[x && y]/p[y]} /. 
        p[x_] :> (qtab.Boole@BooleanTable[x, vars])], qtab], goalvars];
   Thread[goal -> goalvars] /. sols
   );
SyntaxInformation[p] = {"ArgumentsPattern" -> {_, _.}};
SyntaxInformation[bSolve] = {"ArgumentsPattern" -> {_, _}};
Protect[bSolve, p];
End[];

Let's test:

bSolve[{p[a] == a1, p[b] == a2, p[a, b] == a3}, {p[b, a], p[! b], 
  p[a && b], p[Xor[a, b]]}]

(*=> {{p[b, a] -> (a2 a3)/a1, p[! b] -> 1 - a2, p[a && b] -> a2 a3, 
  p[a \[Xor] b] -> a1 + a2 - 2 a2 a3}}*)

bSolve[p[a, b] == a1 && p[a, ! b] == a2 && p[b] == a3, p[a]]

(*=> {{p[a] -> a2 + a1 a3 - a2 a3}}*)

and an example from an 2011 Stackoverflow question on a similar topic:

bSolve[p[Cancer] == 0.01 && p[Test, Cancer] == 0.9 && 
  p[Test, ! Cancer] == 0.2, p[Cancer, Test]]

(*=> {{p[Cancer, Test] -> 0.0434783}}*)

Answer 2

You can create your own discrete probability distributions and then use built-in functions. The following uses this approach if the following are known: p(a),p(b},p(a|b). It can be adapted for other scenario. I present it as motivation.

pr[pa_, pb_, pagb_] := Module[{s},
  s = First@
    Solve[{p00 + p10 + p01 + p11 == 1, p11 == pagb pb, 
      p10 + p11 == pa, p01 + p11 == pb}, {p00, p01, p10, p11}];
  ProbabilityDistribution[Piecewise[{
     {p00 /. s, x == 0 && y == 0},
     {p01 /. s, x == 0 && y == 1},
     {p10 /. s, x == 1 && y == 0},
     {p11 /. s, x == 1 && y == 1}}],
   {x, 0, 1, 1}, {y, 0, 1, 1}]
  ]
full[u_, v_, w_] := Module[{pd, t, r, w1, w2, w3, w4, w5},
  pd = pr[u, v, w];
  t = Tuples[{0, 1}, 2];
  r = Thread[t -> {"Not a,Not b", "Not a, b", "a,Not b", "a,b"}];
  w1 = {# /. r, 
      Probability[{x, y} == #, {x, y} \[Distributed] pd]} & /@ t;
  w2 = {StringForm["`1` given `2`", 
       Sequence @@ (StringSplit[{#1, #2} /. r, ","])], 
      Probability[
       x == #1 \[Conditioned] y == #2, {x, y} \[Distributed] pd]} & @@@
     t;
  w3 = {StringForm["`2` given `1`", 
       Sequence @@ (StringSplit[{#1, #2} /. r, ","])], 
      Probability[
       y == #2 \[Conditioned] x == #1, {x, y} \[Distributed] pd]} & @@@
     t;
  w4 = {"a", Probability[x == 1, {x, y} \[Distributed] pd]};
  w5 = {"b", Probability[y == 1, {x, y} \[Distributed] pd]};
  Row[{Framed[
     TableForm[w1~Join~w2~Join~w3~Join~{w4}~Join~{w5}, 
      TableHeadings -> {None, {"Event", "Probability"}}]],
    DiscretePlot3D[PDF[pd, {x, y}], {x, {0, 1}}, {y, {0, 1}}, 
     ExtentSize -> Full, Ticks -> {{0, 1}, {0, 1}, {0, 1}}, 
     AxesLabel -> {"a", "b", "Probability"}, ImageSize -> 300]}, 
   Frame -> True]]

Examples:

1. Case with dependency:

   full[0.5, 0.3, 0.25]
   

   

2. Independent:

   full[0.6, 0.3, 0.6]