复件 deternum1.m

AIC AND MDL 算法

function [num,seq]=deternum1(w,numl,type)
%依据特征值用(比率)AIC/MDL方法判定目标数目
%用法：
%	[num,seq]=deternum(w,numl,type)
%  或
%	num=deternum(w,numl,type)
%输入参数：
%  w	=== 升序排列的特征值(函数内部会重新排序)
% numl  === 形成相关矩阵所用的快拍数
% type  === 判阶方法
%	        =1	AIC方法；=2 MDL方法 =3 EDC
%	        =11~19 用比率准则，个位数表示方法(同上)
%	        =21~29 用斜率变化率准则，个位数表示方法(同上)
%输出参数：
%  num  === 判定的目标数目
%  seq  === 判定依据的序列

%   作者: 张群飞
%   Copyright 2006-2015, 西北工业大学
%   $Revision: 1.1 $  $Date: 2006/8/7 $
C=1;
amul=1; sum=0; iexp=0;sum1=0;
w=sort(w); % 强制特征值升序排列
k=length(w);  AIC=zeros(k-1,1);
for ii=1:k
    if w(ii)<=0 & w(ii)>-5.e-6,  w(ii)=1.e-7; end
    j=fix( log10(w(ii)) ); temp=1;
    if j~=0,  temp=10.^(-j); end
    amul=amul*w(ii)*temp;
    iexp=iexp+j;
    sum=sum+w(ii);
    sum1=sum1+log10(w(ii));
    AIC(ii)=-ii*numl*log( amul^(1/ii)*10^(iexp/ii)*ii / sum );
    switch mod(type,10)
        case 1
            AIC(ii)=AIC(ii)+(k-ii)*(k+ii+1);		   % AIC
        case 2
            AIC(ii)=AIC(ii)+(k-ii)*(k+ii+1)*log10(numl)/2.;    % MDL
        case 3
            AIC(ii)=-2*ii*num1*(sum1/ii-log10(sum/ii))+((k-ii)*(k+ii+2)+1)*C;
        otherwise
            disp(['指定的判阶方法',num2str(type),'不存在！'])
            return
    end
end
AIC(1:k)=AIC(k:-1:1);
switch fix(type/10)
    case 1,			% Ratio methods
        for ii=length(AIC):-1:2,
            AIC(ii)=AIC(ii)/AIC(ii-1);
        end
        AIC(1)=1;
    case 2,   % 用斜率变化率准则，只适用于目标数目不小于2的情况
        AIC=[0;rad2deg(atan(diff(AIC/max(AIC)*length(AIC))))];
        AIC=-diff(AIC);
end
seq=AIC;
[temp,num]=min(AIC);  % Take the global minimum as the number
num=num-1;