20
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ What do you mean with p(a|!b) ? I've never seen this notation. $\endgroup$ – eldo Jul 3 '14 at 21:18
  • 1
    $\begingroup$ @eldo, I mean p(a| not b) $\endgroup$ – kmace Jul 3 '14 at 23:52
  • $\begingroup$ Like two dices NOT showing the same number - could you give an example ? $\endgroup$ – eldo Jul 3 '14 at 23:59
  • 1
    $\begingroup$ The only example I ever encountered of a symbolic evaluation of Bayes equation is by Gregory Klopper in this link. $\endgroup$ – Romke Bontekoe Jul 4 '14 at 7:27
  • $\begingroup$ thanks @RomkeBontekoe, I've seen that too. Its a great example, but its not as general as it could be. I'm going to assume that Mathematica cannot do such a calculation out of the box. $\endgroup$ – kmace Jul 9 '14 at 0:22
5
$\begingroup$

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}}*)
$\endgroup$
  • $\begingroup$ You are awesome $\endgroup$ – kmace Jun 1 '16 at 20:14
  • $\begingroup$ How could this be extended to include multiple arguments in the probability? For example, in my case I have values like P(A, B | C) and P(D | E, F) (the commas are the intersection / "and" operator) $\endgroup$ – Pro Q Nov 30 '18 at 6:15
2
$\begingroup$

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]
    

    enter image description here

  2. Independent:

    full[0.6, 0.3, 0.6]
    

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.