gaussmix.m

% gaussmix  8/27/20

clear
format compact

load GaussMix

N1 = 12;
N2 = 11;
N = N1 + N2;

[sy,sx] = size(B);

x = -1.4:0.01:1.;
xi0 = where(x,0);

fun1 = zeros(size(x));
for loop = 1:N1
    
    m = B(loop,1);
    s = B(loop,3);
    
    G = gaussprob(x,m,s);
    
    fun1 = fun1 + G;
end


fun2 = zeros(size(x));
for loop = N1+1:N1+N2
    
    m = B(loop,1);
    s = B(loop,3);
    
    G = gaussprob(x,m,s);
    
    fun2 = fun2 + G;
end


figure(1)
plot(x,fun1/N,x,fun2/N)

lim = length(x);
for loop = 1:lim
    
    i1(loop) = trapz(fun1(1:loop))*0.01/N1;
    i2(loop) = trapz(fun2(1:loop))*0.01/N2;
        
end

figure(2)
plot(x,i1,x,i2)


% Find optimal decision point
arg = sqrt((i2.^2) + (1-i1).^2);
[X,I] = min(arg);

ind = I;

TN = N1*i1(ind);
FN = N2*i2(ind);
TP = N2*(1 - i2(ind));
FP = N1*(1 - i1(ind));

Sens = TP/(TP + FN);
Spec = TN/(TN + FP);
PLR = Sens/(1 - Spec);
NLR = (1 - Sens)/Spec;
PPV = TP/(TP + FP);
NPV = TN/(TN + FN);

Acc = (TP + TN)/N;

disp(' ');
disp('ROC optimum (red)')
displine('xthresh = ',x(I))
displine('Acc = ',Acc);
displine('Sens = ',Sens);
displine('Spec = ',Spec);
displine('PLR = ',PLR);
displine('NLR = ',NLR);
displine('PPV = ',PPV);
displine('NPV = ',NPV);


% Zero-threshold
i10 = trapz(fun1(1:xi0))*0.01/N1;
i20 = trapz(fun2(1:xi0))*0.01/N2;

ind = xi0;

TN = N1*i1(ind);
FN = N2*i2(ind);
TP = N2*(1 - i2(ind));
FP = N1*(1 - i1(ind));

Sens = TP/(TP + FN);
Spec = TN/(TN + FP);
PLR = Sens/(1 - Spec);
NLR = (1 - Sens)/Spec;
PPV = TP/(TP + FP);
NPV = TN/(TN + FN);

Acc = (TP + TN)/N;

disp(' ');
disp('Zero threshold (blue)')
displine('Acc = ',Acc);
displine('Sens = ',Sens);
displine('Spec = ',Spec);
displine('PLR = ',PLR);
displine('NLR = ',NLR);
displine('PPV = ',PPV);
displine('NPV = ',NPV);

% Best separation
ixmin = 121;

ind = ixmin;

TN = N1*i1(ind);
FN = N2*i2(ind);
TP = N2*(1 - i2(ind));
FP = N1*(1 - i1(ind));

Sens = TP/(TP + FN);
Spec = TN/(TN + FP);
PLR = Sens/(1 - Spec);
NLR = (1 - Sens)/Spec;
PPV = TP/(TP + FP);
NPV = TN/(TN + FN);

Acc = (TP + TN)/N;

disp(' ');
disp('Min threshold (green)')
displine('Acc = ',Acc);
displine('Sens = ',Sens);
displine('Spec = ',Spec);
displine('PLR = ',PLR);
displine('NLR = ',NLR);
displine('PPV = ',PPV);
displine('NPV = ',NPV);


dx = diff(i2);
AUC = sum(i1(1:lim-1).*dx);
disp(' ');
displine('AUC = ',AUC);



figure(3)
plot(i2,i1,'r','LineWidth',2)
hold on
plot(i2(ixmin),i1(ixmin),'go','MarkerFaceColor','g')
plot(i2(I),i1(I),'ro','MarkerFaceColor','r')
plot(i20,i10,'bo','MarkerFaceColor','b')
hold off



disp(' ')
Printfile5('hovgauss.txt',x,fun1/N,fun2/N,i1,i2);


stat = class2stat(B(1:N1,1),B(N1+1:N,1))

% lowthresh = round(stat.meanx + stat.cimin/2,2);
% hithresh = round(stat.meanx + stat.cimax/2,2);
% 
% lt = where(x,lowthresh);
% ht = where(x,hithresh);

del = stat.mn - stat.cimin;
lowt = stat.meanx + stat.mn/2 - del;
hit = stat.meanx + stat.mn/2 + del;

lowthresh = round(lowt,2);
hithresh = round(hit,2);

lt = where(x,lowthresh);
ht = where(x,hithresh);

% Low confidence limit
ixmin = lt;

ind = ixmin;

TN = N1*i1(ind);
FN = N2*i2(ind);
TP = N2*(1 - i2(ind));
FP = N1*(1 - i1(ind));

Sens = TP/(TP + FN);
Spec = TN/(TN + FP);
PLR = Sens/(1 - Spec);
NLR = (1 - Sens)/Spec;
PPV = TP/(TP + FP);
NPV = TN/(TN + FN);

Acc = (TP + TN)/N;

disp(' ');
disp('Min CI')
displine('Acc = ',Acc);
displine('Sens = ',Sens);
displine('Spec = ',Spec);
displine('PLR = ',PLR);
displine('NLR = ',NLR);
displine('PPV = ',PPV);
displine('NPV = ',NPV);


% High confidence limit
ixmin = ht;

ind = ixmin;

TN = N1*i1(ind);
FN = N2*i2(ind);
TP = N2*(1 - i2(ind));
FP = N1*(1 - i1(ind));

Sens = TP/(TP + FN);
Spec = TN/(TN + FP);
PLR = Sens/(1 - Spec);
NLR = (1 - Sens)/Spec;
PPV = TP/(TP + FP);
NPV = TN/(TN + FN);

Acc = (TP + TN)/N;

disp(' ');
disp('Max CI')
displine('Acc = ',Acc);
displine('Sens = ',Sens);
displine('Spec = ',Spec);
displine('PLR = ',PLR);
displine('NLR = ',NLR);
displine('PPV = ',PPV);
displine('NPV = ',NPV);