From a2465c6ea46be48409c7a059ede1334bf8d21ac2 Mon Sep 17 00:00:00 2001 From: Robin Ince Date: Thu, 28 Sep 2017 11:45:39 +0100 Subject: [PATCH 1/4] James et al Idep 2 predictor PID for Gaussians --- PID_dep_mvn.m | 91 +++++++++++++++++++++++++++++++++++++++++++++++++++ gauss_mi.m | 22 +++++++++++++ 2 files changed, 113 insertions(+) create mode 100644 PID_dep_mvn.m create mode 100644 gauss_mi.m diff --git a/PID_dep_mvn.m b/PID_dep_mvn.m new file mode 100644 index 0000000..8087ab5 --- /dev/null +++ b/PID_dep_mvn.m @@ -0,0 +1,91 @@ +function pid = PID_dep_mvn(C, varsizes); +% 2 predictor PID with Idep for Gaussians + +vs = varsizes; +if length(vs)~=3 + error('only bivaraite PID supported') +end +if sum(vs)~=size(C,1) || size(C,1)~=size(C,2) + error('problem with variable specifications') +end + +mivs = [vs(1)+vs(2) vs(3)]; + +% variable indexes +xidx = 1:vs(1); +yidx = (vs(1)+1):(vs(1)+vs(2)); +zidx = (vs(1)+vs(2)+1):(vs(1)+vs(2)+vs(3)); + +% extract blockwise covariance components +Cxx = C(xidx,xidx); +Cyy = C(yidx,yidx); +Czz = C(zidx,zidx); +Cxy = C(xidx,yidx); +Cxz = C(xidx,zidx); +Cyz = C(yidx,zidx); + +Cdiag = zeros(size(C)); +Cdiag(xidx,xidx) = Cxx; +Cdiag(yidx,yidx) = Cyy; +Cdiag(zidx,zidx) = Czz; + +IXY = gauss_mi(C,mivs); +IY = gauss_mi(C([yidx zidx]:[yidx zidx]), vs(2:3)); +Ct = zeros(vs(1)+vs(3)); +Ct(xidx,xidx) = Cxx; +Ct(zidx,zidx) = Czz; +Ct(xidx,zidx) = Cxz; +Ct(zidx,xidx) = Cxz'; +IX = gauss_mi(Ct, [vs(1) vs(3)]); + + +% unique information in X +Xdelta = zeros(1,4); +Ctest_without = Cdiag; +Ctest_with = Ctest_without; +Ctest_with(xidx,zidx) = Cxz; +Ctest_with(zidx,xidx) = Cxz'; +Xdelta(1) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); + +Ctest_without = Cdiag; +Ctest_without(xidx,yidx) = Cxy; +Ctest_without(yidx,xidx) = Cxy'; +Ctest_with = Ctest_without; +Ctest_with(xidx,zidx) = Cxz; +Ctest_with(zidx,xidx) = Cxz'; +Xdelta(2) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); + +Ctest_without = Cdiag; +Ctest_without(yidx,zidx) = Cyz; +Ctest_without(zidx,yidx) = Cyz'; +Ctest_with = Ctest_without; +Ctest_with(xidx,zidx) = Cxz; +Ctest_with(zidx,xidx) = Cxz'; +Xdelta(3) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); + +Ctest_without = Cdiag; +Ctest_without(xidx,yidx) = Cxy; +Ctest_without(yidx,xidx) = Cxy'; +Ctest_without(yidx,zidx) = Cyz; +Ctest_without(zidx,yidx) = Cyz'; +Ctest_with = Ctest_without; +Ctest_with(xidx,zidx) = Cxz; +Ctest_with(zidx,xidx) = Cxz'; +Xdelta(4) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); + +Xunq = min(Xdelta); + +% fill out PID +Id = zeros(1,4); +Id(2) = Xunq; +Id(1) = IX - Xunq; +Id(3) = IY - Id(1); +Id(4) = IXY - sum(Id(1:3)); + + + + + + + + diff --git a/gauss_mi.m b/gauss_mi.m new file mode 100644 index 0000000..2bc22e0 --- /dev/null +++ b/gauss_mi.m @@ -0,0 +1,22 @@ +function I = gauss_mi(C, varsizes) + +vs = varsizes; +if length(vs)>2 + error('bivariate MI only') +end +if sum(vs) ~= size(C,1) + error('variables incorrectly specified') +end + +CX = C(1:vs(1),1:vs(1)); +CY = C((vs(1)+1):end, (vs(1)+1):end); + +% use closed form expression +chX = chol(CX); +chY = chol(CY); +chXY = chol(C); +% normalisations cancel for information +HX = sum(log(diag(chX))); % + 0.5*Nvarx*log(2*pi*exp(1)); +HY = sum(log(diag(chY))); % + 0.5*Nvary*log(2*pi*exp(1)); +HXY = sum(log(diag(chXY))); % + 0.5*(Nvarx+Nvary)*log(2*pi*exp(1)); +I = (HX + HY - HXY) / log(2); From 4024f2d6d59dc67dd37fcbc4316f109eb21d427a Mon Sep 17 00:00:00 2001 From: Robin Ince Date: Thu, 28 Sep 2017 12:04:19 +0100 Subject: [PATCH 2/4] clean up --- PID_dep_mvn.m | 20 ++++++++++++-------- gauss_mi.m | 1 + 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/PID_dep_mvn.m b/PID_dep_mvn.m index 8087ab5..d1d7c03 100644 --- a/PID_dep_mvn.m +++ b/PID_dep_mvn.m @@ -1,4 +1,4 @@ -function pid = PID_dep_mvn(C, varsizes); +function Id = PID_dep_mvn(C, varsizes); % 2 predictor PID with Idep for Gaussians vs = varsizes; @@ -17,9 +17,12 @@ zidx = (vs(1)+vs(2)+1):(vs(1)+vs(2)+vs(3)); % extract blockwise covariance components -Cxx = C(xidx,xidx); -Cyy = C(yidx,yidx); -Czz = C(zidx,zidx); +% Cxx = C(xidx,xidx); +% Cyy = C(yidx,yidx); +% Czz = C(zidx,zidx); +Cxx = eye(vs(1)); +Cyy = eye(vs(2)); +Czz = eye(vs(3)); Cxy = C(xidx,yidx); Cxz = C(xidx,zidx); Cyz = C(yidx,zidx); @@ -30,12 +33,13 @@ Cdiag(zidx,zidx) = Czz; IXY = gauss_mi(C,mivs); -IY = gauss_mi(C([yidx zidx]:[yidx zidx]), vs(2:3)); +IY = gauss_mi(C([yidx zidx],[yidx zidx]), vs(2:3)); Ct = zeros(vs(1)+vs(3)); +tzidx = (vs(1)+1):(vs(1)+vs(3)); Ct(xidx,xidx) = Cxx; -Ct(zidx,zidx) = Czz; -Ct(xidx,zidx) = Cxz; -Ct(zidx,xidx) = Cxz'; +Ct(tzidx,tzidx) = Czz; +Ct(xidx,tzidx) = Cxz; +Ct(tzidx,xidx) = Cxz'; IX = gauss_mi(Ct, [vs(1) vs(3)]); diff --git a/gauss_mi.m b/gauss_mi.m index 2bc22e0..54c9580 100644 --- a/gauss_mi.m +++ b/gauss_mi.m @@ -1,5 +1,6 @@ function I = gauss_mi(C, varsizes) +% C = C + 0.1*eye(size(C)); vs = varsizes; if length(vs)>2 error('bivariate MI only') From 2bc91b00127cf8ede316611e56f2f778b39ff23d Mon Sep 17 00:00:00 2001 From: Robin Ince Date: Thu, 1 Mar 2018 16:37:10 +0000 Subject: [PATCH 3/4] Idep for Gaussians based on Kay & Ince closed form solutions --- PID_dep_mvn.m | 88 ++++++++++++++++++++++----------------------------- 1 file changed, 37 insertions(+), 51 deletions(-) diff --git a/PID_dep_mvn.m b/PID_dep_mvn.m index d1d7c03..e1c233e 100644 --- a/PID_dep_mvn.m +++ b/PID_dep_mvn.m @@ -3,13 +3,14 @@ vs = varsizes; if length(vs)~=3 - error('only bivaraite PID supported') + error('only bivariate PID supported') end if sum(vs)~=size(C,1) || size(C,1)~=size(C,2) error('problem with variable specifications') end mivs = [vs(1)+vs(2) vs(3)]; +IXY = gauss_mi(C,mivs); % variable indexes xidx = 1:vs(1); @@ -17,22 +18,22 @@ zidx = (vs(1)+vs(2)+1):(vs(1)+vs(2)+vs(3)); % extract blockwise covariance components -% Cxx = C(xidx,xidx); -% Cyy = C(yidx,yidx); -% Czz = C(zidx,zidx); -Cxx = eye(vs(1)); -Cyy = eye(vs(2)); -Czz = eye(vs(3)); +Cxx = C(xidx,xidx); +Cyy = C(yidx,yidx); +Czz = C(zidx,zidx); Cxy = C(xidx,yidx); Cxz = C(xidx,zidx); Cyz = C(yidx,zidx); -Cdiag = zeros(size(C)); -Cdiag(xidx,xidx) = Cxx; -Cdiag(yidx,yidx) = Cyy; -Cdiag(zidx,zidx) = Czz; +chCxx = chol(Cxx); +chCyy = chol(Cyy); +chCzz = chol(Czz); + +% Kay & Ince eq. D2 +P = pinv(chCxx)*Cxy*pinv(chCyy); +Q = pinv(chCxx)*Cxz*pinv(chCzz); +R = pinv(chCyy)*Cyz*pinv(chCzz); -IXY = gauss_mi(C,mivs); IY = gauss_mi(C([yidx zidx],[yidx zidx]), vs(2:3)); Ct = zeros(vs(1)+vs(3)); tzidx = (vs(1)+1):(vs(1)+vs(3)); @@ -42,44 +43,32 @@ Ct(tzidx,xidx) = Cxz'; IX = gauss_mi(Ct, [vs(1) vs(3)]); - -% unique information in X -Xdelta = zeros(1,4); -Ctest_without = Cdiag; -Ctest_with = Ctest_without; -Ctest_with(xidx,zidx) = Cxz; -Ctest_with(zidx,xidx) = Cxz'; -Xdelta(1) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); - -Ctest_without = Cdiag; -Ctest_without(xidx,yidx) = Cxy; -Ctest_without(yidx,xidx) = Cxy'; -Ctest_with = Ctest_without; -Ctest_with(xidx,zidx) = Cxz; -Ctest_with(zidx,xidx) = Cxz'; -Xdelta(2) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); - -Ctest_without = Cdiag; -Ctest_without(yidx,zidx) = Cyz; -Ctest_without(zidx,yidx) = Cyz'; -Ctest_with = Ctest_without; -Ctest_with(xidx,zidx) = Cxz; -Ctest_with(zidx,xidx) = Cxz'; -Xdelta(3) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); - -Ctest_without = Cdiag; -Ctest_without(xidx,yidx) = Cxy; -Ctest_without(yidx,xidx) = Cxy'; -Ctest_without(yidx,zidx) = Cyz; -Ctest_without(zidx,yidx) = Cyz'; -Ctest_with = Ctest_without; -Ctest_with(xidx,zidx) = Cxz; -Ctest_with(zidx,xidx) = Cxz'; -Xdelta(4) = gauss_mi(Ctest_with,mivs) - gauss_mi(Ctest_without,mivs); - -Xunq = min(Xdelta); +% Dependency lattice edges (Kay & Ince Table 9) +% I(Y; Z) +b = IX; + +halflog2det = @(X) sum(log2(diag(chol(X)))); + +inum = halflog2det(eye(vs(1)) - R*Q'*Q*R'); +iden = halflog2det(eye(vs(2))-Q'*Q) + halflog2det(eye(vs(2))-R'*R); +i = inum - iden - IY; + +knum = halflog2det(eye(vs(1)) - P'*P); +Ct = zeros(sum(vs)); +Ct(xidx,xidx) = eye(vs(1)); +Ct(yidx,yidx) = eye(vs(2)); +Ct(zidx,zidx) = eye(vs(3)); +Ct(xidx,yidx) = P; +Ct(yidx,xidx) = P'; +Ct(xidx,zidx) = Q; +Ct(zidx,xidx) = Q'; +Ct(yidx,zidx) = R; +Ct(zidx,yidx) = R'; +kden = halflog2det(Ct); +k = knum - kden - IY; % fill out PID +Xunq = min([b i k]); Id = zeros(1,4); Id(2) = Xunq; Id(1) = IX - Xunq; @@ -90,6 +79,3 @@ - - - From 66d0e556cb5c91a505d2c6085230063676a10b91 Mon Sep 17 00:00:00 2001 From: Robin Ince Date: Thu, 1 Mar 2018 16:43:35 +0000 Subject: [PATCH 4/4] rename and add docstring --- PID_dep_mvn.m => calc_pi_Idep_mvn.m | 30 +++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) rename PID_dep_mvn.m => calc_pi_Idep_mvn.m (68%) diff --git a/PID_dep_mvn.m b/calc_pi_Idep_mvn.m similarity index 68% rename from PID_dep_mvn.m rename to calc_pi_Idep_mvn.m index e1c233e..7dc12cb 100644 --- a/PID_dep_mvn.m +++ b/calc_pi_Idep_mvn.m @@ -1,17 +1,25 @@ -function Id = PID_dep_mvn(C, varsizes); -% 2 predictor PID with Idep for Gaussians +function Id = calc_pi_Idep_mvn(C, varsizes); +% calc_pi_Idep_mvn(C, varsizes) +% Calculate bivariate PID using Idep approach (James et al. 2017) for +% Gaussians (Kay & Ince 2018). +% C is the full joint covariance matrix with variables in order +% X0 (first predictor), X1 (second predictor), S (target) +% varsizes is a length 3 vector containing the dimensionality of each of +% the above. +% Returns a length 4 vector containing +% [Red Unq_1 Unq_2 Syn] +% +% James et al. 2017 http://arxiv.org/abs/1709.06653 +% Kay & Ince 2018 vs = varsizes; if length(vs)~=3 - error('only bivariate PID supported') + error('calc_pi_Idep_mvn: only bivariate PID supported') end if sum(vs)~=size(C,1) || size(C,1)~=size(C,2) - error('problem with variable specifications') + error('calc_pi_Idep_mvn: problem with variable specifications') end -mivs = [vs(1)+vs(2) vs(3)]; -IXY = gauss_mi(C,mivs); - % variable indexes xidx = 1:vs(1); yidx = (vs(1)+1):(vs(1)+vs(2)); @@ -25,6 +33,7 @@ Cxz = C(xidx,zidx); Cyz = C(yidx,zidx); +% cholesky square root chCxx = chol(Cxx); chCyy = chol(Cyy); chCzz = chol(Czz); @@ -34,6 +43,9 @@ Q = pinv(chCxx)*Cxz*pinv(chCzz); R = pinv(chCyy)*Cyz*pinv(chCzz); +% standard mutual informations +mivs = [vs(1)+vs(2) vs(3)]; +IXY = gauss_mi(C,mivs); IY = gauss_mi(C([yidx zidx],[yidx zidx]), vs(2:3)); Ct = zeros(vs(1)+vs(3)); tzidx = (vs(1)+1):(vs(1)+vs(3)); @@ -43,12 +55,10 @@ Ct(tzidx,xidx) = Cxz'; IX = gauss_mi(Ct, [vs(1) vs(3)]); +halflog2det = @(X) sum(log2(diag(chol(X)))); % Dependency lattice edges (Kay & Ince Table 9) -% I(Y; Z) b = IX; -halflog2det = @(X) sum(log2(diag(chol(X)))); - inum = halflog2det(eye(vs(1)) - R*Q'*Q*R'); iden = halflog2det(eye(vs(2))-Q'*Q) + halflog2det(eye(vs(2))-R'*R); i = inum - iden - IY;