conditional probability of dependent random variables

Question

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)

Answer 1

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]

Limit[p, δ -> 0]

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

Answer 2

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.

Answer 3

@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

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