OPTIONS LS=75 PS=60 NODATE;
DATA DIAG;
INPUT n p ajj; h2=0.25; r=.5;
di=( ( ( 1+(n-1)*r) /n ) + (p -1) * ajj * h2 ) / p;
*INPUT n p ajj rel $; CARDS; 2 20 0.25 1 30 0.25 1 135 0.25 3 1 1.00 ; PROC PRINT; RUN;
DATA OD; INPUT G1 G2 G3 G4; CARDS; 0 0.25 0.125 0.25 0.25 0 0.125 0.25 0.125 0.125 0 0 0.25 0.25 0 0 ;
PROC IML; START;
h2=0.25; P={17852 18751 18000 20000}; PMEAN={14000 14000 14000 14000};
P=P-PMEAN;
aij={0.5, 0.5, 0.25, 0.5};
RHS=aij*h2;
USE DIAG; READ all into DII; print DII; di = DII (|,NCOL(DII) |); DI = DIAG(di`); PRINT di;
USE OD; READ all into OD; OD=OD*h2; PRINT OD; VP=OD+DI; PRINT VP; b=INV(VP)*RHS; PRINT b;
INDEX=P*b; PRINT INDEX;
r_IH= SQRT(b`*aij);PRINT r_IH;
FINISH; RUN;