diff --git a/mi.m b/mi.m index f02c7be1338483f8f1162e25c7ef99fc76688155..5711efa621ef2f7e26a7b8a3a06caf71721fa7cd 100644 --- a/mi.m +++ b/mi.m @@ -36,7 +36,7 @@ function varargout=mi(varargin) % Examples: x = sin(0:.2:8*pi)' + .1*randn(126,1); % mi(x,40,100) % -% See also HIST2, HISTN, ENTROPY. +% See also HIST2, HISTN, ENTROPY, MIGRAM. % % References: % Roulston, M. S.: @@ -51,6 +51,9 @@ function varargout=mi(varargin) % $Revision$ % % $Log$ +% 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 *** % @@ -73,7 +76,6 @@ function varargout=mi(varargin) global props init_properties -action=''; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input @@ -85,7 +87,7 @@ if nargout>2, error('Too many output arguments'), end splash_gpl('mi'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% error control -try +%try %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% read the input @@ -95,11 +97,13 @@ delete(findobj('Tag','msgbox')) nogui=0; lag=[]; nbin=10; -x={}; sum_m=0; sm_length=[]; -if nargin & isnumeric(varargin{1}) +x={}; sm_length=[]; +if nargin && isnumeric(varargin{1}) + + sum_m=0; for i=1:nargin - if isnumeric(varargin{i}) & max(size(varargin{i}))>1 + if isnumeric(varargin{i}) && max(size(varargin{i}))>1 y=varargin{i}; if size(y,1)==1, y=y'; end my=size(y,2); sum_m=sum_m+my; @@ -109,7 +113,7 @@ if nargin & isnumeric(varargin{1}) x(i)={y}; end end - m=size(x,2); if sum_m==1; lag=0; x(2)=x(1); sum_m=2; m=2; end + m=size(x,2); if sum_m==1; lag=0; x(2)=x(1); m=2; end if isempty(x), error('Not a valid input vector.'), end i_double=find(cellfun('isclass',varargin,'double')); @@ -139,15 +143,15 @@ if nargin & isnumeric(varargin{1}) nogui=-1; end end - if nargout & (nogui<1) + if nargout && (nogui<1) nogui=1; action='compute'; end if nogui==-1; nogui=0; end -elseif nargin & ischar(varargin{1}) & ~isempty(varargin{1}) +elseif nargin && ischar(varargin{1}) && ~isempty(varargin{1}) - if isempty(findobj('Tag','MI_Fig')) & isempty(findobj('Tag','MI_Fig')) + 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'; @@ -156,17 +160,17 @@ elseif nargin & ischar(varargin{1}) & ~isempty(varargin{1}) h=findobj('Tag','MI_Fig'); x=get(h(1),'UserData'); h=findobj('Tag','nBin'); - nbin=str2num(get(h(1),'String')); + nbin=str2double(get(h(1),'String')); m=size(x,2); h=findobj('Tag','maxLag'); - lag=str2num(get(h(1),'String')); + lag=str2double(get(h(1),'String')); end else action='init'; end if ~nargin - x={sin(0:.2:8*pi)'+.1*randn(126,1)}; x(2)=x; m=2; sum_m=2; + x={sin(0:.2:8*pi)'+.1*randn(126,1)}; x(2)=x; m=2; nbin=10; lag=40; nogui=0; @@ -197,11 +201,11 @@ case 'init' h=findobj('Label','&Help','Type','uimenu'); if isempty(h) h=uimenu('Label','&Help'); - h2=uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi'); + uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi'); else h1=flipud(get(h(1),'Children')); set(h1(1),'Separator','on') - h2=uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi'); + uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi'); copyobj(h1,h(1)) delete(h1) end @@ -219,59 +223,59 @@ case 'init' h2=textwrap(h,{[char(169),' AGNLD'],'University of Potsdam','2002-2006'}); set(h,'String',h2) - h=axes(props.axes,... + axes(props.axes,... 'Tag','mi_axes',... 'Box','On',... 'Position',[9 3.5 72.8333 23.0714]); - h=uicontrol(props.frame,... + uicontrol(props.frame,... 'Tag','frame',... 'Position',[91.5000 1.3571 23.5000 20.7857]); - h=uicontrol(props.text,... + uicontrol(props.text,... 'Tag','text',... 'String','Number of bins',... 'Position',[93.8333 20 17.8333 1.5000]); - h=uicontrol(props.edit,... + uicontrol(props.edit,... 'Tag','nBin',... 'String',num2str(nbin),... 'ToolTip','Number of bins for estimation of the probability function.',... 'Position',[93.8333 18.5714 17.8333 1.5000]); - h=uicontrol(props.text,... + uicontrol(props.text,... 'Tag','text',... 'String','max. Lag',... 'Position',[93.8333 16.5 17.8333 1.5000]); - h=uicontrol(props.edit,... + uicontrol(props.edit,... 'Tag','maxLag',... 'String',num2str(lag),... 'ToolTip','Maximal lag.',... 'Position',[93.8333 15.0714 17.8333 1.5000]); - h=uicontrol(props.text,... + uicontrol(props.text,... 'Tag','text_show',... 'String','Show components',... 'Enable','off',... 'Position',[93.8333 12.5 21.8333 1.5000]); - h=uicontrol(props.popup,... + uicontrol(props.popup,... 'Tag','edit_show1',... 'String','1|2',... 'Enable','off',... 'Callback','mi plot_mi',... 'Position',[93.8333 11.0714 7 1.5000]); - h=uicontrol(props.text,... + uicontrol(props.text,... 'Tag','text_show',... 'String','vs.',... 'Enable','off',... 'HorizontalAlign','center',... 'Position',[101.24995 11.0714-.2 4 1.5000]); - h=uicontrol(props.popup,... + uicontrol(props.popup,... 'Tag','edit_show2',... 'String','1|2',... 'Callback','mi plot_mi',... @@ -279,7 +283,7 @@ case 'init' 'Position',[105.6666 11.0714 7 1.5000]); - h=uicontrol(props.button,... + uicontrol(props.button,... 'String','Store',... 'Tag','button_store',... 'Enable','Off',... @@ -287,21 +291,21 @@ case 'init' 'Callback','mi store',... 'Position',[103.8666 8.2143 8.8 2.2143]); - h=uicontrol(props.button,... + uicontrol(props.button,... 'Tag','button_print',... 'CallBack','mi print',... 'ToolTip','Prints the MI window.',... 'String','Print',... 'Position',[93.8333 8.2143 8.8 2.2143]); - h=uicontrol(props.button,... + uicontrol(props.button,... 'Tag','button_close',... 'CallBack','mi close',... 'ToolTip','Closes the MI window.',... 'String','Close',... 'Position',[93.8333 2.6429 18.8333 2.2143]); - h=uicontrol(props.button,... + uicontrol(props.button,... 'String','Apply',... 'Tag','button_apply',... 'ToolTip','Starts the computation.',... @@ -380,6 +384,9 @@ case 'store' case 'compute' warning off + MI = zeros(m,m,lag); + MI_sigma = zeros(m,m,lag); + if ~nogui h_fig=findobj('tag','MI_Fig'); setptr(gcf,'watch'), @@ -396,82 +403,83 @@ case 'compute' 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 + 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=histn(x{k},x{l},nbin,t,'silent'); % 2D distribution - mP=length(size(P2)); - if sum(~P2(:))/prod(size(P2)) > 0.99 - 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; - 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') + for k=1:m; + for l=1:m; + P2=hist2(x{k},x{l},nbin,t,'silent'); % 2D distribution + mP=length(size(P2)); + if sum(~P2(:))/numel(P2) > 0.99 + 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; + 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; 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')'; % 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))')'; % 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; - end - end, end + end + if nogui~=2; waitbar(t/lag); end - 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') + 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') end - - end if nogui~=2; delete(hw); end if ~nogui @@ -506,13 +514,12 @@ case 'compute' 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 - try, set(0,props.root), end + try set(0,props.root), catch end if nargout==0 mi('plot_mi') end - if nargout==0 & nogui - assignin('base','ans', MI) - ans=MI + if nargout==0 && nogui + varargout{1} = MI; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output @@ -536,6 +543,7 @@ case 'plot_mi' if lag==0 set(h_axes,'Visible','off'),cla set(h_show,'enable','off') + out_str{MI} = []; if m==2; for k=1:length(MI) out_str(k)={sprintf(repmat('%8.3f ą%6.3f ',1,2*length(MI)),reshape([MI(k,:);MI_sigma(k,:)],1,2*length(MI)))}; @@ -549,7 +557,7 @@ case 'plot_mi' 'Tag','mi_result',... 'HorizontalAlignment', 'left',... 'FontWeight', 'bold',... - 'String',['Mutual Information: ']); + 'String','Mutual Information: '); ex=get(h,'Extent'); set(h,'Position',[11 21 ex(3:4)]); h=uicontrol(props.listbox,... 'Tag','mi_result',... @@ -576,14 +584,14 @@ case 'plot_mi' patch(s3,s1,-10,'FaceColor',[.9 .9 1],'EdgeColor',[.85 .85 1]) end hold on - end + end 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 - try, set(0,props.root), end + try set(0,props.root), catch end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the end @@ -596,55 +604,55 @@ set(0,'ShowHidden','off') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% error handling -%if 0 -catch +if 0 +%catch if ~isempty(findobj('Tag','TMWWaitbar')), delete(findobj('Tag','TMWWaitbar')), end - cmd={['mean'];['var'];['std'];['median'];['squmean'];['geomean'];['bias'];['skewness'];['kurtosis']}; + 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 - if ischar(in), in2=in; else, in2=[]; end + if ischar(in), in2=in; else in2=[]; end in=whos('in'); if ~strcmpi(lasterr,'Interrupt') fid=fopen('error.log','w'); - err=fprintf(fid,'%s\n','Please send us the following error report. Provide a brief'); - err=fprintf(fid,'%s\n','description of what you were doing when this problem occurred.'); - err=fprintf(fid,'%s\n','E-mail or FAX this information to us at:'); - err=fprintf(fid,'%s\n',' E-mail: marwan@agnld.uni-potsdam.de'); - err=fprintf(fid,'%s\n',' Fax: ++49 +331 977 1142'); - err=fprintf(fid,'%s\n\n\n','Thank you for your assistance.'); - err=fprintf(fid,'%s\n',repmat('-',50,1)); - err=fprintf(fid,'%s\n',datestr(now,0)); - err=fprintf(fid,'%s\n',['Matlab ',char(version),' on ',computer]); - err=fprintf(fid,'%s\n',repmat('-',50,1)); - err=fprintf(fid,'%s\n',x_lasterr); - err=fprintf(fid,'%s\n',y_lastwarn); - err=fprintf(fid,'%s\n',[' during ==> mi:',action]); - err=fprintf(fid,'%s',[' input ==> ',in.class]); - if ~isempty(in2), err=fprintf(fid,'\t%s\n',[' (',in2,')']); end - err=fprintf(fid,'%s\n',[' errorcode ==> no errorcode available']); - err=fprintf(fid,'%s\n',' workspace dump ==>'); + 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@agnld.uni-potsdam.de'); + fprintf(fid,'%s\n',' Fax: ++49 +331 977 1142'); + 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 ==>'); if ~isempty(z_whos), - err=fprintf(fid,'%s\n',['Name',char(9),'Size',char(9),'Bytes',char(9),'Class']); - for j=1:length(z_whos); - err=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]))); - err=fprintf(fid,'\t%s',content(1:min([length(content),500]))); - elseif strcmp(z_whos(j).class,'cell') - content=eval(z_whos(j).name); - err=fprintf(fid,'\t'); - for j2=1:min([length(content),500]) - content2=content{j2}; - err=fprintf(fid,'{%s} ',content{j2}); + 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(:)'); end - elseif strcmp(z_whos(j).class,'struct') - content=fieldnames(eval(z_whos(j).name)); - content=char(content); content(:,end+1)=' '; content=content'; - err=fprintf(fid,'\t%s',content(:)'); - end - err=fprintf(fid,'\n'); - end, end - err=fclose(fid); + fprintf(fid,'\n'); + end + end + fclose(fid); disp('----------------------------'); disp(' ERROR OCCURED'); disp(' during executing mi'); @@ -663,18 +671,21 @@ catch 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,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 if nargout, varargout={NaN}; end - try, set(0,props.root), end + try set(0,props.root),catch end set(0,'ShowHidden','off') end