How to compute conditional probabilities?

How can I use Mathematica to derive the conditional probability of a given multivariate PDF?

For instance, if we have a bivariate normal:

dist = BinormalDistribution[{mu1, mu2}, {sigma1, sigma2}, r]


How can I ask Mathematica to derive $P[x_1|x_2]$?

(For this simple case, the solution is $\mathcal{N}\left(\mu_1+\frac{\sigma_1}{\sigma_2}\rho( x_2 - \mu_2),\, (1-\rho^2)\sigma_1^2\right)$; see https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Bivariate_case)

• I know very little about statistics and am a little bit too tired to read the wiki, but have you seen, e.g. EstimatedDistribution or the suite of statistical tools? You can also search through the other guides to see if something there is more likely to be what you need. Mar 16, 2017 at 15:38
• EstimatedDistribution is aimed at working with empirical statistics (fitting models to data). My problem is not about data, but about deriving a closed-form starting from other closed forms (finding a conditional probability is akin to deriving the formula for a "slice" of a higher dimensional function). Mathematica has some capabilities to work with probability distributions in closed form (e.g., Probability, TransformedDistribution), but none of these seems to do what I need. Mar 16, 2017 at 15:47
• Alas I am out of my depth. Hopefully someone else can come in and give you a hand. You can also try cross-posting on Wolfram Community. I think a WRI statistician might be more likely to see that. Just be sure to link your posts. Mar 16, 2017 at 15:50

Here's a brute force way: Use the definition of continuous conditional density function:

dist12 = BinormalDistribution[{μ1, μ2}, {σ1, σ2}, ρ]
dist2 = MarginalDistribution[dist12, 2]
conditionalDensity = PDF[dist12, {x1, x2}]/PDF[dist2, x2]


$$\frac{\exp \left(\frac{(\text{x2}-\text{\mu 2})^2}{2 \text{\sigma 2}^2}-\frac{\frac{(\text{x1}-\text{\mu 1})^2}{\text{\sigma 1}^2}-\frac{2 \rho (\text{x1}-\text{\mu 1}) (\text{x2}-\text{\mu 2})}{\text{\sigma 1} \text{\sigma 2}}+\frac{(\text{x2}-\text{\mu 2})^2}{\text{\sigma 2}^2}}{2 \left(1-\rho ^2\right)}\right)}{\sqrt{2 \pi } \sqrt{1-\rho ^2} \text{\sigma 1}}$$

Now...in this case will Mathematica automatically recognize this as a Normal distribution? Maybe not:

d1 = ProbabilityDistribution[conditionalDensity, {x1, -Infinity, Infinity}]


Update: An alternative approach

Here's another way to get the conditional density. First determine the conditional distribution function and then differentiate to get the conditional probability density function:

cdf = Probability[y1 <= x1 \[Conditioned] y2 == x2,
{y1, y2} \[Distributed] BinormalDistribution[{μ1, μ2}, {σ1, σ2}, ρ]]
pdf = D[cdf, x1]


• I would think that you should use dist2 = MarginalDistribution[ dist12, 2] because you want to integrate over the first var to get the marginal distribution of the second var, which is x2.
– gwr
Mar 16, 2017 at 18:00
• @gwr Sorry about that. You are correct. (I clearly couldn't the manual very well. Somehow I got it into my head the the indices given were the ones that were integrated out - rather than the ones that were kept.) I've now edited the answer. Thanks!
– JimB
Mar 16, 2017 at 18:36
• Something does not look right here, as the derived result does not match the result I linked to in wikipedia. That is, conditionalDensity should be identical to PDF[NormalDistribution[\[Mu]1 + \[Sigma]1/\[Sigma]2 \[Rho] (x2 - \ \[Mu]2), (1 - \[Rho]^2) \[Sigma]1^2], x1], and it is not. Mar 16, 2017 at 21:41
• That's because you have the variance placed in NormalDistribution[] and need to use the standard deviation. In other words, Mathematica uses $N(\mu,\sigma)$ rather than $N(\mu,\sigma^2)$ to define a normal distribution. Then my result is uglier but equivalent to the wiki answer.
– JimB
Mar 16, 2017 at 23:36
• @JimBaldwin Yes, you are right. Many thanks! Note for Wolfram Research: Mathematica could provide a ConditionalDistribution function, similarly to how it provides a MarginalDistribution one. This function should know the typical results (e.g., for the multivariate normal) and evaluate the definition (i.e., $f_Y(y \mid X=x) = \frac{f_{X, Y}(x, y)}{f_X(x)}$) for all other cases. Mar 17, 2017 at 2:31

Unfortunately David G. Stork has deleted his answer; I would nevertheless pick up his answer which imo can be "cured" rather easily.

I would believe that the following holds:

\begin{align} p(x_1 | x_2) \propto p(x1,x2|x2) \end{align}

which would mean that, as Dr. Stork has written, the joint probability density $p(x_1, x_2)$ (here the bivariate normal density) evaluated for a variable $x_1$ given a fixed value for $x_2$ will be equal to the conditional probability density $p(x_1|x_2)$ up to a normalizing constant.

We can see this from adding the option Method -> "Normalize" to the former answer given by Dr. Stork:

condDist = ProbabilityDistribution[
PDF[ BinormalDistribution[ {μ1, μ2}, {σ1, σ2}, ρ] ], {x1, x2} ],
{ x1, -∞, +∞ },
Method -> "Normalize"
]


(Note that there are some warnings which likely call for some additional assumptions we should make using the option Assumptions, but this should not worry us here too much.)

Comparing this expression with the answer obtained by @JimBaldwin will show that it is completely equivalent up to some expression which does contain neither $x_1$ nor $x_2$ and thus is a mere constant.

We can use this propotionality for example if we factor a complicated joint probability distribution (for example a Bayesian Network) into a multiplicative term of conditional and marginal probabilities. Here we may use the joint probability densities in the way indicated without the costly normalization, as we could do that only once at the end of our inference.