diff options
Diffstat (limited to 'sourcecodes/bnt-master/SLP/misc/hist_ic.m')
| -rw-r--r-- | sourcecodes/bnt-master/SLP/misc/hist_ic.m | 268 |
1 files changed, 268 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/SLP/misc/hist_ic.m b/sourcecodes/bnt-master/SLP/misc/hist_ic.m new file mode 100644 index 00000000..b16980fd --- /dev/null +++ b/sourcecodes/bnt-master/SLP/misc/hist_ic.m @@ -0,0 +1,268 @@ +function [n,edges,nbedges,xechan] = hist_ic(x,crit) + +%HIST_IC optimal Histogram based on IC information criterion +% +% [N,EDGES,NBEDGES,XECHAN] = HIST-IC(X,CRIT) +% bins the elements of X into an optimal number of bins according +% to a cost function based on Akaike's Criterion. +% +% +% CRIT = 1 | 2 | 3 (choose one of the 3 possible criterium) (default 3) +% 4 (returns the initial histogram instead of the optimal one) +% +% +% N = cell array containing the distribution of each column of X +% (or a vector if X is a column vector) +% EDGES = cell array containing the bin edges of each column of X +% (or a vector if X is a column vector) +% NBEDGES = vector containing the number of bin edges for each column of X +% (or a number if X is a column vector) +% XECHAN = discretized version of X +% +% Ref : O. Colot et al., Information Criteria and Abrupt Changes in +% Probability Laws, Signal Processing VII: Theory and Applications +% pp.1855-18858, September 1994 +% +% F. El-Matouat, O. Colot 2000 (first version) +% Revised 01-06-2001 by Ph. Leray - philippe.leray@univ-nantes.fr +% +% +% Things to do : +% * Call criteron by a name ('aic','xxx', ...) instead of a number +% + + +if nargin == 0 + error('Requires one or two arguments.') +end + +if nargin == 1 + crit = 3; +end; + +if min(size(x))==1, x = x(:); end + +if isstr(x) + error('Input argument must be numeric.') +end + +if isempty(x), + error('No elements to count') +end + + +[nb_l,nb_c]=size(x); + +% Outputs declaration +xechan=zeros(nb_l,nb_c); + +edges=cell(nb_c,1); +% Local variables +maxi = max(x); +mini = min(x); + +%% Erreur ? ancien code : +%% nb_clas_ini=2*round(sqrt(nb_l)-1); % article Fatima + +nb_clas_ini=round(2*sqrt(nb_l)-1); + +pas_ini=(maxi-mini)/nb_clas_ini; % initial step + +for j=1:nb_c, + + % optimal histogram for each column of X + histo_ini =hist(x(:,j),nb_clas_ini); % initial histogram + + if (crit~=4) + [hist_opt,pas_opt]=hist1_ic(histo_ini,nb_l,pas_ini(j),nb_clas_ini,crit); + else + fprintf('Histo initial\n'); + hist_opt=histo_ini; + pas_opt=ones(1,nb_clas_ini)*pas_ini(j); + end; + nbedges(j)=size(hist_opt,2); + edges{j}=mini(j)+cumsum(pas_opt(1:nbedges(j)-1)); %+1e-7; + [n{j} xechan(:,j)]=histc(x(:,j),[-inf edges{j} inf]); + n{j}=n{j}(1:end-1); +end + +if (nb_c==1) + n=n{1}; edges=edges{1}; +end + + +% ============================== subfunctions + +function [hist_opt,step_opt]=hist1_ic(histo,nb,step_ini,m,critere); + +%HIST1_IC optimal Histogram based on IC information criterion +% +% [HIST_OPT, STEP_OPT] = HIST1_IC(HISTO, NB, STEP_INI, NBSTEP_INI, CRIT) +% fusion of an 1D histogramme (HISTO) according to an IC criterion (CRIT) +% +% This function is mainly an internal function used by HIST_IC +% +% Ref : O. Colot et al., Information Criteria and Abrupt Changes in +% Probability Laws, Signal Processing VII: Theory and Applications +% pp.1855-18858, September 1994 +% +% F. El-Matouat, O. Colot 2000 (first version) +% Revised 11-06-2001 by Ph. Leray +% +% +% Things to do : +% * Call criteron by a name ('aic','xxx', ...) instead of a number +% + + +aic=[]; +aic2=[]; + +% Initialisation +histt=histo; +teta = histt/nb; +pas = step_ini*ones(1,m); + +% Calcul de l'ensemble des histogrammes optimaux + +for z=1:m + + % Calcul de AIC pour l'union entre hist(indice,u) et hist(indice,u+1) + aic2 = [aic2 cal_aic(nb,teta,pas,m+1-z,critere)]; + + if (z~=m) + % Calcul des couples de classes adjacentes + if critere==1 + penalite=(2*(m-z)-1)/nb; + elseif critere==2 + penalite=(m-z-1)*(1+log(nb))/nb; + else + penalite=(m-z)*(1+log(log(nb)))/nb; + end + aic=cla_adj(nb,teta,pas,step_ini,m-z+1,histt,penalite,aic); + % Recherche de la valeur min du crit�re pour les classes adjacentes + [min_aic classe]=min(aic(1:(m-z))); + + + % Fusion de hist(classe) et hist(classe+1) + nb_pas1=pas(classe)/step_ini; + nb_pas2=pas(classe+1)/step_ini; + + ess=round( nb_pas1*histt(classe)+nb_pas2*histt(classe+1) ); + teta(classe)=ess / nb; + histt(classe)=ess / (nb_pas1+nb_pas2); + pas(classe)=pas(classe)+pas(classe+1); + + + % Cr�ation du nouvel histogramme + itemp = setdiff(1:m+1-z,classe+1); + histt = histt(itemp); + pas = pas(itemp); + teta = teta(itemp); + end +end + +% Recherche du crit�re minimun AIC +[min_AIC fusion]=min(aic2(1:m)); + +% Initialisation de histo +histt=histo; +teta = histt/nb; +pas = step_ini*ones(1,m); + +% Calcul de l'histogramme optimal + +for z=1:fusion-1 + + % Calcul des couples de classes adjacentes + if critere==1 + penalite=(2*(m-1)-1)/nb; + elseif critere==2 + penalite=(m-2)*(1+log(nb))/nb; + else + penalite=(m-1)*(1+log(log(nb)))/nb; + end + + aic=cla_adj(nb,teta,pas,step_ini,m,histt,penalite,aic); + + + % Recherche de la valeur min du crit�re pour les classes adjacentes + [min_aic classe]=min(aic(1:m-1)); + + % Fusion de hist(indice,classe) et hist(indice,classe+1) + + nb_pas1=pas(classe)/step_ini; + nb_pas2=pas(classe+1)/step_ini; + + teta(classe)=(round(nb_pas1*histt(classe)+nb_pas2*histt(classe+1)))/nb; + histt(classe)=(nb_pas1*histt(classe)+nb_pas2*histt(classe+1))/(nb_pas1+nb_pas2); + pas(classe)=pas(classe)+pas(classe+1); + + % Cr�ation du nouvel histogramme + + itemp=setdiff(1:m,classe+1); + histt = histt(itemp); + pas = pas(itemp); + teta = teta(itemp); + %aic=zeros(1,m-1); + + m=m-1; +end +hist_opt=histt; +step_opt=pas; + + +%===================================================== +% Calcul du Critere pour l'ensemble des classes + +function akaike=cal_aic(size_ech,teta,pas,m,critere); + + +if critere==1 + a=(2*m-1)/size_ech; +elseif critere==2 + a=(m-1)*(1+log(size_ech))/size_ech; +else + a=m*(1+log(log(size_ech)))/size_ech; +end + +indu = find(teta); +akaike = a - 2*sum(teta(indu).*log(teta(indu)./pas(indu))); + + +%===================================================== +% Cla_adj.m +% aic=cla_adj(taille,indice,teta,pas,pas_ini,m,hist,penalite,aic) +% taille=nombre d'�l�ments dans chacune des hypotheses; +% indice=numero de la classe; +% Calcul du critere de Akaike pour l'histogramme totale avec +% fusion de deux classes adjacentes u et (u+1). + +function aic=cla_adj(size_ech,teta,pas,pas_ini,m,hist,penalite,aic); + + +for u=1:m-1 + + cumul=0; + + % This loop is faster than a sum of a vectorised computation ! + for x=1:m + if x~=u & x~=u+1 & teta(x)~=0 + cumul=cumul+teta(x)*log(teta(x)/pas(x)); + end + end + + + nb_pas1=pas(u)/pas_ini; + nb_pas2=pas(u+1)/pas_ini; + + b=( round(nb_pas1*hist(u)+nb_pas2*hist(u+1) ) ) / size_ech; + + if b~=0 + c=2*b*log( b / ( pas(u) + pas(u+1) ) ); + else + c=0; + end + + aic(u)=penalite-2*cumul-c; +end |