4
$\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$
  • 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 '19 at 14:23
10
$\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.

|improve this answer|||||
$\endgroup$
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.

|improve this answer|||||
$\endgroup$
  • 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 '19 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$ – Seth Chandler Aug 5 '19 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 '19 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)}$$

|improve this answer|||||
$\endgroup$
  • 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 '19 at 20:15
  • $\begingroup$ I'll update it, thanks $\endgroup$ – Conor Cosnett Aug 5 '19 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 '19 at 20:18

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.