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