Newer
Older
function varargout=mi(varargin)
%MI Histogram based mutual information.
% I=MI(X) computes the auto mutual information using 10 bins,
% where X can be a multi-column vector.
%
% I=MI(X1,X2,...,Xn) computes the pair-wise mutual information
% between all pairs of X1,X2,...,Xn using 10 bins. The result I
% is a N x N matrix.
%
% I=MI(X1,X2,...,Xn,N), where N is a scalar, determines the number
% of bins N.
%
% I=MI(X1,X2,...,Xn,L) or I=MI(X1,X2,...,Xn,N,L), where L is a
% scalar, computes mutual informations for shifting of the
% components of each pair of X1,X2,...,Xn until a maximal lag L.
%
% [I S]=MI(...) computes the mutual information and the
% standard error (only for one- and two-dimensional data).
%
% MI(...) without any output arguments opens a GUI for interactively
%
% By using the GUI, the mutual information can be stored into the
% workspace. If their standard error is available, they will be
% appended to the mutual information matrix as the last two columns.
%
% MI(...,param) additional parameters according to the GUI are
% available:
% gui - Creates the GUI.
% nogui - Suppresses the GUI.
% silent - Suppresses all output.
%
% MI without any arguments calls a demo (the same as the example
% below).
%
% Remark
% Please note that the mutual information derived with MI slightly
% differs from the results derived with MIGRAM. The reason is that
% MI also considers estimation errors.
%
%
% References:
% Roulston, M. S.:
% Estimating the errors on measured entropy and mutual
% information, Physica D, 125, 1999.
% Copyright (c) 2008-2009
% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
% http://www.pik-potsdam.de
%
% Copyright (c) 2002-2008
% Andre Sitz/ Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
% $Date$
% $Revision$
%
% $Log$
% Revision 3.15 2010/06/29 12:49:57 marwan
% debugging line removed
%
% Revision 3.14 2009/03/24 08:32:57 marwan
% copyright address changed
%
% Revision 3.13 2007/12/20 16:26:06 marwan
% changed gpl splash behaviour
%
% Revision 3.12 2007/07/17 09:04:38 marwan
% bug for single data series resolved
%
% Revision 3.11 2007/06/19 12:26:28 marwan
% speed enhanced
%
% Revision 3.10 2006/03/29 13:07:55 marwan
% problems regarding OPRPs and embedding resolved
%
% Revision 3.9 2006/02/14 11:46:15 marwan
% *** empty log message ***
%
% Revision 3.8 2005/03/16 11:19:02 marwan
% help text modified
%
% Revision 3.7 2004/11/10 07:05:09 marwan
% initial import
%
%
% This program is part of the new generation XXII series.
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or any later version.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties
global props
init_properties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input
error(nargchk(0,99,nargin));
if nargout>2, error('Too many output arguments'), end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% splash the GPL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% read the input
set(0,'ShowHidden','on')
delete(findobj('Tag','msgbox'))
nogui=0;
lag=[];
nbin=10;
x={}; sm_length=[];
if nargin && isnumeric(varargin{1})
sum_m=0;
y=varargin{i};
if size(y,1)==1, y=y'; end
my=size(y,2); sum_m=sum_m+my;
mx=size(y,1);
if isempty(sm_length), sm_length=mx; end
if mx < sm_length, sm_length=mx; end
x(i)={y};
end
end
if isempty(x), error('Not a valid input vector.'), end
i_double=find(cellfun('isclass',varargin,'double'));
i_char=find(cellfun('isclass',varargin,'char'));
if length(i_double)>1
if max(size(varargin{i_double(end-1)}))==1
nbin=varargin{i_double(end-1)};
end
end
if max(size(varargin{i_double(end)}))==1
lag=varargin{i_double(end)}; else if isempty(lag), lag=0; end
end
action='init';
if ~isempty(i_char)
if findstr(lower(varargin{i_char(1)}(1)),'n')
if findstr(lower(varargin{i_char(1)}(1)),'s')
if findstr(lower(varargin{i_char(1)}(1)),'g')
nogui=1;
action='compute';
end
if nogui==-1; nogui=0; end
if isempty(findobj('Tag','MI_Fig')) && isempty(findobj('Tag','MI_Fig'))
% h=errordlg('Start the programme with mi(vector,vector).');
% set(h,'Tag','msgbox')
action='end';
else
action=varargin{1};
h=findobj('Tag','MI_Fig');
x=get(h(1),'UserData');
h=findobj('Tag','nBin');
nbin=10;
lag=40;
nogui=0;
else
if ischar(varargin{1})
action=varargin{1};
end
end
switch(action)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% create gui
case 'init'
h8=figure(props.window,... % Plot Figure
'Tag','MI_Fig',...
'MenuBar','Figure',...
'Position',[69.5000 36.6429 120.0000 30.0714],...
'PaperPosition',[0.25 0.25 7.7677 11.193],...
'PaperOrientation','portrait',...
'UserData',x,...
'Name','Mutual Information');
set(0,'showhidden','on')
h=findobj('Label','&Help','Type','uimenu');
if isempty(h)
h=uimenu('Label','&Help');
uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi');
else
h1=flipud(get(h(1),'Children'));
set(h1(1),'Separator','on')
uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi');
copyobj(h1,h(1))
delete(h1)
end
h=axes(props.axes,...
'Position',[89 24.8 6.8 3.5]);
logo=load('logo');
h2=imagesc([logo.logo fliplr(logo.logo)]);
set(h2,'Tag','uniLogo')
set(h,props.logo,'Tag','axes_logo')
h=uicontrol(props.text,...
'Tag','text_logo',...
h2=textwrap(h,{[char(169),' AGNLD'],'University of Potsdam','2002-2006'});
'Tag','mi_axes',...
'Box','On',...
'Position',[9 3.5 72.8333 23.0714]);
'Tag','frame',...
'Position',[91.5000 1.3571 23.5000 20.7857]);
'Tag','text',...
'String','Number of bins',...
'Position',[93.8333 20 17.8333 1.5000]);
'Tag','nBin',...
'String',num2str(nbin),...
'ToolTip','Number of bins for estimation of the probability function.',...
'Position',[93.8333 18.5714 17.8333 1.5000]);
'Tag','text',...
'String','max. Lag',...
'Position',[93.8333 16.5 17.8333 1.5000]);
'Tag','maxLag',...
'String',num2str(lag),...
'ToolTip','Maximal lag.',...
'Position',[93.8333 15.0714 17.8333 1.5000]);
'Tag','text_show',...
'String','Show components',...
'Enable','off',...
'Position',[93.8333 12.5 21.8333 1.5000]);
'Tag','edit_show1',...
'String','1|2',...
'Enable','off',...
'Callback','mi plot_mi',...
'Position',[93.8333 11.0714 7 1.5000]);
'Tag','text_show',...
'String','vs.',...
'Enable','off',...
'HorizontalAlign','center',...
'Position',[101.24995 11.0714-.2 4 1.5000]);
'Tag','edit_show2',...
'String','1|2',...
'Callback','mi plot_mi',...
'Enable','off',...
'Position',[105.6666 11.0714 7 1.5000]);
'String','Store',...
'Tag','button_store',...
'Enable','Off',...
'ToolTip','Stores the mutual information into a variable in the workspace.',...
'Callback','mi store',...
'Position',[103.8666 8.2143 8.8 2.2143]);
'Tag','button_print',...
'CallBack','mi print',...
'ToolTip','Prints the MI window.',...
'String','Print',...
'Position',[93.8333 8.2143 8.8 2.2143]);
'Tag','button_close',...
'CallBack','mi close',...
'ToolTip','Closes the MI window.',...
'String','Close',...
'Position',[93.8333 2.6429 18.8333 2.2143]);
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
'String','Apply',...
'Tag','button_apply',...
'ToolTip','Starts the computation.',...
'Callback','mi compute',...
'Position',[93.8333 5.4286 18.8333 2.2143]);
set(h8, 'HandleVis','CallBack')
tags={'MI_Fig';'axes_logo';'text_logo';'frame';'text';'mi_axes';'nBin';'maxLag';'text_show';'edit_show1';'edit_show2';'button_store';'button_print';'button_close';'button_apply';};
h=[];
for i=1:length(tags); h=[h; findobj('Tag',tags{i})]; end
set(h,'Units','Norm')
mi('compute')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close windows
case 'close'
set(0,props.root)
h=findobj('Tag','MI_Fig');
if ~isempty(h), close(h(1)), end
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% store
case 'store'
if ~isempty(findobj('Tag','button_store'))
h=findobj('Tag','button_store');
h1=findobj('Tag','button_close');
if ~isempty(h1), vname_old=get(h1(1),'UserData'); else vname_old=''; end
if isempty(vname_old), vname_old=''; end
vname=char(inputdlg('Choose a variable name.','Store output',1,{vname_old}));
if isempty(vname)
return
else
MI=get(h(1),'UserData');
assignin('base',vname, [MI{1:end}])
h=helpdlg(['Mutual information has been assigned to the workspace variable ''',vname,'''.',10,...
'If available, the standard errors are given in the second column.'],'Store output');
set(h,'Tag','msgbox')
set(h1(1),'UserData',vname)
end
end
set(0,props.root)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print
case 'print'
h=findobj('Tag','uniLogo');
h_axes=findobj('Tag','mi_axes','Parent',gcf); h_axes=h_axes(1);
h=[h; findobj('Tag','text','Parent',gcf)];
h=[h; findobj('Tag','text_logo','Parent',gcf)];
h=[h; findobj('Tag','frame','Parent',gcf)];
h=[h; findobj('Tag','nBin','Parent',gcf)];
h=[h; findobj('Tag','maxLag','Parent',gcf)];
h=[h; findobj('Tag','text_show','Parent',gcf)];
h=[h; findobj('Tag','edit_show1','Parent',gcf)];
h=[h; findobj('Tag','edit_show2','Parent',gcf)];
h=[h; findobj('Tag','button_store','Parent',gcf)];
h=[h; findobj('Tag','button_print','Parent',gcf)];
h=[h; findobj('Tag','button_close','Parent',gcf)];
h=[h; findobj('Tag','button_apply','Parent',gcf)];
set(h,'Visible','Off')
set(h_axes,'Units','Character')
old_pos=get(h_axes,'Position');
set(h_axes,'Units','normalize','Position',[0.1300 0.1100 0.7750 0.8150])
h_dlg=printdlg;
waitfor(h_dlg)
set(h_axes,'Units','Character','Position',old_pos)
set(h_axes,'Units','normalize')
set(h,'Visible','On')
set(0,props.root)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% computation
case 'compute'
warning off
if ~nogui
h_fig=findobj('tag','MI_Fig');
setptr(gcf,'watch'),
obj=({'maxLag','nBin','button_store','button_print','button_close','text','text_show','edit_show1','edit_show2'});
for j=1:length(obj);
h=findobj('Tag',obj{j},'Parent',h_fig(1));
if ~isempty(h)
set(h,'Enable','Off')
end
end
h=findobj('tag','button_apply');
set(h(1),'ToolTip','Stops the computation.','String','Stop','Callback','set(0,''ShowHidden'',''on'');h=findobj(''tag'',''button_apply'');set(h(1),''String'',''Stopped'');set(0,''ShowHidden'',''off'')')
end
if nogui~=2; hw=waitbar(0,'Estimation progress'); end
for t=0:lag,
if ~nogui
set(0,'ShowHidden','on')
h=findobj('tag','button_apply','Parent',h_fig(1));
if strcmpi(get(h(1),'string'),'stopped')
MI((t:lag)+1,1)=NaN;
MI_sigma((t:lag)+1,1)=NaN;
break
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the distributions
for k=1:m;
for l=1:m;
P2=hist2(x{k},x{l},nbin,t,'silent'); % 2D distribution
mP=length(size(P2));
%disp(['k: ',num2str(k),' l: ',num2str(l),' test: ',num2str(sum(~P2(:))/numel(P2))])
if sum(~P2(:))/numel(P2) > 0.9995
warning on
warning('Too less data points for the estimation of the joint distribution.');
MI(k,l,(t:lag)+1)=NaN;
MI_sigma(k,l,(t:lag)+1)=NaN;
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
if ~nogui,
h=findobj('tag','button_apply');
set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','mi compute')
for j=1:length(obj);
h=findobj('Tag',obj{j},'Parent',h_fig(1));
if ~isempty(h)
set(h,'Enable','On')
end
end
setptr(h_fig(1),'arrow')
end
warning on
if nogui~=2, delete(hw); end
return
end
P2=P2/sum(P2(:)); % normalization
clear P1x
if mP==2
P1x=sum(P2,2); % distribution for x
P1y=sum(P2); % distribution for y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the mutual information
I1=[-sum((P1x(P1x~=0)).*log(P1x(P1x~=0))), -sum((P1y(P1y~=0)).*log(P1y(P1y~=0)))]; % entropies of Px and Py
I2=-sum(P2(P2~=0).*log(P2(P2~=0))); % entropy of joint distribution Pxy
I2_syserr=(length(P1x(P1x~=0))+length(P1y(P1y~=0))-length(P2(P2~=0).*log(P2(P2~=0)))-1)/(2*length(x{k})); % standard error for estimation of entropy
MI(k,l,t+1)=I1(1)+I1(2)-I2+I2_syserr;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the standard errors
P2xy=P1x*P1y;
i=(P2xy~=0 & P2~=0);
MI_sigma(k,l,t+1)=sqrt( sum( (log(P2xy(i)./P2(i))+MI(t+1)).^2 .*(P2(i).*(1-P2(i)))) /length(x{k}));
else
I1=[];
for i=1:mP;
P21=permute(P2,[i,i+1:mP, 1:(i-1)]);
P1x(:,i)=sum(reshape(P21,size(P21,2),size(P21,2)^(mP-1)),2); % distribution for x_i
I1 = -sum((P1x(P1x~=0)).*log(P1x(P1x~=0))); % sum of entropies of Px_i
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the mutual information
I2=-sum(P2(P2~=0).*log(P2(P2~=0))); % entropy of joint distribution Pxy
I2_syserr=(length(P1x(P1x~=0))-length(P2(P2~=0).*log(P2(P2~=0)))-1)/(2*length(x{k})); % standard error for estimation of entropy
MI(k,l,t+1)=I1-I2+I2_syserr;
if ~nogui && t*10/lag == ceil(t*10/lag)
h=findobj('Tag','button_store');
if m==2
out_str{1}=MI; out_str{2}=MI_sigma;
else
out_str{1}=MI;
end
set(h(1),'UserData',out_str)
mi('plot_mi')
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
h=findobj('tag','button_apply');
set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','mi compute')
for j=1:length(obj);
h=findobj('Tag',obj{j},'Parent',h_fig(1));
if ~isempty(h)
set(h,'Enable','On')
end
end
end
if ~isempty(findobj('Tag','MI_Fig'))
h=findobj('Tag','button_store');
if m==2
out_str{1}=MI; out_str{2}=MI_sigma;
else
out_str{1}=MI;
end
set(h(1),'Enable','On',...
'UserData',out_str)
end
if nargout==1
varargout(1)={MI};
elseif nargout==2
varargout(1)={MI};
varargout(2)={MI_sigma};
end
if ~nogui, setptr(h_fig(1),'arrow'), end
warning on
if isnan(MI(1,1)), warning('Mutual information is NaN. Use smaller bin size.'), end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output
case 'plot_mi'
h_fig = findobj('Tag','MI_Fig');
if ~isempty(h_fig)
h=findobj('Tag','button_store','Parent',h_fig(1));
if isempty(h), return, end
x=get(h(1),'UserData');
MI=x{1}; if m==2; MI_sigma=x{2}; end
set(0,'current',h_fig(1))
h_axes=findobj('Tag','mi_axes','Parent',h_fig(1)); h_axes=h_axes(1);
h_res=findobj('Tag','mi_result','Parent',h_fig(1)); delete(h_res);
h_show=[findobj('Tag','text_show','Parent',h_fig(1));...
findobj('Tag','edit_show1','Parent',h_fig(1));...
findobj('Tag','edit_show2','Parent',h_fig(1))];
if lag==0
set(h_axes,'Visible','off'),cla
set(h_show,'enable','off')
if m==2;
for k=1:length(MI)
out_str(k)={sprintf(repmat(['%8.3f ',char(177),'%6.3f '],1,2*length(MI)),reshape([MI(k,:);MI_sigma(k,:)],1,2*length(MI)))};
end
else
for k=1:length(MI)
out_str(k)={sprintf(repmat('%8.3f ',1,length(MI)),MI(k,:))};
end
end
h=uicontrol(props.text,...
'Tag','mi_result',...
'HorizontalAlignment', 'left',...
'FontWeight', 'bold',...
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
ex=get(h,'Extent'); set(h,'Position',[11 21 ex(3:4)]);
h=uicontrol(props.listbox,...
'Tag','mi_result',...
'Position',[9 7 72.8333 3.0714],...
'HorizontalAlignment', 'left',...
'String',out_str);
ex1=get(h,'Extent'); if ex1(3)>75, ex1(3)=75; end
if m<=10
set(h,'Position',[11 21-m*ex(4) ex1(3) ex1(4)*m]);
else
set(h,'Position',[11 21-10*ex(4) ex1(3) ex1(4)*10]);
end
else
set(h_axes,'Visible','on'); tx=sprintf('%i|',1:m);
set(h_show(end-1:end),'String',tx(1:end-1))
plot(0)
dimx=get(h_show(end-1),'Value');
dimy=get(h_show(end),'Value');
if m==2;
s1=[permute(MI(dimx,dimy,:),[3,1,2])+3*permute(MI_sigma(dimx,dimy,:),[3,1,2]); flipud(permute(MI(dimx,dimy,:),[3,1,2])-3*permute(MI_sigma(dimx,dimy,:),[3,1,2]))]';
s3=[0:size(MI,3)-1, size(MI,3)-1:-1:0];
if all(size(s3) == size(s1))
patch(s3,s1,-10,'FaceColor',[.9 .9 1],'EdgeColor',[.85 .85 1])
end
hold on
plot(0:size(MI,3)-1,permute(MI(dimx,dimy,:),[3,1,2])), grid on
xlabel('Lag'), ylabel('Mutual Information')
hold off
set(gca,'Tag','mi_axes','layer','top','xlim',[0 lag])
end
drawnow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the end
case 'end'
end
set(0,'ShowHidden','off')
if ~isempty(findobj('Tag','TMWWaitbar')), delete(findobj('Tag','TMWWaitbar')), end
cmd={'mean';'var';'std';'median';'squmean';'geomean';'bias';'skewness';'kurtosis'};
z_whos=whos;x_lasterr=lasterr;y_lastwarn=lastwarn;if nargin;in=varargin{1};else in.class='no input given';end
in=whos('in');
if ~strcmpi(lasterr,'Interrupt')
fid=fopen('error.log','w');
fprintf(fid,'%s\n','Please send us the following error report. Provide a brief');
fprintf(fid,'%s\n','description of what you were doing when this problem occurred.');
fprintf(fid,'%s\n','E-mail or FAX this information to us at:');
fprintf(fid,'%s\n',' E-mail: marwan@pik-potsdam.de');
fprintf(fid,'%s\n',' Fax: ++49 +331 288 2640');
fprintf(fid,'%s\n\n\n','Thank you for your assistance.');
fprintf(fid,'%s\n',repmat('-',50,1));
fprintf(fid,'%s\n',datestr(now,0));
fprintf(fid,'%s\n',['Matlab ',char(version),' on ',computer]);
fprintf(fid,'%s\n',repmat('-',50,1));
fprintf(fid,'%s\n',x_lasterr);
fprintf(fid,'%s\n',y_lastwarn);
fprintf(fid,'%s\n',[' during ==> mi:',action]);
fprintf(fid,'%s',[' input ==> ',in.class]);
if ~isempty(in2), fprintf(fid,'\t%s\n',[' (',in2,')']); end
fprintf(fid,'%s\n',' errorcode ==> no errorcode available');
fprintf(fid,'%s\n',' workspace dump ==>');
fprintf(fid,'%s\n',['Name',char(9),'Size',char(9),'Bytes',char(9),'Class']);
for j=1:length(z_whos);
fprintf(fid,'%s',[z_whos(j).name,char(9),num2str(z_whos(j).size),char(9),num2str(z_whos(j).bytes),char(9),z_whos(j).class]);
if ~strcmp(z_whos(j).class,'cell') && ~strcmp(z_whos(j).class,'struct')
content=eval(z_whos(j).name);
content=mat2str(content(1:min([size(content,1),500]),1:min([size(content,2),500])));
fprintf(fid,'\t%s',content(1:min([length(content),500])));
elseif strcmp(z_whos(j).class,'cell')
content=eval(z_whos(j).name);
fprintf(fid,'\t');
for j2=1:min([length(content),500])
fprintf(fid,'{%s} ',content{j2});
end
elseif strcmp(z_whos(j).class,'struct')
content=fieldnames(eval(z_whos(j).name));
content=char(content); content(:,end+1)=' '; content=content';
fprintf(fid,'\t%s',content(:)');
disp('----------------------------');
disp(' ERROR OCCURED');
disp(' during executing mi');
disp('----------------------------');
disp(x_lasterr);
disp([' during ',action]);
disp('----------------------------');
disp(' Please send us the error report. For your convenience, ')
disp(' this information has been recorded in: ')
disp([' ',fullfile(pwd,'error.log')]), disp(' ')
disp(' Provide a brief description of what you were doing when ')
disp(' this problem occurred.'), disp(' ')
disp(' E-mail or FAX this information to us at:')
disp(' E-mail: marwan@pik-potsdam.de')
disp(' Fax: ++49 +331 288 2640'), disp(' ')
disp(' Thank you for your assistance.')
warning('on')
end
try
if ~nogui
setptr(h_fig(1),'arrow')
h=findobj('tag','button_apply');
set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','mi compute')
for j=1:length(obj);
h=findobj('Tag',obj{j},'Parent',h_fig(1));
if ~isempty(h)
set(h,'Enable','On')
end
end
end
catch
end