function [RAD,CONV,PCRET,PBRET,VARRET,EM,FIG] = exp_max(AS,CK,QUAL,S,PLOT,ASR)

if nargin < 5 | isempty(PLOT)
   PLOT = 0;
end

E=size(AS,1);
D=E-2;
R=sqrt(E-1);

RD=dist_misval(AS,CK,2);
RD=RD(find(isfinite(RD)));
samples = size(RD,2);
if nargin==6 & ~isempty(ASR)
    RDR=dist_misval(ASR,CK,2);
    RDR=RDR(find(isfinite(RD)));
end

MAXITER=200;

CDIF=0.001;

PC=sum(RD<QUAL)/size(RD,2);
PB=1-PC;
VAR=(1/D)*sum(RD(find(RD<QUAL)).^2)/sum(RD<QUAL);

CONV=0;
FIG=0;
for em=1:MAXITER
   prc=clusterdistrib(RD,VAR,D,R);
   prb=background(RD,D,R);
   prcpc=prc*PC;
   prbpb=prb*PB;
   pr=sum([prcpc;prbpb],1);
   pcr=prcpc./pr;
   SM=sum(pcr);
   if SM==0 | ~isfinite(SM)
       break
   end
   VAR_new=(1/D)*sum((RD.^2).*pcr)/SM;
   PC_new=SM/samples;
   PB_new=1-PC_new;
   if abs(VAR_new-VAR)<CDIF & abs(PC_new-PC)<CDIF
      CONV=1;
      break
   end
   PC=PC_new;
   PB=PB_new;
   VAR=VAR_new;     
end

if CONV==1
   PC=PC_new;
   PB=PB_new;
   VAR=VAR_new;
   
   if PC==0 | PB==0 
      CONV=0;
      disp(['             EM: No or incorrect convergence'])
      RAD=0;
      VARRET=0;
      PCRET=0;
      PBRET=0;
      EM=0;
      return
   end
   
   SD=(2*(pi^(D/2)))/GAMMA(D/2);
   SD1=(2*(pi^((D+1)/2)))/GAMMA((D+1)/2);
   CC=SD*(1/((2*pi*VAR)^(D/2)));
   CB=(SD/(SD1*(R^D)));
   LO=(S/(1-S))*((PB*CB)/(PC*CC));
   if LO <= 0
      CONV=0;
      disp(['             EM: Not possible to calculate radius'])
      RAD=0;
      VARRET=0;
      PCRET=0;
      PBRET=0;
      EM=0;
      return
   end
   DIS=-2*VAR*log(LO);

   if DIS < 0
      CONV=0;
      disp(['             EM: Not possible to calculate radius'])
      RAD=0;
      VARRET=0;
      PCRET=0;
      PBRET=0;
      EM=0;
      return
   end
   
   RAD=sqrt(DIS);
   VARRET=VAR;
   PCRET=PC;
   PBRET=PB;
   EM=em;
   
   if PLOT==1
       SZ=size(RD,2);
       FIG=figure;
       if nargin==6 & ~isempty(ASR)
           subplot(1,2,1)
       end
       hist(RD,100);
       [NH XH]=hist(RD,100);
       DELTA=(max(XH)-min(XH))/99;
       hold on
       INT=DELTA:DELTA:R;
       prc=clusterdistrib(INT,VAR,D,R);
       prb=background(INT,D,R);
       plot(INT,prc*PC*SZ*DELTA)
       plot(INT,prb*PB*SZ*DELTA,'r')    
       plot(INT,(prc*PC+prb*PB)*SZ*DELTA,'g-')
       
       if nargin==6 & ~isempty(ASR)
           subplot(1,2,2)
           hist(RDR,100)
           [NHR XHR]=hist(RDR,100);
           DELTAR=(max(XHR)-min(XHR))/99;
           hold on
           INTR=DELTAR:DELTAR:R;
           prbr=background(INTR,D,R);
           plot(INTR,prbr*SZ*DELTAR,'r')    
       end
   end
      
else
   disp(['             EM: No or incorrect convergence'])
   RAD=0;
   VARRET=0;
   PCRET=0;
   PBRET=0;
   EM=0;
end


   