Calculating probabilities symbolically

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}}*)

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]