5
$\begingroup$

Suppose I have 3 random variables:

$$X \sim \mbox{Bernoulli}(1/2)$$ $$Z \sim \mbox{Normal}(0,1)$$ $$Y = X+Z$$

How do I compute the conditional probability density:

$$P(X=1 | Y=y)$$

Attempt1:

Probability[ X == 1 \[Conditioned] X + Z == y, 
           {
            X \[Distributed] BernoulliDistribution[1/2]
           ,Z \[Distributed] NormalDistribution[]
           }
         ]

Attempt2:

D[Probability[ X == 1 \[Conditioned] X + Z >= y, 
           {
            X \[Distributed] BernoulliDistribution[1/2]
           ,Z \[Distributed] NormalDistribution[]
           }
         ],y]

Attempt3:

Likelihood[
      TransformedDistribution[X + Z, 
                   {
                    X \[Distributed]BernoulliDistribution[1/2], 
                   Z \[Distributed] NormalDistribution[]}]
           , {y}]

Pencil and Paper attempt:

$$P(X=1 | Y=y) = \frac{P(X=1 , Y=y)}{P(Y=y)}$$ $$= \frac{P(X=1 , X+Z=y)}{P(Y=y)}$$ $$= \frac{P(X=1)P(Z=y-1)}{P(Y=y)}$$ $$= \frac{P(X=1)P(Z=y-1)}{P(X=1)P(Z=y-1)+P(X=0)P(Z=y-0)}$$


$$P(Z=y)=\frac{e^{-\frac{y^2}{2}}}{\sqrt{2 \pi }}$$ $$P(Z=y-0)=\frac{e^{-\frac{y^2}{2}}}{\sqrt{2 \pi }}$$ $$P(Z=y-1)=\frac{e^{-\frac{1}{2} (y-1)^2}}{\sqrt{2 \pi }}$$ $$P(X=1)=\frac{1}{2}$$ $$P(X=0)=\frac{1}{2}$$


$$P(X=1 | Y=y) = \frac{e^{-\frac{1}{2} (y-1)^2}}{2 \sqrt{2 \pi } \left(\frac{e^{-\frac{y^2}{2}}}{2 \sqrt{2 \pi }}+\frac{e^{-\frac{1}{2} (y-1)^2}}{2 \sqrt{2 \pi }}\right)}$$

$$P(X=1|Y=y) = \frac{e^y}{e^y+\sqrt{e}}$$

(please excuse my sloppy notation, please interpret $P(x)$ as "probability density" where appropriate)

$\endgroup$
1
  • 4
    $\begingroup$ The probability that a random variable with a continuous distribution (e.g., normal) takes a specific value is zero. The PDF is a probability density and a probability only arises when integrating over an interval. $\endgroup$
    – Bob Hanlon
    Aug 5, 2019 at 14:23

3 Answers 3

11
$\begingroup$

As @BobHanlon points out that the probability is zero for a specific value. But probabilities are not (necessarily) 0 in intervals. So we can get a probability statement for an interval and then take the limit as the size of that interval goes to zero.

p = Probability[X == 1 \[Conditioned] Abs[X + Z - y] <= δ, 
  {X \[Distributed] BernoulliDistribution[1/2], Z \[Distributed] NormalDistribution[]},
  Assumptions -> δ > 0 && y ∈ Reals]

Probability associated with an interval around y

Limit[p, δ -> 0]

Limit of above probability statement as delta goes to zero

So we end up with what you did with pencil and paper.

$\endgroup$
0
3
$\begingroup$

WARNING: The answer I provided earlier is wrong because, perhaps among other reasons, it equates a PDF with a probability. I provide a better answer immediately below that slightly extends @JimB's correct answer.

  prob = Limit[
  FullSimplify[
   Probability[
    X == 1 \[Conditioned] 
     Abs[X + Z - y] < \[Delta], {X \[Distributed] 
      BernoulliDistribution[f], 
     Z \[Distributed] NormalDistribution[m, s]}, 
    Assumptions -> {y \[Element] Reals, \[Delta] > 0}]], \[Delta] -> 
   0]

The output is:

 (1 - (E^((1 + 2*m - 2*y)/(2*s^2))*(-1 + f))/f)^(-1)

If we now do the following substitution, we get the special case considered in the original post.

prob /. {f -> 1/2, m -> 0, s -> 1}

We get the same answer as the OP did using pencil and paper and @JimB did using his code. It is in a slightly different form, but if you test for equality, you get True.

Also, if anyone wants a numeric answer, I get the following:

prob /. {f -> 1/2, m -> 0, s -> 1, y -> 1} // N

0.622459


The following earlier answer is not correct. Sorry for sowing confusion.

Let's try something harder and then do your problem as a special case:

    ydist[f_, m_, s_] := ProbabilityDistribution[{"CDF", 
    FullSimplify[
    Probability[
     X == 1 \[Conditioned] 
      X + Z >= y, {X \[Distributed] BernoulliDistribution[f], 
      Z \[Distributed] 
       NormalDistribution[m, s]}]]}, {y, -\[Infinity], \[Infinity]}]

You can then write

 ydist[f,m,s]

and get an elaborate formula. Or you can plug in the values from your example:

 ydist[1/2,0,1]

From there you can compute a PDF and bunch of other stuff. I would note that the answer produced by this method is more complicated than your paper and pencil version. So, either you have erred, I have erred, or the Wolfram engine has erred. Hmm.

$\endgroup$
3
  • 1
    $\begingroup$ Would you elaborate on how your approach gets $Pr(X=1|Y=y)$? (Even for a particular value of $y$?) $\endgroup$
    – JimB
    Aug 5, 2019 at 15:44
  • $\begingroup$ It doesn't. I was wrong. I have attempted to atone by editing my earlier incorrect effort. Again, sorry for sowing confusion. $\endgroup$ Aug 5, 2019 at 17:40
  • 2
    $\begingroup$ No worries. Join the club. You might have seen one of my recent mess-ups. But the good thing is that there are lots of folks on this forum who will not hesitate to point out my mistakes. I think that's a good thing (although I don't like it at all when I do mess up). $\endgroup$
    – JimB
    Aug 5, 2019 at 18:12
1
$\begingroup$

@Seth Chandler's original approach succeeds when calculating the inverse conditional probability. It is easy to apply Bayes theorem to this expression to get my desired result, thanks.

pYgX = PDF[ProbabilityDistribution[{"SF", 

Probability[  

         X + Z >= y \[Conditioned] X == 1

       , {
         X \[Distributed] BernoulliDistribution[1/2]
       , Z \[Distributed] NormalDistribution[]
         }

                       ]} ,{y, -\[Infinity], \[Infinity]}] , y];

pX = 1/2;


pY = PDF[TransformedDistribution[X + Z, {X \[Distributed] BernoulliDistribution[1/2], Z \[Distributed] NormalDistribution[]}] , y];

pXgY = (pYgX* pX)/pY // Simplify

enter image description here

$$P(X=1 | Y=y) = \frac{P(Y=y| X=1)P(X=1)}{P(Y=y)}$$

$\endgroup$
3
  • 1
    $\begingroup$ If you use X + Z <= y in place of X + Z >= y and Simplify instead of FullSimplify, you get the desired answer. $\endgroup$
    – JimB
    Aug 5, 2019 at 20:15
  • $\begingroup$ I'll update it, thanks $\endgroup$
    – Conor
    Aug 5, 2019 at 20:17
  • 1
    $\begingroup$ Or you could change CDF to SF and FullSimplify to Simplify (because you did give it a survival function rather than a cdf. $\endgroup$
    – JimB
    Aug 5, 2019 at 20:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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