function BEST_DATA=MISSING_classification(R,IHD,time_series_names)


num_params=size(R,2);
DATA=zeros(size(R,1),3);

include=zeros(size(R,2),1);









for par_idx=1:num_params
    
    X=R(:,par_idx);
    B = mnrfit(X,IHD+1);
    
    PHAT = mnrval(B,X);
    
    
    
    
    max_score1=0;
    max_score2=0;
    best_th=0;
    
    for class_th=.1:.01:0.5
        
        
        
        DATA(:,1)=str2double('0000');
        
        DATA(:,2)=PHAT(:,2);
        
        DATA(:,3)=PHAT(:,2)> class_th;
        
        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;
                max_score2=H;
                %  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);
   
   if max_score1>0
       
      % include(par_idx)=1;
       
   disp(' ---')
   
   disp([' non NaNs:' num2str(deaths1/(deaths1+survivals1)) ])
   disp([' NaNs:' num2str(deaths2/(deaths2+survivals2)) ])
   
   ratio=(deaths1/(deaths1+survivals1)) / (deaths2/(deaths2+survivals2))
   
  name= time_series_names{par_idx}
  score= max_score1
   
   disp(' ---')
   
     if ratio >1.5 || ratio <0.6
       
       include(par_idx)=1;
     end
   
   end
    
end


num_params=size(R,2);
DATA=zeros(size(R,1),3);

%for par_idx=1:num_params


    X=R(:,include==1);
    
    
    
% max_N=size(X,2);
% X=X(:,1:max_N);
% 
% X_orig=X;
% 
% scores=zeros(max_N,1);
% 
% for idx=10:max_N
% 
%     X=X_orig(:,1:idx);
% %[best_th0,max_score_0,D] = opt_th(X,IHD);
% % 
% % [best_th1,max_score_1,D0] = opt_th(X(Gender_a==0,:),IHD(Gender_a==0));
% % [best_th2,max_score_2,D1] = opt_th(X(Gender_a==1,:),IHD(Gender_a==1));
% 
% [best_th0,max_score_0,D] = opt_th_valid(X,IHD);
% % [best_th1,max_score_1,D0] = opt_th_valid(X(Gender_a==0,:),IHD(Gender_a==0));
% % [best_th2,max_score_2,D1] = opt_th_valid(X(Gender_a==1,:),IHD(Gender_a==1));
% 
% scores(idx,:)=[max_score_0];
% end
% 
% figure(1)
% plot(scores(:,1),'g')
% % hold on
% % plot(scores(:,2),'r')
% % plot(scores(:,3),'b')
% pause(.1)
    
    
    %0.2708
    %X=[R(:,~strcmp(time_series_names,'MechVent'))];   
    
    %0.277
%     X=[R(:,strcmp(time_series_names,'TroponinT'))];    
%     X=[R(:,strcmp(time_series_names,'TroponinI')) X];
%     X=[R(:,strcmp(time_series_names,'RespRate')) X];
%     X=[R(:,strcmp(time_series_names,'Mg')) X];
%     X=[R(:,strcmp(time_series_names,'Lactate')) X];
%     X=[R(:,strcmp(time_series_names,'K')) X];
%     X=[R(:,strcmp(time_series_names,'FiO2')) X];
%     X=[R(:,strcmp(time_series_names,'Bilirubin')) X];
%     X=[R(:,strcmp(time_series_names,'AST')) X];
%     X=[R(:,strcmp(time_series_names,'ALT')) X];
%     X=[R(:,strcmp(time_series_names,'ALP')) X];
%     X=[R(:,strcmp(time_series_names,'Albumin')) X];
%     
    B = mnrfit(X,IHD+1);
    
    PHAT = mnrval(B,X);
    
    
    
    
    max_score1=0;
    max_score2=0;
    best_th=0;
    
    for class_th=.1:.01:0.5
        
        
        
        DATA(:,1)=str2double('0000');
        
        DATA(:,2)=PHAT(:,2);
        
        DATA(:,3)=PHAT(:,2)> class_th;
        
        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;
                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)) ])
%    
   
   max_score1
   

function [best_th,max_score1,BEST_DATA]=opt_th_valid(X_,IHD_)
        
  %  num_params=size(time_series_names,1);
    
    
    n=size(X_,1);
    s=zeros(n,1);
    s(1:round(n/2))=1;
    s=boolean(s);
    B = mnrfit(X_(s,:),IHD_(s)+1);
    
    PHAT = mnrval(B,X_(~s,:));
    
    BEST_DATA=[];
    DATA_=zeros(size(PHAT,1),3);
    
    max_score1=0;
    max_score2=0;
    best_th=0;
    
    for class_th=.1:.01:0.5
        
        
        
        DATA_(:,1)=str2double('0000');
        
        DATA_(:,2)=PHAT(:,2);
        
        DATA_(:,3)=PHAT(:,2)> class_th;
        
        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_(~s)==1,3));
            FN=sum(~DATA_(IHD_(~s)==1,3));
            FP=sum(DATA_(IHD_(~s)==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;
               % 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)) ])
%    
%    
%    max_score1
   
end

end