mfbox_mdl.m

来自「toolbox for spm 5 for data, model free a」· M 代码 · 共 201 行

function [K,erg,vK]=mfbox_mdl(e,N,gamma,verbose)
%MDL - find number of principal components such that it is optimal
%      with respect to the minimum description length
%      first in the input come the expected signal components 
%
% [K,erg,vK] = mdl(e,N,gamma,verbose);
%
%  e       - signal strengths (row vector)
%  N       - number of samples
%  gamma   - gamma in MDL estimator (>0 mdl(vetter), =0 aic, <0 mdl)
%
%  erg     - array of mdl values
%  K       - argmin(erg(i))
%  vK      - erg(K)
%
% 1.8.03 - initial version
% 1.9.04 - rewrite
% Peter Gruber
%
% Copyright by Peter Gruber and Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

error(nargchk(2,4,nargin))     % check no. of input args
if (nargin<3), gamma = 32; end
if (nargin<4), verbose = 0; end

if (gamma>0), mdltype = 'p';
elseif (gamma==0), mdltype = 's'; 
elseif (isnan(gamma)), mdltype = 'a'; 
else mdltype = 'r';
end;
gamma = abs(gamma);
THRESH = 10^(-3); % or sqrt(eps);

if (size(e,2)==1), e = e'; end

% normalize the signal strengths
e(e<0)=0;
for i=1:size(e,1)
    e(i,:) = size(e,2)*real(e(i,:)/norm(e(i,:)));
end

m = size(e,2);

maxi = find(any(e<THRESH,1),1)-1;
if (isempty(maxi)), maxi = m; end
% trivial cases
if (m==1 || maxi==0), K = 1; erg = K; return; end
if (sum(e)==0), K = m; erg = ones(m)*K; return; end

e = e(1:min(maxi,size(e,1)),1:maxi);
loge = log(e);

% initialize mdl
erg = ones(3,maxi)*NaN;

% get mdl
M = ((1:maxi)*maxi)-(((1:maxi).*((1:maxi)-1))/2)+1;
coeffs = getfittingcoeffs(maxi,N);

erg(1,maxi) = 0;
r = zeros(1,maxi-1);
for i=1:(maxi-1)
    if (size(e,1)==1), f = loge(1:maxi);
    else f = loge(i,1:maxi);
    end
    j = maxi-i;
    [c,r(i)] = getnoisefitting(f,j,coeffs);
    f = f-c;
    f = f((i+1):maxi);
    erg(1,i) = -N*j.*(sum(f)./j-log(sum(exp(f))./j));
end

erg(2,:) = M*log(N)/2;

switch mdltype
    case 'p'
        f = cumsum(loge,2)-[loge(:,2:end),repmat(0,size(loge,1),1)].*repmat(1:size(loge,2),size(loge,1),1);
        f(f>gamma) = 0; f(:,end) = 0;
        [i,j] = max(f,[],2);
        f = loge;
        for i=1:size(f,1), f(i,:) = f(i,:)-f(i,j(i)+1); end
        f(f<=0) = 0;
        if (size(f,1)==1)
            erg(3,:) = -M.*(cumsum(f)./i);
        else
            x = tril(f);
            erg(3,:) = -M.*(sum(x,2)'./i);
        end
        erg(3,:) = erg(3,:)+M*gamma;
    case 'r'
        erg(3,:) = M.*gamma;
    case 'a'
        erg(3,:) = M*log(N)/2;
    otherwise
        erg(3,:) = 0;
end

e = sum(erg(:,1:maxi),1);
serg = ones(1,m)*NaN;
serg(1:maxi) = e;
% find minimum
[vK,K] = min(serg);

if (verbose>1)
    warning(['MDL: ' sprintf('%d ',N,maxi) '-> ' sprintf('%f ',serg)]);
    if (verbose>7), cplot('combined',erg(1,:)/max(abs(erg(1,:))),erg(2,:)/max(abs(erg(2,:))),erg(3,:)/max(abs(erg(3,:)))); title('MDL parts'); pause; end
    if (verbose>5), cplot('combined',e/max(e),serg/max(serg)); title('MDL with e'); pause; end
end
erg = serg;

return


function f=getfittingcoeffs(T,N)

p = mfbox_mdlcoeffs();
f = evalpolcoeffs3(p,T,log(N),0,1);


function [c,r]=getnoisefitting(oc,T,f)

l = length(oc);
oc = oc-mean(oc);
%e = getmdlnoisecoeffs(oc,T);
is = repmat((0:(T-1))/(l-1),size(f,1),1); % at positions
for i=1:(size(f,1))
    is(i,:) = is(i,:).^(size(f,1)-i);
end
pe = is'*f; % evaluated
voc = oc(l:-1:(l-T+1))';
mpe = sum(pe-[zeros(size(pe,1),size(pe,2)-1),voc],1)/size(pe,1);
pe(:,3) = pe(:,3)-voc;
pe = pe-repmat(mpe,size(pe,1),1); % remove mean distance
p = zeros(1,2*size(f,2)-1);
for i=1:size(pe,1)
    p = p+conv(pe(i,:),pe(i,:)); %squared
end
if (all(abs(p(:))<eps))
    lp = 0.05;
else
    lp = [0.05;0.7;roots(polyder(p))]; % estimation only stable between 0.05 and 0.7
    lp = real(lp);
end
lp = lp(lp>=0.05 & lp<=0.7);
a = polyval(p,lp);
[v,r] = min(a);
r = lp(r);
p = zeros(1,size(f,1));
for i=1:size(f,1)
    p(i) = polyval(f(i,:),r);
end
if (l==1)
    c = polyval(p,0.5);
else
    c = polyval(p,((l-1):-1:0)/(l-1));
end


function v=evalpolcoeffs3(pc,a,b,c,n)

v = zeros(size(pc,1),size(pc,3),size(pc,4));
for p=1:size(pc,1)
    for i=1:size(pc,3)
        for j=1:size(pc,4)
            x = pc(p,:,i,j);
            v(p,i,j) = polyval(x(:),a);
        end
    end
end
if (n>0), v = evalpolcoeffs2(v,b,c,n-1); end 


function v=evalpolcoeffs2(pc,a,b,n)

v = zeros(size(pc,1),size(pc,3));
for p=1:size(pc,1)
    for i=1:size(pc,3)
        x = pc(p,:,i);
        v(p,i) = polyval(x(:),a);
    end
end
if (n>0), v = evalpolcoeffs1(v,b,n-1); end


function v=evalpolcoeffs1(pc,a,n)

v = zeros(size(pc,1),1);
for p=1:size(pc,1)
    x = pc(p,:);
    v(p) = polyval(x(:),a);
end