[argout1 argout2 argout3 argout4]=sdindex2(data,threshold)

function [argout1 argout2 argout3 argout4]=sdindex2(data,threshold)
%calculation of statistical dependence index
%INPUTS
%data: a two column matrix, each column is a variable, each row an observation
%thresholds: reference thresholds for calculation of probability
%OUTPUTS
%statistical dependence index: SDI = Pab / Pa*Pb
%joint probability: Pab=Pa*Pa/b where Pa/b is conditional probability

if size(data,2)~=2
    disp('data matrix should contain two columns')
    return
end

if length(threshold)>1
    for tc=1:length(threshold)
    %find Pa, assuming site A column 1 of data
        Pa(tc)=length(find(data(:,1)>threshold(tc)))/length(data(:,1)); %wanted events/all possible events
      %find Pb, assuming site B column 2 of data
        Pb(tc)=length(find(data(:,2)>threshold(tc)))/length(data(:,2)); %wanted events/all possible events
       %find the conditional probability Pab
       Xa=data(:,1)*0;Xa(find(data(:,1)>threshold(tc)))=1;
       Xb=data(:,2)*0;Xb(find(data(:,2)>threshold(tc)))=1;
       Pab(tc)=sum(Xa.*Xb)/length(Xa);
       %Pab(tc)=length(find(and(data(:,1)>threshold(tc),data(:,2)>threshold(tc))))/size(data,1); %simultaneous exceedance of threshold/all possible outcomes
       %calculate index
        sdi(tc)=Pab(tc)/(Pa(tc)*Pb(tc));
    end
    
else
        Pa=length(find(data(:,1)>threshold))/length(data(:,1));
        Pb=length(find(data(:,2)>threshold))/length(data(:,2));
        Xa=data(:,1)*0;Xa(find(data(:,1)>threshold))=1;
        Xb=data(:,2)*0;Xb(find(data(:,2)>threshold))=1;
        Pab=sum(Xa.*Xb)/length(Xa);
        %Pab=length(find(and(data(:,1)>threshold,data(:,2)>threshold)))/size(data,1);
        sdi=Pab/(Pa*Pb);
end
    
switch nargout
    case 1
        argout1=sdi;
    case 2
        argout1=sdi;
        argout2=Pab;
    case 3
        argout1=sdi;
        argout2=Pab;
        argout3=Pb;
    case 4
        argout1=sdi;
        argout2=Pab;
        argout3=Pb;
        argout4=Pa;
    otherwise
        disp('error in number of specified output arguments')
end
return