mfbox_cordes.m

toolbox for spm 5 for data, model free analysis

function [K,erg,vK]=mfbox_cordes(e,N,tol,verbose)
%MFBOX_CORDES - find number of principal components such that it is optimal
%      with respect to the criterion introduced by cordes et al.
%
% [K,erg] = cordes(e,N,tol,verbose);
%
%  e       - signal strengths (row vector)
%  N       - number of samples
%  tol     - tolerance
%
%  erg     - noise variance estimate
%  K       - min(|erg(i)-e(i)|<(tol/e(i)))
%
% 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), tol = 0.1; end
if (nargin<4), verbose = 0; end

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);

coeffs = getfittingcoeffs(maxi,N);

loge = log(e);

vK = ones(1,maxi);
erg = zeros(maxi,maxi);
for i=1:(maxi-1)
    if (size(e,1)==1), f = loge(1:maxi); o = e(1:maxi);
    else f = loge(i,1:maxi); o = e(i,1:maxi);
    end
    j = maxi-i;
    c = getnoisefitting(f,j,coeffs);
    f = exp(c);
    erg(i,:) = f;
    l = max(find((abs(f-o)./o)>=tol));
    vK(i) = l;
end

K = min(find(vK<=[1:length(vK)]));

if (verbose>1)
    warning(['CORDES: ' sprintf('%d ',N,maxi) '-> ' sprintf('%f ',vK)]);
end

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
c = c-polyval(mpe,r);


function p=getmdlnoisecoeffs(e,T)

l = length(e);
e = e(l:-1:(l-T+1));
v = min(8,length(e));
T = T-1;
p = polyfit((0:T)/l,e,v-1);
if (v<8), p = [zeros(1,8-v),p]; 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


function e=correctvariance(e,N)

T = length(e);
x = linspace(1,-1,T);
x = asin(x)+(x.*sqrt(1-x.^2));
x = (2*x*T)/(N*pi);
e = e.*(1+x);


function [e,D]=simnoisevar(N,r,T)

D = zeros(T,T);
for j=1:T, D = D+j*diag(ones(1,T-j+1),j-1); end
D = r.^(D+D'-diag(diag(D))-1);
v = real(eig(D));
e = sort(v,'descend');
e = correctvariance(e',N)';
e(e<=0) = 0;


function [e,D]=getnoisevar(N,r,T)

D = zeros(N,T);
D(:,1) = randn(N,1);
for i=2:T, D(:,i) = (r*D(:,i-1)+randn(N,1))/sqrt(1+r^2); end
v = real(eig(cov(D)));
v(v<=0) = 0;
e = sort(v,'descend');


function [Vars,Params,Polcoeffs]=genmdlcoeffs()

Params = zeros(0,3);
Vars = struct();
l = 1;
for i=1:5
    for j=1:30
        for k=1:10
            Vars(l) = struct();
            l = l+1;
        end
    end
end
l = 1;
for i=1:5
    for j=1:30
        for k=1:10
            T = 10*j;
            N = 5*10^(i+1);
            r = 0.05*(k-1);
            v = sprintf('%d,%d,%f',T,N,r);
            Params(l,:) = [T,N,r];
            Vars(l).name = v;
            Vars(l).vars = simnoisevar(N,r,T)';
            l = l+1;
        end
    end
end
l = 1;
for i=1:5
    for j=1:30
        for k=1:10
            T = 10*j;
            N = 10^(i+1);
            r = 0.1*(k-1);
            v = real(Vars(l).vars);
            v(v<=0) = min(v(v>0));
            v = real(log(v));
            w = ((T-1):-1:0)/(T-1);
            Vars(l).pol3 = polyfit(w,v,3);
            Vars(l).pol5 = polyfit(w,v,5);
            Vars(l).pol7 = polyfit(w,v,7);
            l = l+1;
        end
    end
end

L = (cat(1,Vars.pol7));
Lpot = reshape(L,10,30,5,8);
Lpar = reshape(Params,10,30,5);
Lpar = Lpar(:,:,:,[3,1,2]);
sx = [6,5,3];
Lpp = zeros(8,30,5,sx(3));
for p=1:8
    for i=1:30
        for j=1:5
            Lpp(p,i,j,:) = polyfit(Lpar(:,i,j,1),Lpot(:,i,j,p),sx(3)-1);
        end
    end
end
Lppp = zeros(8,30,sx(2),sx(3));
for p=1:8
    for i=1:30
        for j=1:sx(3)
            Lppp(p,i,:,j) = polyfit(log(Lpar(1,i,:,3)(:)),Lpp(p,i,:,j)(:),sx(2)-1);
        end
    end
end
Lpppp = zeros(8,sx(1),sx(2),sx(3));
for p=1:8
    for i=1:sx(2)
        for j=1:sx(3)
            Lpppp(p,:,i,j) = polyfit(Lpar(1,:,1,2)(:),Lppp(p,:,i,j)(:),sx(1)-1);
        end
    end
end
Polcoeffs = Lpppp;