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)]
