function [tpfpReg, tpfpOLS, ROC] = tpfpcurves(D)

D.Reg = squeeze(D.Reg);


for ii = 1:size(D.Reg,3)
    k=0;
    
    for q = 0:0.025:1
        k = k+1;
        Dreal = D.Struct;
        Dhat = D.estimated;
        DReg = D.Reg(:,:,ii);
        
%         Dreal(Dreal<quantile(Dreal(:),q)) = 0;
        P = sum(sum(Dreal>0));
        N2 = sum(sum(Dreal<=0));
        DReg(DReg<quantile(DReg(:),q)) = 0;
        Dhat(Dhat<quantile(Dhat(:),q)) = 0;
        
        Test = DReg>0;
        TestOls = Dhat>0;
        Real = Dreal>0;
        
        Aux = Test-2*Real;
        
        
        tpfpReg(1,k,ii) = sum(sum(Aux==-1))/P;
        tpfpReg(2,k,ii) = sum(sum(Aux==1))/N2;
        
        Aux = TestOls-2*Real;
        
        tpfpOLS(1,k) = sum(sum(Aux==-1))/P;
        tpfpOLS(2,k) = sum(sum(Aux==1))/N2;
    end
    
end

if length(size(tpfpReg))>2
    for jj=1:size(tpfpReg,3)
        AUROC(jj)=-trapz(tpfpReg(2,:,jj),fliplr(tpfpReg(1,:,jj)));
    end
    [~,Imax]=max(AUROC);
    tpfpReg = tpfpReg(:,:,Imax);
end

%     AUROCols = -trapz(tpfpOLS(2,:),fliplr(tpfpOLS(1,:)));
% plot(tpfpReg(2,:),tpfpReg(1,:));
% hold on
% plot(tpfpOLS(2,:),tpfpOLS(1,:))
% hold off

end