clc
clear all
load('pre_product.mat');
load('pre_gauge.mat');
for m=1:51;%%(number of gauge stations) change accordingly
aa=0;
bb=0;
cc=0;
dd=0;
pre_g_hit=[];
pre_p_hit=[];
pre_g_false=[];
pre_p_false=[];
pre_g_miss=[];
pre_p_miss=[];
for k=1:365 %%(days, since I used daily precip data) change accordingly
if (pre_gauge(m,k)>=0.1)&&(pre_product(m,k)>=0.1)
aa=aa+1;
pre_g_hit=[pre_g_hit;pre_gauge(m,k)];
pre_p_hit=[pre_p_hit;pre_product(m,k)];
elseif (pre_gauge(m,k)<0.1)&&(pre_product(m,k)>=0.1)
bb=bb+1;
pre_g_false=[pre_g_false;pre_gauge(m,k)];
pre_p_false=[pre_p_false;pre_product(m,k)];
elseif (pre_gauge(m,k)>=0.1)&&(pre_product(m,k)<0.1)
cc=cc+1;
pre_g_miss=[pre_g_miss;pre_gauge(m,k)];
pre_p_miss=[pre_p_miss;pre_product(m,k)];
elseif (pre_gauge(m,k)<0.1)&&(pre_product(m,k)<0.1)
dd=dd+1;
end
end
r0=corrcoef(pre_product(m,:),pre_gauge(m,:),'rows','complete');
Corr(m)=r0(2);
POD(m)=aa/(aa+cc); % this calculates the pod at each stations
FAR(m)=bb/(aa+bb);
ACC(m)=(aa+dd)/(aa+bb+cc+dd);
CSI(m)=aa/(aa+bb+cc);
FBI(m)=(bb+aa)/(aa+cc);
TS(m)=((aa.*dd)-(bb.*cc))/((aa+cc).*(bb+dd));
BIAS(m)=(aa+bb)/(aa+cc);
HSS(m)=2*((aa*dd)-(bb*cc))/[(aa+cc)*(cc+dd)+(aa+bb)*(bb+dd)]
AA(m)=aa;
BB(m)=bb;
CC(m)=cc;
DD(m)=dd;
end
ALL_POD=AA/(AA+CC) % average of all staitons' pod
ALL_FAR=BB/(AA+BB)
ALL_ACC=(AA+DD)/(AA+BB+CC+DD)
ALL_FBI=(AA+BB)/(AA+CC)
ALL_CSI=AA/(AA+BB+CC)
ALL_TS=((AA.*DD)-(BB.*CC))/((AA+CC).*(BB+DD));
ALL_BIAS=(AA+BB)/(AA+CC);
ALL_HSS=2*((AA.*DD)-(BB.*CC))/[(AA+CC).*(CC+DD)+(AA+BB).*(BB+DD)]
