%----- Compare 2 different models on the same dataset by Akaike's AICc  ----
% by Evgenii Kovrigin    3/27/08
% For details of model testing see  [1]
%
% [1].	Motulsky, H. and Christopoulos, A., Fitting Models to Biological Data Using 
% Linear and Nonlinear Regression: A Practical Guide to Curve Fitting. 1 edition ed. 2004: 
% Oxford University Press, USA. 352 
%
% Here we have SAME NUMBER OF POINTS and variable number of parameters
%
% To use: 
%      - Run  both models fitting for this dataset 
%         NOTE: make sure fits are reasonable!
%
%      - LineShapeKin saves the sum of resudual squares and number of points in
%         the files resname_SS.txt. Look up composite name for the global
%         fitting session in fit_results/  file.
%
%      -  put necessary file names into AIC_setup on local folder with
%      ...SS.txt files.
%
%     - enter  number of fitting parameter as well
%
%      -  execute 'compare.m'. Results are displayed and saved into plain
%      text file.

AIC_setup_2models


resname=deblank(resname);

%--- reading 1st fit result --- (former global)
 
fname=sprintf('%s.SS.txt',resname);
fp=fopen(fullfile(Folder_1,'results','SS',fname),'r');
data=fscanf(fp,'%d%e');
SS1=data(2);
N1=data(1);
disp(sprintf('%s\t%d\t%e\t\n\n',resname,N1,SS1));


%--- reading 2nd fit result --- (former individual)

fname=sprintf('%s.SS.txt',resname);
fp=fopen(fullfile(Folder_2,'results','SS',fname),'r');
data=fscanf(fp,'%d%e');
SS2=data(2);
N2=data(1);
disp(sprintf('%s\t%d\t%e\t\n\n',resname,N2,SS2));


%% Akaike criterion testing (see [1])
disp('Performing evaluation of likelyhood for each model by Akaike Information Criterion');


K2=Num_Fit_Pars_2 + 1;
K1=Num_Fit_Pars_1 + 1;


%-- check applicability
if K2 > (N2-2)
    error('Too few experimental points for Akaike test! N  must be at least 2 larger than K value!:\nK2=%f, N2=%f', K2, N2)
end
%-- check applicability
if K1 > (N1-2)
    error('Too few experimental points for Akaike test! N  must be at least 2 larger than K value!:\nK1=%f, N1=%f', K1, N1)
end



%--model 1:  AICc  
AICc_1= AICc(N1, K1, SS1);

%-- model 2: AICc  
AICc_2= AICc(N2, K2, SS2);


if AICc_1 > AICc_2
    conclusion='Model 2 is more correct';
else
    conclusion='Model 1 is more correct';
end


%-- Evidence ratio : probability of correct 1/ probability of correct 2
EvRat=1/exp(-0.5*(AICc_2-AICc_1));

%--- output for statistical model testing ----
fp=fopen(outfname,'W');
outstring1=sprintf('1st Model  \nN = %d\nK=%d\nSS = %g\nAICc=%g\n',N1, K1, SS1, AICc_1);
outstring2=sprintf('2nd Model  \nN = %d\nK=%d\nSS = %g\nAICc=%g\n',N2, K2, SS2, AICc_2);
outstring3=sprintf('\nEvidence ratio (model 1 is more likely than model 2 by factor of): %.1e\n',EvRat);
outstring4=sprintf('\n1/Evidence ratio (model 2 is more likely than model 1 by factor of): %.1e\n',1/EvRat);





disp(outstring1);disp(outstring2);disp(conclusion);disp(outstring3);disp(outstring4);
fprintf(fp,'%s\n%s\n%s\n%s\n%s ',outstring1,outstring2,conclusion,outstring3,outstring4);




fclose(fp);
disp(sprintf('Results are written to %s', outfname));