function [mobs,mest,sdobs,sdest,R2,imae,rmse,E,Elog,E1,d,d1] = goodfit2(obs,est)

% this is a function I developed to calculate

% different indices for checking the goodness-of-fit

% of any model (not valid for binary outputs)

%

% inputs are:

% est - estimated values by the model (should be a vector (ntv,1))

% obs - values the model should be compared to (should be a vector (ntv,1))

%

% outputs are:

% R2 - square of the correlation btw observed and estimated values

% MAE - mean absolut error

% RMSE - root mean square error

% E - coefficient of efficiency - Nash and Sutcliffe (1970)

% Elog - log Nash-Sutcliffe coefficient of efficiency

% E1 - modified Nash-Sutcliffe coefficient of efficiency

% d - index of agreement (Willmott 1981)

% d1 - modified index of agreement

[nz,ntv] = size(obs);

% means

disp(' ')

mobs = mean(obs);

disp (['mean observed values = ' num2str(mobs)])

mest = mean(est);

disp (['mean estimated values = ' num2str(mest)])

% standart deviations

sdobs = std(obs);

disp (['std dev observed values = ' num2str(sdobs)])

sdest = std(est);

disp (['std dev estimated values = ' num2str(sdest)])

% assessment of goodness-of-fit

% correlation and coefficient of determination (explained variance)

[corre,p] = corrcoef([obs est]);R2=corre.^2;

disp (['R2 corre^2 = ' num2str(R2(1,2)) ' p = ' num2str(p(1,2))])

% mean absolute error

imae = sum(abs(obs-est))/ntv;

disp (['mae = ' num2str(imae)])

% root mean square error

rmse = sqrt(sum((obs-est).^2)/ntv);

disp (['rmse = ' num2str(rmse)])

disp(' ')

% Nash-Sutcliffe coeff of efficiency

E = 1-(sum((obs-est).^2)/sum((obs-mobs).^2));

disp (['N-Sut efficiency (E) = ' num2str(E)])

stat = bootstrp(500,@Nash_Sut,obs,est);

[a b]=ci95(stat);

disp ([ 'E 95% ci = ' num2str(a) ' ' num2str(b)])

disp (' ');clear stat

% log Nash-Sutcliffe coeff of efficiency

Elog = 1-(sum((log(obs)-log(est)).^2)/sum((log(obs)-log(mobs)).^2));

disp (['Elog = ' num2str(Elog)])

stat = bootstrp(500,@log_Nash_Sut,obs,est);

[a b]=ci95(stat);

disp ([ 'E 95% ci = ' num2str(a) ' ' num2str(b)])

disp (' ');clear stat

% Modified Nash-Sutcliffe coeff of efficiency

E1 = 1-(sum(abs(obs-est))/sum(abs(obs-mobs)));

disp (['N-Sut efficiency modified (E1) = ' num2str(E1)])

stat = bootstrp(500,@m_Nash_Sut,obs,est);

[a b]=ci95(stat);

disp ([ 'E1 95% ci = ' num2str(a) ' ' num2str(b)])

disp (' '); clear stat

% index of agreement

d = 1-(sum((obs-est).^2)/sum((abs(obs-mobs)+abs(obs-mobs)).^2));

disp (['Index of agreement (d) = ' num2str(d)])

stat = bootstrp(500,@I_agree,obs,est);

[a b]=ci95(stat);

disp ([ 'd 95% ci = ' num2str(a) ' ' num2str(b)])

disp (' '); clear stat

% Modified index of agreement

d1 = 1-(sum(abs(obs-est))/sum((abs(obs-mobs)+abs(obs-mobs))));

disp (['Index of agreement modified (d1) = ' num2str(d1)])

stat = bootstrp(500,@m_I_agree,obs,est);

[a b]=ci95(stat);

disp ([ 'd1 95% ci = ' num2str(a) ' ' num2str(b)])

disp (' '); clear stat