/* rettelse 11/12/91                                     USIJB  */
/* linie 38    Do Q = 1 to Cp+1;  rettet til   Do Q = 1 to Cp;  */
option nosource;
/***************************************************************
 ***                                                         ***
 ***   Program til beregning af den multiple Cooks afstand   ***
 ***                                                         ***
 ***   BEM#RK ||| P} grund af afruningsfejl i maskinen n}r   ***
 ***   Binomial koeficienten beregnes, kan programmet slutte ***
 ***   Med fejlmeddelelsen INDEXING PROBLEM. Dette er uden   ***
 ***   betydning, da alle v{rdier er beregnet.               ***
 ***                                                         ***
 ***   Programmet afvikles med programmet EXECCOOK           ***
 ***                                                         ***
 ***************************************************************/;

%Macro Makcook(Data=INDDATA, Out=UDDATA,Y=Respons, X=X1 X2 X3, M=3
outliere, sort=n, print = all);

Title3 "Data : &data" "   Model : &Y =F(&X)";

Proc Iml;

Start Binom(Cp,N,M);
   CpT = Sum(log(1:N));
   CpN = sum(log(1:M))+Sum(log(1:(N-M)));
   Cp  = Exp(Cpt-Cpn);
Finish;

Start COOKD(Y,X,P,N,M);
Run Binom(Cp,N,M);
XX = X`*X;
 B = INV(XX)*X`*Y;
 S = Ssq(Y-X*B)/(N-P);
Ps = P#S;
V  = J(M,1,"I");
Vname = concat(V,Char((1:M)`,1,0))//{"COOKD"};
Do Q = 1 to Cp;
   IF Q = 1 then do;
      I = 1 : M;
      goto BEREGN;
   End;

   IF I(|1,M|) = N Then Goto S1;
      I(|1,M|) = I(|1,M|) + 1;
      GOTO BEREGN;

S1:
   IF I(|1,M-1|) = N -1 Then Goto S2;
      I(|1,M-1|) = I(|1,M-1|) + 1;
      I(|1,M|)   = I(|1,M-1|) + 1;
      GOTO BEREGN;

S2:
   IF I(|1,M-2|) = N -2 Then Goto S3;
      I(|1,M-2|) = I(|1,M-2|) + 1;
      I(|1,M-1|) = I(|1,M-2|) + 1;
      I(|1,M|)   = I(|1,M-1|) + 1;
      GOTO BEREGN;

S3:
   IF I(|1,M-3|) = N-3 Then Goto S4;
      I(|1,M-3|) = I(|1,M-3|) + 1;
      I(|1,M-2|) = I(|1,M-3|) + 1;
      I(|1,M-1|) = I(|1,M-2|) + 1;
      I(|1,M|)   = I(|1,M-1|) + 1;
      GOTO BEREGN;

S4:
   IF I(|1,M-4|) = N-4 Then Goto S5;
      I(|1,M-4|) = I(|1,M-4|) + 1;
      I(|1,M-3|) = I(|1,M-4|) + 1;
      I(|1,M-2|) = I(|1,M-3|) + 1;
      I(|1,M-1|) = I(|1,M-2|) + 1;
      I(|1,M|)   = I(|1,M-1|) + 1;
      GOTO BEREGN;


S5:
   IF I(|1,M-5|) = N-5 Then Goto S6;
      I(|1,M-5|) = I(|1,M-5|) + 1;
      I(|1,M-4|) = I(|1,M-5|) + 1;
      I(|1,M-3|) = I(|1,M-4|) + 1;
      I(|1,M-2|) = I(|1,M-3|) + 1;
      I(|1,M-1|) = I(|1,M-2|) + 1;
      I(|1,M|)   = I(|1,M-1|) + 1;
      GOTO BEREGN;

S6:
   IF I(|1,M-6|) = N-6 Then Goto BEREGN;
      I(|1,M-6|) = I(|1,M-6|) + 1;
      I(|1,M-5|) = I(|1,M-6|) + 1;
      I(|1,M-4|) = I(|1,M-5|) + 1;
      I(|1,M-3|) = I(|1,M-4|) + 1;
      I(|1,M-2|) = I(|1,M-3|) + 1;
      I(|1,M-1|) = I(|1,M-2|) + 1;
      I(|1,M|)   = I(|1,M-1|) + 1;
      GOTO BEREGN;


BEREGN:
    V  = J(N,1,0);
    V(|I,|)=1;
    Xi = X(|Loc(V=0),|);
    Yi = Y(|Loc(V=0),|);
    Bi = (Inv(Xi`*Xi))*Xi`*Yi;
    CD = (B-Bi)`*XX*(B-Bi)/ps;
    UD = I||CD;
    If Q = 1 then do;
       create &Out from UD (|Colname = Vname|);
       append from UD;
    end;
    IF q > 1 then do;
    append from ud;
    end;
    End;
finish;

Reset noname;

Res  = { &Y };
Fork = { &X };

* Indl{s data;
Use &Data;
Read  All Var Res Into Y;
Read  All Var Fork Into X;

* Beregn Hj{lpevariable;
M = { &M };
N = Nrow(X);
X = J(N,1,1)||X;
P = Ncol(X);
Run COOKD(Y,X,P,N,M);
quit;

%if &sort=y %then %do;
    proc sort data=&out;
    by descending cookd;
%end;
%let pri=%upcase(&print);
%If %str(&pri) ne %str(ALL) %then %do;
    Data &Out;
    Set &Out;
    IF _N_ > &Print then delete;
%end;
Proc print Data=&out;
run;
%Mend;