function [best_th,max_score1,BEST_DATA]=opt_th_multi(D1,D2,D3,D4,IHD_)

% num_params=size(time_series_names,1);
DATA=zeros(size(IHD_,1),3);
DATA2=zeros(size(IHD_,1),3);
%B = mnrfit(X_,IHD_+1);

PHAT = D1(:,1);
PHAT2 = D2(:,1);
PHAT3 =  D3(:,1);
PHAT4 =  D4(:,1);

max_score1=0;
max_score2=0;
best_th=0;

for class_th=.1:.01:0.5
    for class_th2=.1:.01:0.5
        for class_th3=.1:.01:0.5
            for class_th4=.1:.01:0.5

                DATA(:,1)=str2double('0000');

                DATA(:,2)=PHAT(:,1);

                %0.4991
                DATA(:,3)=PHAT(:,1)> class_th | PHAT2(:,1)> class_th2 | PHAT3(:,1) > class_th3| PHAT4(:,1) > class_th4 ;

                %0.4838
                %  DATA(:,3)=mode(double([PHAT(:,1)> class_th  PHAT2(:,1)> class_th2  PHAT3(:,1) > class_th3])')';

                DATA(DATA(:,2)<0.01,2)=0.01;
                DATA(DATA(:,2)>0.99,2)=0.99;


                %   if(~isempty(results))

                % Calculate sensitivity (Se) and positive predictivity (PPV)
                TP=sum(DATA(IHD_==1,3));
                FN=sum(~DATA(IHD_==1,3));
                FP=sum(DATA(IHD_==0,3));
                Se=TP/(TP+FN);
                PPV=TP/(TP+FP);

                show=0; % if show is 1, the decile graph will be displayed by lemeshow()
                % H=lemeshow([IHD_ DATA(:,2)],show);

                % Use the title of figure to display the results
                % title(['H= ' num2str(H) ' Se= ' num2str(Se) ' PPV= ' num2str(PPV) '. ' num2str(class_th) ])

                % The event 1 score is the smaller of Se and PPV.
                score1 = min(Se, PPV);
                if score1>max_score1
                    max_score1=score1;
                    best_th=[class_th class_th2 class_th3 class_th4];
                    % best_th2=class_th2;
                    %                max_score2=H;
                    BEST_DATA=DATA;
                    %  display(['Unofficial Event 1 score: ' num2str(score1)]);
                end

                % end

            end

            %    survivals1=sum(X==1 & IHD==0);
            %    deaths1=sum(X==1 & IHD==1);
            %    survivals2=sum(X==0 & IHD==0);
            %    deaths2=sum(X==0 & IHD==1);
            %
            %    disp([' non NaNs:' num2str(deaths1/(deaths1+survivals1)) ])
            %    disp([' NaNs:' num2str(deaths2/(deaths2+survivals2)) ])
            %


        end
    end
end
max_score1
end