# conditional probability of dependent random variables

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)

• 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. – Bob Hanlon Aug 5 '19 at 14:23

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.

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.

• Would you elaborate on how your approach gets $Pr(X=1|Y=y)$? (Even for a particular value of $y$?) – JimB Aug 5 '19 at 15:44
• It doesn't. I was wrong. I have attempted to atone by editing my earlier incorrect effort. Again, sorry for sowing confusion. – Seth Chandler Aug 5 '19 at 17:40
• 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). – JimB Aug 5 '19 at 18:12

@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)}$$

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