%macro lissage(index=);

proc iml;
    use &index;
    read all var _all_ into index;

        n=nrow(index);
        cvb1=j(1000,1,0);
        liss=j(1000,1,0);
        plot=j(1000,2,0);

            do l=1 to 1000 ;

                b=l/1000;
                c1=1/(n*b);
                c2=1/n;
                c3=2/(n-1);
                somme1=0;
                somme2=0;

                    do i=1 to n;

                        do j=1 to n;
                            som1=1/(2*sqrt(3.14))*exp(-1/4*((index[j]-index[i])/b)##2);
                            somme1=somme1+som1;
                            som2=1/sqrt(2*3.14)*exp(-1/2*((index[j]-index[i])/b)##2);
                            somme2=somme2+som2;

                        end;

                    end;

                cvb=c1*((c2*somme1)-(c3*somme2)+(c3*(n/sqrt(2*3.14))));
                cvb1[l,1]=cvb;
                liss[l,1]=l/1000;
                plot[l,1]=cvb;
                plot[l,2]=l/1000;

            end;
        
        resulta=min(cvb1);

        fin=j(1000,2,0);

            do j=1 to 1000;

                if cvb1[j,1]=resulta then fin[j,1]=cvb1[j,1];
                else fin[j,1]=0;
                if cvb1[j,1]=resulta then fin[j,2]=liss[j,1];
                else fin[j,2]=0;

            end;

    create sortie FROM fin ;
    append FROM fin ;
    close sortie ;
    create graph FROM plot ;
    append FROM plot ;
    close graph ;

quit;

data sortie; set sortie;
    label COL1='valeur de la lscv';
    label COL2='paramétre de lissage';
run;

proc print data=sortie label noobs;
    where COL1^=0;
run;

proc sort data=graph; 
    by COL2;
run;

symbol i=join v=none c=blue;
axis1 label=('h' justify=right);
axis2 label=('lscv' justify=right);
proc gplot data=graph;
    plot COL1*COL2 / vref=0 haxis=axis1 vaxis=axis2;
    title 'La fonction lscv';
run;

%mend;

%lissage(index=score);