function plot_means_comparison2(X1,X2,IHD,ALL_CATEGORIES)

    num_params=size(X1,2);

    for param_idx=1:num_params
%         mean_all=X(:,param_idx);
%         mean_1=X(IHD==1,param_idx);
%         mean_0=X(IHD==0,param_idx);
%         
%         mean_1_abs=abs(X(IHD==1,param_idx));
%         mean_0_abs=abs(X(IHD==0,param_idx));
%         
%         isnan1=sum(~isnan(X(IHD==1,param_idx)))/(sum(IHD==1));
%         isnan0=sum(~isnan(X(IHD==0,param_idx)))/(sum(IHD==0));
%         
%         M1 = mode(mean_1);
%         M0 = mode(mean_0);
%         
        selector=~isnan(X1(:,param_idx))& ~isnan(X2(:,param_idx));
        IHD_Test=IHD(selector);
        
        x1=[X1(selector,param_idx) X1(selector,param_idx)-X2(selector,param_idx) ];
        x2=[X1(selector,param_idx) X2(selector,param_idx) ];
        x3=[X1(selector,param_idx) abs(X1(selector,param_idx)-X2(selector,param_idx)) ];
        
     
        
        max_score1_dif= th_search(x1,IHD_Test);
        max_score1_norm= th_search(x2,IHD_Test);
        max_score1_abs= th_search(x3,IHD_Test);
        
        
        
        if max_score1_abs >max_score1_dif
            ALL_CATEGORIES{param_idx} 
            max_score1_dif
            max_score1_norm
            max_score1_abs
        end
%         figure(param_idx)
%         
%         subplot(1,3,1)
%         hist(mean_all)
%         title([ALL_CATEGORIES{param_idx} ' all'])
%         
%         
%         subplot(1,3,2)
%         hist(mean_1)   
%         title([ALL_CATEGORIES{param_idx} ' 1'])
%        
%         
%         subplot(1,3,3)
%         hist(mean_0)    
%         title([ALL_CATEGORIES{param_idx} ' 0'])
        
%         figure(param_idx+num_params+5)
%         hist(mean_1,20)
%         h = findobj(gca,'Type','patch');
%         set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75)
%         hold on;
%         hist(mean_0,20)
%         h1 = findobj(gca,'Type','patch');
%         set(h1,'facealpha',0.75);
%         title([ALL_CATEGORIES{param_idx} ' mode0=' num2str(M0) '. '  ' mode1=' num2str(M1)]);
% 
%      
% % %         
% % %         figure(param_idx+num_params+5)
% % %         title([ALL_CATEGORIES{param_idx} ' norm=' num2str(max_score1_norm) '. '  ' abs= ' num2str(max_score1_abs)  ]);
% % % 
% % %         subplot(1,2,1)
% % % %         [N, X] =  hist(mean_1,25);
% % % %         bar(X, N./sum(N), 1);
% % % %         h = findobj(gca,'Type','patch');
% % % %         set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75)
% % % %         hold on;
% % % %         [N, X] =  hist(mean_0, X);
% % % %         bar(X, N./sum(N), 1);
% % % %         h1 = findobj(gca,'Type','patch');
% % % %         set(h1,'facealpha',0.75);
% % % %         title([ALL_CATEGORIES{param_idx} ' mode0=' num2str(M0) '. '  ' mode1= ' num2str(M1) '. isnan0= ' num2str(isnan0)  '. isnan1= ' num2str(isnan1)]);
% % % 
% % %         scatter(mean_1,ones(size(mean_1)),'xr')
% % %         mean1=mean(mean_1(~isnan(mean_1)));
% % %        % vari=var(mean_1(~isnan(mean_1)));
% % %         medi=median(mean_1(~isnan(mean_1)));
% % %         line([mean1 mean1],[1 0.5],'Color','r')
% % %         
% % %         line([medi medi],[1 0.5],'Color','r','LineStyle','--')
% % %         
% % %         hold on
% % %         scatter(mean_0,zeros(size(mean_0)),'ob')
% % %         mean0=mean(mean_0(~isnan(mean_0)));
% % %         %vari=var(mean_0(~isnan(mean_0)));
% % %          medi=median(mean_0(~isnan(mean_0)));
% % %         line([mean0 mean0],[0 0.5],'Color','b')
% % %         line([medi medi],[0 0.5],'Color','r','LineStyle','--')
% % %       %  title([ALL_CATEGORIES{param_idx} ' mode0=' num2str(M0) '. '  ' mode1= ' num2str(M1) '. isnan0= ' num2str(isnan0)  '. isnan1= ' num2str(isnan1)]);
% % % 
% % %                 title([ALL_CATEGORIES{param_idx} ' norm=' num2str(max_score1_norm) '. '  ' abs= ' num2str(max_score1_abs)  ]);
% % % 
% % %         
% % %         subplot(1,2,2)
% % %         
% % %         scatter(mean_1_abs,ones(size(mean_1_abs)),'xr')
% % %         mean1_abs=mean(mean_1_abs(~isnan(mean_1_abs)));
% % %        % vari=var(mean_1(~isnan(mean_1)));
% % %         medi_abs=median(mean_1_abs(~isnan(mean_1_abs)));
% % %         line([mean1_abs mean1_abs],[1 0.5],'Color','r')
% % %         
% % %         line([medi_abs medi_abs],[1 0.5],'Color','r','LineStyle','--')
% % %         
% % %         hold on
% % %         scatter(mean_0_abs,zeros(size(mean_0_abs)),'ob')
% % %         mean0_abs=mean(mean_0_abs(~isnan(mean_0_abs)));
% % %         %vari=var(mean_0(~isnan(mean_0)));
% % %          medi_abs=median(mean_0_abs(~isnan(mean_0_abs)));
% % %         line([mean0_abs mean0_abs],[0 0.5],'Color','b')
% % %         line([medi_abs medi_abs],[0 0.5],'Color','r','LineStyle','--')
% % %       %  title([ALL_CATEGORIES{param_idx} ' mode0=' num2str(M0) '. '  ' mode1= ' num2str(M1) '. isnan0= ' num2str(isnan0)  '. isnan1= ' num2str(isnan1)]);
% % % 
% % %                 title([ALL_CATEGORIES{param_idx} ' norm=' num2str(max_score1_norm) '. '  ' abs= ' num2str(max_score1_abs)  ]);
% % % 
% % %         
% % % 
% % %         pause(.1)
     
        
        
        
        
        
%         data = rand(1, 100);
% [N, X] = hist(data, 0:0.1:1);
% 
% subplot(2, 1, 1);
% hist(data, 0:0.1:1);
% 
% subplot(2, 1, 2);
% bar(X, N./sum(N), 1);
        
    end

    function best_th= th_search(X_Test,IHD_Test)

        DATA=zeros(size(X_Test,1),3);
        
        best_lam=0;
        max_score1_lam=0;
        max_score2_lam=0;
        best_lam_score=[];
        lam=1;



        for lam_idx=1:length(lam)

            
            
            B = mnrfit(X_Test,IHD_Test+1);
            PHAT = mnrval(B,X_Test);
            
            
            if size(PHAT,2)==1
                max_score1_lam=-1;
                best_th=-1;
                return
            end
            
%         Y=(IHD-1)+IHD; % classes y ={-1,1}
%         X=X_Test;
%         R=~isnan(X);
%         PHAT=teach_log_reg(X,Y,R,lam(lam_idx));

        %[IDX,C,sumd] = kmeans(double(R),5);

        %X_trunc_saps=[SAPS_SCORES(isnan(PHAT(:,2)))-SAPS_SCORES_48(isnan(PHAT(:,2))) SAPS_SCORES(isnan(PHAT(:,2))) ];
        %B_saps = mnrfit(X_trunc_saps,IHD(isnan(PHAT(:,2)))+1);
        %PHAT_saps = mnrval(B_saps,X_trunc_saps);

        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(isnan(PHAT(:,2)),2)=PHAT_saps(:,2);
            DATA(:,3)=PHAT(:,2)> class_th;
           % DATA(isnan(PHAT(:,2)),3)=PHAT_saps(:,2)> 0.33;

        % % %     DATA(:,2)=PHAT(:,2);
        % % %     DATA(isnan(PHAT(:,2)),2)=PHAT_saps(:,2);
        % % %     DATA(~isnan(PHAT(:,2)),3)=PHAT(~isnan(PHAT(:,2)),2)> class_th;
        % % %     DATA(isnan(PHAT(:,2)),3)=PHAT_saps(:,2)> 0.33;

            DATA(DATA(:,2)<0.01,2)=0.01;
            DATA(DATA(:,2)>0.99,2)=0.99;
            %DATA(:,3)=PHAT1(:,2)>0.17 & PHAT2(:,2)>class_th ;

        %     if max_score1==0
        %           savefile = 'DATA_learn.mat';
        %         save(savefile, 'DATA')
        %    
        %     end


       %  if(~isempty(results))


            % Calculate sensitivity (Se) and positive predictivity (PPV)
            TP=sum(DATA(IHD_Test==1,3));
            FN=sum(~DATA(IHD_Test==1,3));
            FP=sum(DATA(IHD_Test==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_Test 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

            % The event 2 score is the Hosmer-Lemeshow statistic (H).
            %display(['Unofficial Event 2 score: ' num2str(H)]);

%             figure(1000)
%             hold on
%             scatter(class_th,score1,'x')
       %  end


        %B'

        end
        best_lam_score(lam_idx)=max_score1;
        if max_score1 > max_score1_lam
         best_lam=0;
         max_score1_lam=max_score1;
         max_score2_lam=max_score2;
%         best_th
        end

        end
    end


end