matlab几种计算互信息的方法

里面包含几种用于matlab计算两个序列之间的互信息的方法，用于评价变量之间的耦合程度，可以参考和学习。

function [ I1 I2 ] = KraskovMI（ X Y k varargin ）
%KraskovMI computes the Kraskov estimator for the mutual information.
%   1. Input: X Y
%             k: nearest neighbour
%             zeroFix （optional）: fix the negative estimation to 0 （default
%                                 false）;
%
%   univariate: X Y （n x 1） vector
%   multivariate: X Y （n x m） matrix （rows=observations
%   columns=variables）
%
%   2. Output: I1 I2: the two estimator of MI I（1） I（2） （see Ref.）
%
% Ref: Kraskov Alexander Harald St鰃bauer and Peter Grassberger.
%      “Estimating mutual information.“ Physical review E 69.6 （2004）: 066138.
%
% Author: Paolo Inglese 
% Last revision: 17-05-2015

if nargin < 3 || nargin > 4
    error（‘Wrong input number.‘）;
end
if nargin == 3
    zeroFix = false;
end
if nargin == 4
    if ~islogical（varargin{1}）
        error（‘zeroFix must be true or false‘）;
    else
        zeroFix = varargin{1};
    end
end
    

if size（X 1） ~= size（Y 1）
    error（‘X and Y must contain the same number of samples‘）;
end

nObs = size（X 1）;

% compute distance between each sample and its k-th nearest neighbour
dz = zeros（nObs nObs）;
dx = zeros（nObs nObs）;
dy = zeros（nObs nObs）;
for i = 1:nObs
    for j = 1:nObs
        dx（ij） = sqrt（sum（（X（i :） - X（j :））.^2））;
        dy（ij） = sqrt（sum（（Y（i :） - Y（j :））.^2））;
        dz（ij） = max（[dx（i j） dy（i j）]）;
    end
end

% find nx（i） and ny（i）
Eps = zeros（nObs 1）;
Nn = zeros（nObs 1）;

nx1 = zeros（nObs 1）;
ny1 = zeros（nObs 1）;
nx2 = zeros（nObs 1）;
ny2 = zeros（nObs 1）;
for i = 1:nObs
    
    dxSample = dx（i :）;
    dxSample（i） = [];
    
    dySample = dy（i :）;
    dySample（i） = [];
    
    dzSample = dz（i :）;
    dzSample（i） = [];
    [EpsSample NnSample] = sort（dzSample ‘ascend‘）;
    Eps（i） = EpsSample（k）;
    Nn（i） = NnSample（k）;
    
    nx1（i） = sum（dxSample < Eps（i））;
    ny1（i） = sum（dySample < Eps（i））;
    
    nx2（i） = sum（dxSample <= Eps（i））;
    ny2（i） = sum（dySample <= Eps（i））;
    
end

% mutual information estimators
I1 = psi（k） - sum（psi（nx1 + 1） + psi（ny1 + 1）） / nObs + psi（nObs）;
I2 = psi（k） - 1/k - sum（psi（nx2） + psi（ny2）） / nObs + psi（nObs）;

if （zeroFix）
    if I1 < 0
        warning（‘First estimator is negative -> 0‘）;
        I1 = 0;
    end
    if I2 < 0
        warning（‘Second estimator is negative -> 0‘）;
        I2 = 0;
    end
end

end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       5843  2004-06-04 04:01  互信息1\1ICAtests.C

     文件       2766  2004-06-04 10:23  互信息1\ICAtests.m

     文件        502  2005-11-23 14:34  互信息1\makefile

     文件        286  2005-09-26 21:36  互信息1\mamari.m

     文件       7918  2004-06-04 04:00  互信息1\MIClustering.C

     文件        819  2004-06-04 10:23  互信息1\MIClustering.m

     文件       3254  2004-06-04 03:59  互信息1\MIhigherdim.C

     文件        882  2004-06-04 10:23  互信息1\MIhigherdim.m

     文件       9370  2004-06-04 04:17  互信息1\milca.C

     文件       1505  2018-01-20 23:03  互信息1\milca.m

     文件      10459  2005-11-23 14:32  互信息1\milcadelay.C

     文件      84992  2005-11-23 15:25  互信息1\milcadelay.exe

     文件       1872  2004-06-04 10:23  互信息1\milcadelay.m

     文件     320832  2004-06-04 10:23  互信息1\milcaECG8ch.mat

     文件      60784  2004-06-04 04:13  互信息1\miutils.C

     文件       7957  2004-06-04 03:21  互信息1\miutils.h

     文件       3417  2004-06-04 04:00  互信息1\MIxnyn.C

     文件        900  2004-06-04 10:23  互信息1\MIxnyn.m

     文件       6185  2004-06-04 10:23  互信息1\Tutorial.m

     文件       1034  2006-01-06 11:29  互信息2\mi\condentropy.m

     文件       1237  2006-01-06 11:29  互信息2\mi\condmutualinfo.m

     文件        724  2006-01-06 11:29  互信息2\mi\demo_mi.m

     文件        477  2007-08-23 09:49  互信息2\mi\elementmexheader.h

     文件        761  2006-01-06 11:29  互信息2\mi\entropy.m

     文件       1800  2006-01-06 11:29  互信息2\mi\estcondentropy.cpp

     文件      40960  2006-01-06 11:29  互信息2\mi\estcondentropy.dll.old

     文件       8901  2006-01-06 11:29  互信息2\mi\estcondentropy.mexglx

     文件      54608  2007-08-23 09:49  互信息2\mi\estcondentropy.mexmac

     文件      20480  2009-04-18 17:58  互信息2\mi\estcondentropy.mexw32

     文件       1249  2006-01-06 11:29  互信息2\mi\estentropy.cpp

............此处省略97个文件信息