Dear Rick, Thanks for contributing such a wonderful book on SAS IML. I also bout the book and expecting it within coming days. I need your support to sove the following errors from SAS IML Studio. ========================================================================= »ERROR: NLPNRR call: Module GRAD must compute a 6 element row vector. »ERROR: Execution error as noted previously. (rc=100) » » operation : NLPNRR at line 212 column 1 » operands : *LIT1085, beta0, optn, , , , , *LIT1086, *LIT1087 » »*LIT1085 1 row 1 col (character, size 2) » » LL » »beta0 6 rows 1 col (numeric) » » 0 » 0 » 0 » 0 » 0 » 0 » »optn 1 row 2 cols (numeric) » » 0 2 » »*LIT1086 1 row 1 col (character, size 4) » » GRAD » »*LIT1087 1 row 1 col (character, size 4) » » HESS » » statement : CALL at line 212 column 1 ============================================================================== Please find below my SAS IML programs: ============================================================================== SUBMIT; %macro IPC(dat=_last_, y=, x=, beta=0); PROC IML; RESET noname; use &dat; read all var {&x} into X; read all var {&y} into Y; /* Design matrix with intercept */ X = j(nrow(X),1,1)||X; P = NCOL(X); START FMat(Y); N1 = NROW(Y); G1 = NCOL(UNIQUE(Y)); F = j(N1, G1, 0); DO i = 1 TO N1; g1 = Y; F[i, g1] = 1; END; RETURN(F); FINISH FMat; /* Caling transformaion module */ F = FMat(Y); N = NROW(F); G = NCOL(F); D = G - 1; gamma = {0 0 0, 1 0 0, 0 1 0, 0 0 1}; * Define log likelihood function; START LL(beta) global(F, X, gamma, N, P, G, D); loglik = 0; eta = j(N, D, 0); /* Ideal points */ DO i=1 TO N; DO j=1 TO D; tmp = 0; DO k=1 TO P; tmpjk = X[i, k] * beta[(j-1)*P + k]; tmp = tmp + tmpjk; END; eta[i, j] = tmp; END; END; /* Unfolding/Eucliean Distance */ delta = vecdiag(eta * t(eta)) * j(1, nrow(gamma), 1) + j(nrow(eta), 1, 1) * t(vecdiag(gamma * t(gamma))) - 2 * eta * t(gamma); /* Category probability */ pi = j(N, G, 0); DO i = 1 to N; DO j=1 TO G; pi[i,j] = exp(-0.5 * delta[i,j]) / exp(-0.5 * delta[i, +]); END; END; DO i = 1 TO N; DO j = 1 TO G; ijloglik = F[i, j] * log(pi[i, j]); loglik = loglik + ijloglik; END; END; RETURN(loglik); FINISH LL; * Define gradient; START GRAD(beta) global(F, X, gamma, N, P, G, D); dl = j(P*D, 1, 0); eta = j(N, D, 0); /* Ideal points */ DO i=1 TO N; DO j=1 TO D; tmp = 0; DO k=1 TO P; tmpjk = X[i, k] * beta[(j-1)*P + k]; tmp = tmp + tmpjk; END; eta[i, j] = tmp; END; END; /* Unfolding/Eucliean Distance */ delta = vecdiag(eta * t(eta)) * j(1, nrow(gamma), 1) + j(nrow(eta), 1, 1) * t(vecdiag(gamma * t(gamma))) - 2 * eta * t(gamma); /* Category probability */ pi = j(N, G, 0); DO i = 1 to N; DO j=1 TO G; pi[i,j] = exp(-0.5 * delta[i,j]) / exp(-0.5 * delta[i, +]); END; END; DO j = 1 to D; DO k = 1 TO P; dl[(j-1)*P + k, 1] = (F[, j] - pi[, j])` * X[, k]; END; END; RETURN(dl); FINISH GRAD; * Define Hessian matrix; START HESS(beta) global(F, X, gamma, N, P, G, D); d2l = j(P*D, P*D, 0); eta = j(N, D, 0); /* Ideal points */ DO i=1 TO N; DO j=1 TO D; tmp = 0; DO k=1 TO P; tmpjk = X[i, k] * beta[(j-1)*P + k]; tmp = tmp + tmpjk; END; eta[i, j] = tmp; END; END; /* Unfolding/Eucliean Distance */ delta = vecdiag(eta * t(eta)) * j(1, nrow(gamma), 1) + j(nrow(eta), 1, 1) * t(vecdiag(gamma * t(gamma))) - 2 * eta * t(gamma); /* Category probability */ pi = j(N, G, 0); DO i = 1 to N; DO j=1 TO G; pi[i,j] = exp(-0.5 * delta[i,j]) / exp(-0.5 * delta[i, +]); END; END; jj = 1; /* jj maintains a running counter of the current row of beta_tmp and xtwx_tmp */ DO j = 1 TO D; w1 = pi[ ,j] * (1 - pi[ ,j])`; DO k = 1 TO P; kk = jj - 1; /* kk maintains a running counter of the current column index xtwx_tmp */ /* second derivative when jprime=j */ DO kprime = k TO P; kk = kk + 1; d2tempj = d2l[jj, kk]; d2l[jj, kk] = d2tempj + (X[ ,k])` * w1 * X[ ,kprime]; d2l[kk,jj] = d2l[jj,kk]; END; /* second derivative when jprime^=j */ DO jprime = j+1 TO D; w2 = - pi[ ,j] * (pi[ ,jprime])`; DO kprime = 1 TO P; kk = kk + 1; d2tempdprime = d2l[jj, kk]; d2l[jj, kk] = d2tempdprime + (X[ ,k])` * w2 * X[ ,kprime]; d2l[kk,jj] = d2l[jj,kk]; END; END; jj = jj + 1; END; END; RETURN(d2l); FINISH HESS; optn = {0 2}; beta0 = β if beta0=0 then beta0 = j(P*D, 1, 0); else print 'User-supplied starting values for beta = ' beta0; CALL NLPNRR(rc, beta, 'LL', beta0, optn) grd='GRAD' hes='HESS'; CALL NLPFDD(f, g, h, 'LL', beta); var = inv(-h); sd = sqrt(vecdiag(var)); print 'Hessian = ' h ' Covariance matrix = ' var 'Standard errors = ' sd; run; %mend; %IPC(dat=nesda2, y=resp, x=sex); ENDSUBMIT; ============================================================== Thnk you very much for your help. Best, Haile NL
... View more
