Information flow using entropy does not give me an output

Here are few equations to calculate the quantity of information flow (i..e., mutual information): I wrote the following code (a primitive one) to calculate x(21):

ClearAll[dist, prX1aX2, prX2, prX1gX2, H1, H12, x21];

dist = MultinormalDistribution[{Subscript[μ, 1], Subscript[μ,
2]}, {{Subscript[σ,1]^2, ρ Subscript[σ, 1] Subscript[σ, 2]}, {ρ Subscript[σ, 1] Subscript[σ, 2],
Subscript[σ, 2]^2}}];

prX1aX2 =
PDF[dist, {x1, x2}];  (* joint probability distribution: p(x1,x2) *)

prX2 = PDF[MarginalDistribution[dist, 2], x2];  (* marginal probability distribution of x2: p(x2) *)

prX1gX2 = (prX1aX2/prX2);  (* conditional density of x1 given x2 *)

(* entropy *)
H1 = Expectation[-prX1*Log[prX1], x1 \[Distributed]
NormalDistribution[Subscript[μ, 1], Subscript[σ, 1]]]

H12 = Expectation[-prX1gX2*Log[prX1gX2], {x1 \[Distributed]
NormalDistribution[Subscript[μ, 1], Subscript[σ, 1]],
x2 \[Distributed] NormalDistribution[Subscript[μ, 2], Subscript[σ, 2]]}]

x21 = H1 - H12;

This 'code' does not give me any output but keeps running. Where do I do wrong? I need help to find out why it is not working? Also, if possible, make this 'code' more elegant...

• I think there are two main issues. (1) You're multiplying by the associated density in $H1$ and $H12$ when Expectation already includes that. And (2) H12 treats x2 as if it were random: you've conditioned on a specific value so it's not a random variable. In fact I think the $dx_2$ in H12 should be $dx_1$. A more minor issue is that I'd avoid the use of subscripts until you're ready to make the final result look pretty.
– JimB
Nov 11 '18 at 1:00

I've made a few small changes in the variable names (used "pdf" instead of "pr") and removed the extra use of the probability density functions in Expectation as the pdf's are automatically used when the associated distribution is defined.

(* Marginal density of x1 *)
pdfx1 = PDF[NormalDistribution[μ1, σ1], x1]

(* Log of pdfx1 *)
logpdfx1 = Log[pdfx1] //.
{Log[a_ b_] -> Log[a] + Log[b], Log[Exp[a_]] -> a,  Log[1/a_] -> -Log[a]}
(* -((x1 - μ1)^2/(2 σ1^2)) + 1/2 (-Log - Log[π]) - Log[σ1] *)

(* H1 two ways *)
H1 = Expectation[-logpdfx1, x1 \[Distributed] NormalDistribution[μ1, σ1],
Assumptions -> σ1 > 0]
(* 1/2 (1+Log+Log[π]+2 Log[σ1]) *)
H1 = Integrate[-logpdfx1 pdfx1, {x1, -∞, ∞}, Assumptions -> σ1 > 0]
(* 1/2 (1+Log[2 π σ1^2]) *)

Now construct H12:

(* Marginal density of x2 *)
pdfx2 = PDF[NormalDistribution[μ2, σ2], x2]
(* Joint density *)
pdfx1x2 = PDF[BinormalDistribution[{μ1, μ2}, {σ1, σ2}, ρ], {x1, x2}]
(* Conditional density *)
pdfx1gx2 = FullSimplify[pdfx1x2/pdfx2]

(* Find log of conditional density *)
logpdfx1gx2 = Log[pdfx1gx2] //.
{Log[a_ b_] -> Log[a] + Log[b], Log[Exp[a_]] -> a, Log[1/a_] -> -Log[a]}
(* ((x2 - μ2) ρ σ1 + (-x1 + μ1) σ2)^2/(2 (-1 + ρ^2) σ1^2 σ2^2) +
1/2 (-Log - Log[π]) + Log[1/Sqrt[1 - ρ^2]] - Log[σ1] *)

(* H12 two ways *)
H12 = Integrate[-logpdfx1gx2 pdfx1gx2, {x1, -∞, ∞},
Assumptions -> {σ1 > 0, σ2 > 0, -1 < ρ < 1}]
(* 1/2 (1+Log[π]+Log[2-2 ρ^2]+2 Log[σ1]) *)
H12 = Expectation[-logpdfx1gx2,
x1 \[Distributed] NormalDistribution[μ1 + ((x2 - μ2) ρ σ1)/σ2, σ1 Sqrt[1 - ρ^2]],
Assumptions -> {σ1 > 0, σ2 > 0, -1 < ρ < 1}]
(* 1/2 (1+Log+Log[π]+2 Log[σ1]) *)

x21 = FullSimplify[H1 - H12] /. Log[a_^k_] -> k Log[a]
(* -(1/2) Log[1-ρ^2] *)

If one is willing to believe that x1 given x2 has a normal distribution, then one just needs to find the mean and variance of that distribution to be able to construct H12.

assumptions = {σ1 > 0, σ2 > 0, -1 < ρ < 1}
(* Direct calculation of mean and variance of conditional distribution *)
meanx1gx2 = Integrate[x1 pdfx1gx2, {x1, -∞, ∞}, Assumptions -> assumptions]
varx1gx2 = FullSimplify[Integrate[x1^2 pdfx1gx2, {x1, -∞, ∞},
Assumptions -> assumptions] - meanx1gx2^2]
(* Alternative method of calculating mean and variance of  conditional distribution *)
distx1gx2 = ProbabilityDistribution[pdfx1gx2, {x1, -∞, ∞}, Assumptions -> assumptions];
meanx1gx2 = Mean[distx1gx2]
(* μ1+((x2-μ2) ρ σ1)/σ2 *)
varx1gx2 = Variance[distx1gx2]
(* -(-1+ρ^2) σ1^2 *)

H12 = Expectation[-logpdfx1gx2, x1 \[Distributed] NormalDistribution[meanx1gx2, varx1gx2^(1/2)],
Assumptions -> assumptions]
(* 1/2 (1+Log+Log[π]+Log[1-ρ^2]+2 Log[σ1]) *)
• It is well-known about the mean, variance, and distribution of $x_1|x_2$ when $x_1$ and $x_2$ are from a bivariate normal distribution. I'll add in something that finds the mean and variance using basic Mathematica commands but I'll leave out showing that $x_1|x_2$ has a normal distribution.
– JimB
Nov 11 '18 at 17:39