function [dens0,dens,plage,C,S]=estimcatSJgr6(X,nbplage)

%construit une plage
X=sort(X);
D=size(X);
N=D(1);
mini=min(X);
maxi=max(X);
taille=maxi-mini;
mm=mini-taille/5;
MM=maxi+taille/5;
plage=[mm:(MM-mm)/nbplage:MM]';
pas=(MM-mm)/nbplage;
Dp=size(plage);
Np=Dp(1);
ZZ=pdist(X);
Md=squareform(ZZ);
Dist2=Md.^2;
%[logp, s] = kde1(X, plage);
s=bandwidth_SJ(X,'norm');
s=real(s);
dens=zeros(Np,1);
for i=1:Np
    for j=1:N
        u=(plage(i)-X(j))/s;
        dens(i)=dens(i)+(1/s)*normpdf(u,0,1);
    end
end
dens=dens/N;

dens0=dens;
dansa=dens0;
C=ones(N,1);
H=zeros(1,1);
H(1,1)=s;
LCV0=LSCVgauss(Dist2,H,C);
figure
plot(plage,dens)
hold on
pause
Ca=C;
densa=dens;
nbcoupe=1
nbcoupe=input('how many bandwidthes ?')
nbcoupe=nbcoupe-1;
Sa=s;
LCVa=LCV0;
if nbcoupe>0
    coup=0;
    testcont=1
    while testcont==1
        coup=coup+1;

        %recherche les minimum local de la densité
        minim=zeros(Np,1);
        indicem=zeros(Np,1);
        c=0;
        for i=2:Np-1
            if dens(i)<=dens(i+1)
                if dens(i)<=dens(i-1)
                    c=c+1;
                    minim(c)=plage(i);
                    indicem(c)=i;
                end
            end
        end
        minimum=minim(1:c,1);
        indicemin=indicem(1:c,1);
        %regroupe pour éviter trop de minimum
        cl=ones(c,1);
        for i=2:c
            if indicemin(i)==indicemin(i-1)+1
                cl(i)=cl(i-1);
            else
                cl(i)=cl(i-1)+1;
            end
        end
        %fait des min par classes
        nb=max(cl);
        mini=zeros(nb,1);
        efm=zeros(nb,1);
        for i=1:c
            mini(cl(i))=mini(cl(i))+minim(i);
            efm(cl(i))=efm(cl(i))+1;
        end
        mini=mini./efm;
        %teste toute les coupures
        LCV=zeros(nb,1);
        for ii=1:nb
            seuil=mini(ii);
            C=ones(N,1);
            seuilsmM=ones(2,2);
            seuilsmM(:,2)=N*ones(2,1);
            for i=1:N
                if X(i)>seuil
                    C(i)=2;
                end
            end
            for j=1:2
                for i=2:N-1
                    if C(i)==j
                        if C(i+1)==j+1
                            seuilsmM(j,2)=i;
                            seuilsmM(j+1,1)=i+1;
                        end
                    end
                end
            end
            S=zeros(2,1);
            for i=1:2;
                sh=bandwidth_SJ(X(seuilsmM(i,1):seuilsmM(i,2)),'norm');
                S(i)=real(sh);
            end
            LCV(ii)=LSCVgauss(Dist2,S,C);
        end

        %recupere les (au max nbcoupe) coupures qui diminue le LCV
        [LCVt,It]=sort(LCV);
        nblcvinf=0;
        for i=1:nb
            if LCVt(i)<LCV0
                nblcvinf=nblcvinf+1;
            end
        end
        sicoupe=zeros(nb,1);

        if nblcvinf>=1
            nbcoupure=min(nbcoupe,nblcvinf);
            indicecoupfin=zeros(nbcoupure,1);
            for i=1:nbcoupure
                sicoupe(It(i))=1;
            end
        end
        %fait une estimation de densité finale avec tout le monde
        NbCl=sum(sicoupe);
        cc=0;
        seuils=zeros(NbCl,1);
        for i=1:nb
            if sicoupe(i)==1
                cc=cc+1;
                seuils(cc)=mini(i);
            end
        end
        C=ones(N,1);
        seuilsmM=ones(NbCl+1,2);
        seuilsmM(:,2)=N*ones(NbCl+1,1);
        for i=1:N
            for j=1:NbCl
                if X(i)>seuils(j)
                    C(i)=j+1;
                end
            end
        end
        for j=1:NbCl+1
            for i=2:N-1
                if C(i)==j
                    if C(i+1)==j+1
                        seuilsmM(j,2)=i;
                        seuilsmM(j+1,1)=i+1;
                    end
                end
            end
        end
        S=zeros(NbCl+1,1);
        for i=1:NbCl+1
            sh=bandwidth_SJ(X(seuilsmM(i,1):seuilsmM(i,2)),'norm');
            S(i)=real(sh);
        end
        LCVap=LSCVgauss(Dist2,S,C);
        LCVap
        LCVa
        if LCVap<LCVa
            dens=zeros(Np,1);
            for i=1:Np
                for j=1:N
                    u=(plage(i)-X(j))/S(C(j));
                    dens(i)=dens(i)+(1/S(C(j)))*normpdf(u,0,1);
                end
            end
            dens=dens/N;
            if sum(abs(Ca-C))==0
                testcont=0;
            end
            %plot(plage,dens)
            %hold on
            %pause
            densa=dens;
            Ca=C;
            Sa=S;
            LCVa=LCVap;
            if coup==10
                testcont=0;
            end
        else
            testcont=0;
            dens=densa;
            C=Ca;
            S=Sa;
        end 
    end
end