Here are few equations to calculate the quantity of information flow (i..e., mutual information):

enter image description here

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...

  • 1
    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 at 1:00
up vote 5 down vote accepted

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[2] - Log[π]) - Log[σ1] *)

(* H1 two ways *)
H1 = Expectation[-logpdfx1, x1 \[Distributed] NormalDistribution[μ1, σ1],
  Assumptions -> σ1 > 0]
(* 1/2 (1+Log[2]+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[2] - 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[2]+Log[π]+2 Log[σ1]) *)

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

Addendum:

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[2]+Log[π]+Log[1-ρ^2]+2 Log[σ1]) *)
  • 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 at 17:39

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.