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